pml_helmholtz_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for Helmholtz elements
27#ifndef OOMPH_PML_HELMHOLTZ_ELEMENTS_HEADER
28#define OOMPH_PML_HELMHOLTZ_ELEMENTS_HEADER
29
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36#include "math.h"
37#include <complex>
38
39// OOMPH-LIB headers
40#include "../generic/nodes.h"
41#include "../generic/Qelements.h"
42#include "../generic/oomph_utilities.h"
43#include "../generic/pml_meshes.h"
44#include "../generic/projection.h"
45#include "../generic/pml_mapping_functions.h"
46
47/// /////////////////////////////////////////////////////////////////
48/// /////////////////////////////////////////////////////////////////
49/// /////////////////////////////////////////////////////////////////
50
51namespace oomph
52{
53 //=============================================================
54 /// A class for all isoparametric elements that solve the
55 /// Helmholtz equations with pml capabilities.
56 /// This contains the generic maths. Shape functions, geometric
57 /// mapping etc. must get implemented in derived class.
58 //=============================================================
59 template<unsigned DIM>
60 class PMLHelmholtzEquations : public virtual PMLElementBase<DIM>,
61 public virtual FiniteElement
62 {
63 public:
64 /// Function pointer to source function fct(x,f(x)) --
65 /// x is a Vector!
66 typedef void (*PMLHelmholtzSourceFctPt)(const Vector<double>& x,
67 std::complex<double>& f);
68
69 /// Constructor
71 {
72 // Provide pointer to static method (Save memory)
75 }
76
77
78 /// Broken copy constructor
80
81 /// Broken assignment operator
82 // Commented out broken assignment operator because this can lead to a
83 // conflict warning when used in the virtual inheritence hierarchy.
84 // Essentially the compiler doesn't realise that two separate
85 // implementations of the broken function are the same and so, quite
86 // rightly, it shouts.
87 /*void operator=(const PMLHelmholtzEquations&) = delete;*/
88
89 /// Return the index at which the unknown value
90 /// is stored.
91 virtual inline std::complex<unsigned> u_index_helmholtz() const
92 {
93 return std::complex<unsigned>(0, 1);
94 }
95
96 /// Get pointer to k_squared
97 double*& k_squared_pt()
98 {
99 return K_squared_pt;
100 }
101
102
103 /// Get the square of wavenumber
104 double k_squared()
105 {
106#ifdef PARANOID
107 if (K_squared_pt == 0)
108 {
109 throw OomphLibError(
110 "Please set pointer to k_squared using access fct to pointer!",
111 OOMPH_CURRENT_FUNCTION,
112 OOMPH_EXCEPTION_LOCATION);
113 }
114#endif
115 return *K_squared_pt;
116 }
117
118 /// Alpha, wavenumber complex shift
119 const double& alpha() const
120 {
121 return *Alpha_pt;
122 }
123
124 /// Pointer to Alpha, wavenumber complex shift
125 double*& alpha_pt()
126 {
127 return Alpha_pt;
128 }
129
130
131 /// Number of scalars/fields output by this element. Reimplements
132 /// broken virtual function in base class.
133 unsigned nscalar_paraview() const
134 {
135 return 2;
136 }
137
138 /// Write values of the i-th scalar field at the plot points. Needs
139 /// to be implemented for each new specific element type.
140 void scalar_value_paraview(std::ofstream& file_out,
141 const unsigned& i,
142 const unsigned& nplot) const
143 {
144 // Vector of local coordinates
145 Vector<double> s(DIM);
146
147 // Loop over plot points
148 unsigned num_plot_points = nplot_points_paraview(nplot);
149 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
150 {
151 // Get local coordinates of plot point
152 get_s_plot(iplot, nplot, s);
153 std::complex<double> u(interpolated_u_pml_helmholtz(s));
154
155 // Paraview need to ouput the fields separately so it loops through all
156 // the elements twice
157 switch (i)
158 {
159 // Real part first
160 case 0:
161 file_out << u.real() << std::endl;
162 break;
163
164 // Imaginary part second
165 case 1:
166 file_out << u.imag() << std::endl;
167 break;
168
169 // Never get here
170 default:
171#ifdef PARANOID
172 std::stringstream error_stream;
173 error_stream << "PML Helmholtz elements only store 2 fields (real "
174 "and imaginary) "
175 << "so i must be 0 or 1 rather "
176 << "than " << i << std::endl;
177 throw OomphLibError(error_stream.str(),
178 OOMPH_CURRENT_FUNCTION,
179 OOMPH_EXCEPTION_LOCATION);
180#endif
181 break;
182 }
183 }
184 }
185
186 /// Write values of the i-th scalar field at the plot points. Needs
187 /// to be implemented for each new specific element type.
189 std::ofstream& file_out,
190 const unsigned& i,
191 const unsigned& nplot,
192 FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
193 {
194 // Vector of local coordinates
195 Vector<double> s(DIM);
196
197 // Vector for coordinates
198 Vector<double> x(DIM);
199
200 // Exact solution Vector
201 Vector<double> exact_soln(2);
202
203 // Loop over plot points
204 unsigned num_plot_points = nplot_points_paraview(nplot);
205 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
206 {
207 // Get local coordinates of plot point
208 get_s_plot(iplot, nplot, s);
209
210 // Get x position as Vector
211 interpolated_x(s, x);
212
213 // Get exact solution at this point
214 (*exact_soln_pt)(x, exact_soln);
215
216 // Paraview need to output the fields separately so it loops through all
217 // the elements twice
218 switch (i)
219 {
220 // Real part first
221 case 0:
222 file_out << exact_soln[0] << std::endl;
223 break;
224
225 // Imaginary part second
226 case 1:
227 file_out << exact_soln[1] << std::endl;
228 break;
229
230 // Never get here
231 default:
232#ifdef PARANOID
233 std::stringstream error_stream;
234 error_stream << "PML Helmholtz elements only store 2 fields (real "
235 "and imaginary) "
236 << "so i must be 0 or 1 rather "
237 << "than " << i << std::endl;
238 throw OomphLibError(error_stream.str(),
239 OOMPH_CURRENT_FUNCTION,
240 OOMPH_EXCEPTION_LOCATION);
241#endif
242 break;
243 }
244 }
245 }
246
247 /// Name of the i-th scalar field. Default implementation
248 /// returns V1 for the first one, V2 for the second etc. Can (should!)
249 /// be overloaded with more meaningful names in specific elements.
250 std::string scalar_name_paraview(const unsigned& i) const
251 {
252 switch (i)
253 {
254 case 0:
255 return "Real part";
256 break;
257
258 case 1:
259 return "Imaginary part";
260 break;
261
262 // Never get here
263 default:
264#ifdef PARANOID
265 std::stringstream error_stream;
266 error_stream << "PML Helmholtz elements only store 2 fields (real "
267 "and imaginary) "
268 << "so i must be 0 or 1 rather "
269 << "than " << i << std::endl;
270 throw OomphLibError(error_stream.str(),
271 OOMPH_CURRENT_FUNCTION,
272 OOMPH_EXCEPTION_LOCATION);
273#endif
274
275 // Dummy return for the default statement
276 return " ";
277 break;
278 }
279 }
280
281 /// Output with default number of plot points
282 void output(std::ostream& outfile)
283 {
284 const unsigned n_plot = 5;
285 output(outfile, n_plot);
286 }
287
288 /// Output FE representation of soln: x,y,u_re,u_im or
289 /// x,y,z,u_re,u_im at n_plot^DIM plot points
290 void output(std::ostream& outfile, const unsigned& n_plot);
291
292
293 /// Output function for real part of full time-dependent solution
294 /// constructed by adding the scattered field
295 /// u = Re( (u_r +i u_i) exp(-i omega t)
296 /// at phase angle omega t = phi computed here, to the corresponding
297 /// incoming wave specified via the function pointer.
299 std::ostream& outfile,
300 FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt,
301 const double& phi,
302 const unsigned& nplot);
303
304
305 /// Output function for real part of full time-dependent solution
306 /// u = Re( (u_r +i u_i) exp(-i omega t))
307 /// at phase angle omega t = phi.
308 /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
309 /// direction
310 void output_real(std::ostream& outfile,
311 const double& phi,
312 const unsigned& n_plot);
313
314 /// Output function for imaginary part of full time-dependent
315 /// solution u = Im( (u_r +i u_i) exp(-i omega t) ) at phase angle omega t =
316 /// phi. x,y,u or x,y,z,u at n_plot plot points in each coordinate
317 /// direction
318 void output_imag(std::ostream& outfile,
319 const double& phi,
320 const unsigned& n_plot);
321
322 /// C_style output with default number of plot points
323 void output(FILE* file_pt)
324 {
325 const unsigned n_plot = 5;
326 output(file_pt, n_plot);
327 }
328
329 /// C-style output FE representation of soln: x,y,u_re,u_im or
330 /// x,y,z,u_re,u_im at n_plot^DIM plot points
331 void output(FILE* file_pt, const unsigned& n_plot);
332
333 /// Output exact soln: x,y,u_re_exact,u_im_exact
334 /// or x,y,z,u_re_exact,u_im_exact at n_plot^DIM plot points
335 void output_fct(std::ostream& outfile,
336 const unsigned& n_plot,
338
339 /// Output exact soln: (dummy time-dependent version to
340 /// keep intel compiler happy)
341 virtual void output_fct(
342 std::ostream& outfile,
343 const unsigned& n_plot,
344 const double& time,
346 {
347 throw OomphLibError(
348 "There is no time-dependent output_fct() for PMLHelmholtz elements.",
349 OOMPH_CURRENT_FUNCTION,
350 OOMPH_EXCEPTION_LOCATION);
351 }
352
353
354 /// Output function for real part of full time-dependent fct
355 /// u = Re( (u_r +i u_i) exp(-i omega t)
356 /// at phase angle omega t = phi.
357 /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
358 /// direction
359 void output_real_fct(std::ostream& outfile,
360 const double& phi,
361 const unsigned& n_plot,
363
364 /// Output function for imaginary part of full time-dependent fct
365 /// u = Im( (u_r +i u_i) exp(-i omega t))
366 /// at phase angle omega t = phi.
367 /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
368 /// direction
369 void output_imag_fct(std::ostream& outfile,
370 const double& phi,
371 const unsigned& n_plot,
373
374
375 /// Get error against and norm of exact solution
376 void compute_error(std::ostream& outfile,
378 double& error,
379 double& norm);
380
381
382 /// Dummy, time dependent error checker
383 void compute_error(std::ostream& outfile,
385 const double& time,
386 double& error,
387 double& norm)
388 {
389 throw OomphLibError(
390 "There is no time-dependent compute_error() for PMLHelmholtz elements.",
391 OOMPH_CURRENT_FUNCTION,
392 OOMPH_EXCEPTION_LOCATION);
393 }
394
395
396 /// Compute norm of fe solution
397 void compute_norm(double& norm);
398
399 /// Access function: Pointer to source function
401 {
402 return Source_fct_pt;
403 }
404
405 /// Access function: Pointer to source function. Const version
407 {
408 return Source_fct_pt;
409 }
410
411
412 /// Get source term at (Eulerian) position x. This function is
413 /// virtual to allow overloading in multi-physics problems where
414 /// the strength of the source function might be determined by
415 /// another system of equations.
416 inline virtual void get_source_helmholtz(const unsigned& ipt,
417 const Vector<double>& x,
418 std::complex<double>& source) const
419 {
420 // If no source function has been set, return zero
421 if (Source_fct_pt == 0)
422 {
423 source = std::complex<double>(0.0, 0.0);
424 }
425 else
426 {
427 // Get source strength
428 (*Source_fct_pt)(x, source);
429 }
430 }
431
432 /// Pure virtual function in which we specify the
433 /// values to be pinned (and set to zero) on the outer edge of
434 /// the pml layer. All of them! Vector is resized internally.
436 Vector<unsigned>& values_to_pin)
437 {
438 values_to_pin.resize(2);
439
440 values_to_pin[0] = 0;
441 values_to_pin[1] = 1;
442 }
443
444
445 /// Get flux: flux[i] = du/dx_i for real and imag part
447 Vector<std::complex<double>>& flux) const
448 {
449 // Find out how many nodes there are in the element
450 const unsigned n_node = nnode();
451
452
453 // Set up memory for the shape and test functions
454 Shape psi(n_node);
455 DShape dpsidx(n_node, DIM);
456
457 // Call the derivatives of the shape and test functions
458 dshape_eulerian(s, psi, dpsidx);
459
460 // Initialise to zero
461 const std::complex<double> zero(0.0, 0.0);
462 for (unsigned j = 0; j < DIM; j++)
463 {
464 flux[j] = zero;
465 }
466
467 // Loop over nodes
468 for (unsigned l = 0; l < n_node; l++)
469 {
470 // Cache the complex value of the unknown
471 const std::complex<double> u_value(
472 this->nodal_value(l, u_index_helmholtz().real()),
473 this->nodal_value(l, u_index_helmholtz().imag()));
474 // Loop over derivative directions
475 for (unsigned j = 0; j < DIM; j++)
476 {
477 flux[j] += u_value * dpsidx(l, j);
478 }
479 }
480 }
481
482
483 /// Add the element's contribution to its residual vector (wrapper)
485 {
486 // Call the generic residuals function with flag set to 0
487 // using a dummy matrix argument
490 }
491
492
493 /// Add the element's contribution to its residual vector and
494 /// element Jacobian matrix (wrapper)
496 DenseMatrix<double>& jacobian)
497 {
498 // Call the generic routine with the flag set to 1
500 }
501
502
503 /// Return FE representation of function value u_helmholtz(s)
504 /// at local coordinate s
505 inline std::complex<double> interpolated_u_pml_helmholtz(
506 const Vector<double>& s) const
507 {
508 // Find number of nodes
509 const unsigned n_node = nnode();
510
511 // Local shape function
512 Shape psi(n_node);
513
514 // Find values of shape function
515 shape(s, psi);
516
517 // Initialise value of u
518 std::complex<double> interpolated_u(0.0, 0.0);
519
520 // Get the index at which the helmholtz unknown is stored
521 const unsigned u_nodal_index_real = u_index_helmholtz().real();
522 const unsigned u_nodal_index_imag = u_index_helmholtz().imag();
523
524 // Loop over the local nodes and sum
525 for (unsigned l = 0; l < n_node; l++)
526 {
527 // Make a temporary complex number from the stored data
528 const std::complex<double> u_value(
529 this->nodal_value(l, u_nodal_index_real),
530 this->nodal_value(l, u_nodal_index_imag));
531 // Add to the interpolated value
532 interpolated_u += u_value * psi[l];
533 }
534 return interpolated_u;
535 }
536
537
538 /// Self-test: Return 0 for OK
539 unsigned self_test();
540
541
542 /// Compute pml coefficients at position x and integration point ipt.
543 /// pml_laplace_factor is used in the residual contribution from the laplace
544 /// operator, similarly pml_k_squared_factor is used in the contribution
545 /// from the k^2 of the Helmholtz operator.
547 const unsigned& ipt,
548 const Vector<double>& x,
549 Vector<std::complex<double>>& pml_laplace_factor,
550 std::complex<double>& pml_k_squared_factor)
551 {
552 /// Vector which points from the inner boundary to x
553 Vector<double> nu(DIM);
554 for (unsigned k = 0; k < DIM; k++)
555 {
556 nu[k] = x[k] - this->Pml_inner_boundary[k];
557 }
558
559 /// Vector which points from the inner boundary to the edge of the
560 /// boundary
561 Vector<double> pml_width(DIM);
562 for (unsigned k = 0; k < DIM; k++)
563 {
564 pml_width[k] =
565 this->Pml_outer_boundary[k] - this->Pml_inner_boundary[k];
566 }
567
568 // Declare gamma_i vectors of complex numbers for PML weights
569 Vector<std::complex<double>> pml_gamma(DIM);
570
571 if (this->Pml_is_enabled)
572 {
573 // Cache k_squared to pass into mapping function
574 double k_squared_local = k_squared();
575
576 for (unsigned k = 0; k < DIM; k++)
577 {
578 // If PML is enabled in the respective direction
579 if (this->Pml_direction_active[k])
580 {
581 pml_gamma[k] = Pml_mapping_pt->gamma(
582 nu[k], pml_width[k], k_squared_local, alpha());
583 }
584 else
585 {
586 pml_gamma[k] = 1.0;
587 }
588 }
589
590 /// for 2D, in order:
591 /// g_y/g_x, g_x/g_y for Laplace bit and g_x*g_y for Helmholtz bit
592 /// for 3D, in order: g_y*g_x/g_x, g*x*g_z/g_y, g_x*g_y/g_z for Laplace
593 /// bit and g_x*g_y*g_z for Helmholtz bit
594 if (DIM == 2)
595 {
596 pml_laplace_factor[0] = pml_gamma[1] / pml_gamma[0];
597 pml_laplace_factor[1] = pml_gamma[0] / pml_gamma[1];
598 pml_k_squared_factor = pml_gamma[0] * pml_gamma[1];
599 }
600 else // if (DIM == 3)
601 {
602 pml_laplace_factor[0] = pml_gamma[1] * pml_gamma[2] / pml_gamma[0];
603 pml_laplace_factor[1] = pml_gamma[0] * pml_gamma[2] / pml_gamma[1];
604 pml_laplace_factor[2] = pml_gamma[0] * pml_gamma[1] / pml_gamma[2];
605 pml_k_squared_factor = pml_gamma[0] * pml_gamma[1] * pml_gamma[2];
606 }
607 }
608 else
609 {
610 /// The weights all default to 1.0 as if the propagation
611 /// medium is the physical domain
612 for (unsigned k = 0; k < DIM; k++)
613 {
614 pml_laplace_factor[k] = std::complex<double>(1.0, 0.0);
615 }
616
617 pml_k_squared_factor = std::complex<double>(1.0, 0.0);
618 }
619 }
620
621 /// Return a pointer to the PML Mapping object
623 {
624 return Pml_mapping_pt;
625 }
626
627 /// Return a pointer to the PML Mapping object (const version)
629 {
630 return Pml_mapping_pt;
631 }
632
633 /// Static so that the class doesn't need to instantiate a new default
634 /// everytime it uses it
636
637 /// The number of "DOF types" that degrees of freedom in this element
638 /// are sub-divided into: real and imaginary part
639 unsigned ndof_types() const
640 {
641 return 2;
642 }
643
644 /// Create a list of pairs for all unknowns in this element,
645 /// so that the first entry in each pair contains the global equation
646 /// number of the unknown, while the second one contains the number
647 /// of the "DOF type" that this unknown is associated with.
648 /// (Function can obviously only be called if the equation numbering
649 /// scheme has been set up.) Real=0; Imag=1
651 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
652 {
653 // temporary pair (used to store dof lookup prior to being added to list)
654 std::pair<unsigned, unsigned> dof_lookup;
655
656 // number of nodes
657 unsigned n_node = this->nnode();
658
659 // loop over the nodes
660 for (unsigned n = 0; n < n_node; n++)
661 {
662 // determine local eqn number for real part
663 int local_eqn_number =
664 this->nodal_local_eqn(n, u_index_helmholtz().real());
665
666 // ignore pinned values
667 if (local_eqn_number >= 0)
668 {
669 // store dof lookup in temporary pair: First entry in pair
670 // is global equation number; second entry is dof type
671 dof_lookup.first = this->eqn_number(local_eqn_number);
672 dof_lookup.second = 0;
673
674 // add to list
675 dof_lookup_list.push_front(dof_lookup);
676 }
677
678 // determine local eqn number for imag part
680
681 // ignore pinned values
682 if (local_eqn_number >= 0)
683 {
684 // store dof lookup in temporary pair: First entry in pair
685 // is global equation number; second entry is dof type
686 dof_lookup.first = this->eqn_number(local_eqn_number);
687 dof_lookup.second = 1;
688
689 // add to list
690 dof_lookup_list.push_front(dof_lookup);
691 }
692 }
693 }
694
695
696 protected:
697 /// Shape/test functions and derivs w.r.t. to global coords at
698 /// local coord. s; return Jacobian of mapping
700 const Vector<double>& s,
701 Shape& psi,
702 DShape& dpsidx,
703 Shape& test,
704 DShape& dtestdx) const = 0;
705
706
707 /// Shape/test functions and derivs w.r.t. to global coords at
708 /// integration point ipt; return Jacobian of mapping
710 const unsigned& ipt,
711 Shape& psi,
712 DShape& dpsidx,
713 Shape& test,
714 DShape& dtestdx) const = 0;
715
716 /// Compute element residual Vector only (if flag=and/or element
717 /// Jacobian matrix
719 Vector<double>& residuals,
720 DenseMatrix<double>& jacobian,
721 const unsigned& flag);
722
723 /// Pointer to wavenumber complex shift
724 double* Alpha_pt;
725
726 /// Pointer to source function:
728
729 /// Pointer to wave number (must be set!)
731
732 /// Pointer to class which holds the pml mapping function, also known as
733 /// gamma
735
736 /// Static default value for the physical constants (initialised to zero)
738 };
739
740
741 /// ////////////////////////////////////////////////////////////////////////
742 /// ////////////////////////////////////////////////////////////////////////
743 /// ////////////////////////////////////////////////////////////////////////
744
745
746 //======================================================================
747 /// QPMLHelmholtzElement elements are linear/quadrilateral/
748 /// brick-shaped PMLHelmholtz elements with isoparametric
749 /// interpolation for the function.
750 //======================================================================
751 template<unsigned DIM, unsigned NNODE_1D>
752 class QPMLHelmholtzElement : public virtual QElement<DIM, NNODE_1D>,
753 public virtual PMLHelmholtzEquations<DIM>
754 {
755 private:
756 /// Static int that holds the number of variables at
757 /// nodes: always the same
758 static const unsigned Initial_Nvalue;
759
760 public:
761 /// Constructor: Call constructors for QElement and
762 /// PMLHelmholtz equations
764 : QElement<DIM, NNODE_1D>(), PMLHelmholtzEquations<DIM>()
765 {
766 }
767
768 /// Broken copy constructor
770 delete;
771
772 /// Broken assignment operator
773 /*void operator=(const QPMLHelmholtzElement<DIM,NNODE_1D>&) = delete;*/
774
775 /// Required # of `values' (pinned or dofs)
776 /// at node n
777 inline unsigned required_nvalue(const unsigned& n) const
778 {
779 return Initial_Nvalue;
780 }
781
782 /// Output function:
783 /// x,y,u or x,y,z,u
784 void output(std::ostream& outfile)
785 {
787 }
788
789
790 /// Output function:
791 /// x,y,u or x,y,z,u at n_plot^DIM plot points
792 void output(std::ostream& outfile, const unsigned& n_plot)
793 {
794 PMLHelmholtzEquations<DIM>::output(outfile, n_plot);
795 }
796
797 /// Output function for real part of full time-dependent solution
798 /// u = Re( (u_r +i u_i) exp(-i omega t)
799 /// at phase angle omega t = phi.
800 /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
801 /// direction
802 void output_real(std::ostream& outfile,
803 const double& phi,
804 const unsigned& n_plot)
805 {
806 PMLHelmholtzEquations<DIM>::output_real(outfile, phi, n_plot);
807 }
808
809 /// Output function for imaginary part of full time-dependent
810 /// solution u = Im( (u_r +i u_i) exp(-i omega t)) at phase angle omega t =
811 /// phi. x,y,u or x,y,z,u at n_plot plot points in each coordinate
812 /// direction
813 void output_imag(std::ostream& outfile,
814 const double& phi,
815 const unsigned& n_plot)
816 {
817 PMLHelmholtzEquations<DIM>::output_imag(outfile, phi, n_plot);
818 }
819
820
821 /// C-style output function:
822 /// x,y,u or x,y,z,u
823 void output(FILE* file_pt)
824 {
826 }
827
828
829 /// C-style output function:
830 /// x,y,u or x,y,z,u at n_plot^DIM plot points
831 void output(FILE* file_pt, const unsigned& n_plot)
832 {
833 PMLHelmholtzEquations<DIM>::output(file_pt, n_plot);
834 }
835
836
837 /// Output function for an exact solution:
838 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
839 void output_fct(std::ostream& outfile,
840 const unsigned& n_plot,
842 {
843 PMLHelmholtzEquations<DIM>::output_fct(outfile, n_plot, exact_soln_pt);
844 }
845
846
847 /// Output function for real part of full time-dependent fct
848 /// u = Re( (u_r +i u_i) exp(-i omega t)
849 /// at phase angle omega t = phi.
850 /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
851 /// direction
852 void output_real_fct(std::ostream& outfile,
853 const double& phi,
854 const unsigned& n_plot,
856 {
858 outfile, phi, n_plot, exact_soln_pt);
859 }
860
861 /// Output function for imaginary part of full time-dependent fct
862 /// u = Im( (u_r +i u_i) exp(-i omega t))
863 /// at phase angle omega t = phi.
864 /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
865 /// direction
866 void output_imag_fct(std::ostream& outfile,
867 const double& phi,
868 const unsigned& n_plot,
870 {
872 outfile, phi, n_plot, exact_soln_pt);
873 }
874
875
876 /// Output function for a time-dependent exact solution.
877 /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
878 /// (Calls the steady version)
879 void output_fct(std::ostream& outfile,
880 const unsigned& n_plot,
881 const double& time,
883 {
885 outfile, n_plot, time, exact_soln_pt);
886 }
887
888 protected:
889 /// Shape, test functions & derivs. w.r.t. to global coords. Return
890 /// Jacobian.
892 Shape& psi,
893 DShape& dpsidx,
894 Shape& test,
895 DShape& dtestdx) const;
896
897
898 /// Shape, test functions & derivs. w.r.t. to global coords. at
899 /// integration point ipt. Return Jacobian.
901 const unsigned& ipt,
902 Shape& psi,
903 DShape& dpsidx,
904 Shape& test,
905 DShape& dtestdx) const;
906 };
907
908
909 // Inline functions:
910
911
912 //======================================================================
913 /// Define the shape functions and test functions and derivatives
914 /// w.r.t. global coordinates and return Jacobian of mapping.
915 ///
916 /// Galerkin: Test functions = shape functions
917 //======================================================================
918 template<unsigned DIM, unsigned NNODE_1D>
921 Shape& psi,
922 DShape& dpsidx,
923 Shape& test,
924 DShape& dtestdx) const
925 {
926 // Call the geometrical shape functions and derivatives
927 const double J = this->dshape_eulerian(s, psi, dpsidx);
928
929 // Set the test functions equal to the shape functions
930 test = psi;
931 dtestdx = dpsidx;
932
933 // Return the jacobian
934 return J;
935 }
936
937
938 //======================================================================
939 /// Define the shape functions and test functions and derivatives
940 /// w.r.t. global coordinates and return Jacobian of mapping.
941 ///
942 /// Galerkin: Test functions = shape functions
943 //======================================================================
944 template<unsigned DIM, unsigned NNODE_1D>
947 Shape& psi,
948 DShape& dpsidx,
949 Shape& test,
950 DShape& dtestdx) const
951 {
952 // Call the geometrical shape functions and derivatives
953 const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
954
955 // Set the pointers of the test functions
956 test = psi;
957 dtestdx = dpsidx;
958
959 // Return the jacobian
960 return J;
961 }
962
963 /// /////////////////////////////////////////////////////////////////////
964 /// /////////////////////////////////////////////////////////////////////
965 /// /////////////////////////////////////////////////////////////////////
966
967
968 //=======================================================================
969 /// Face geometry for the QPMLHelmholtzElement elements: The spatial
970 /// dimension of the face elements is one lower than that of the
971 /// bulk element but they have the same number of points
972 /// along their 1D edges.
973 //=======================================================================
974 template<unsigned DIM, unsigned NNODE_1D>
975 class FaceGeometry<QPMLHelmholtzElement<DIM, NNODE_1D>>
976 : public virtual QElement<DIM - 1, NNODE_1D>
977 {
978 public:
979 /// Constructor: Call the constructor for the
980 /// appropriate lower-dimensional QElement
981 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
982 };
983
984 /// /////////////////////////////////////////////////////////////////////
985 /// /////////////////////////////////////////////////////////////////////
986 /// /////////////////////////////////////////////////////////////////////
987
988
989 //=======================================================================
990 /// Face geometry for the 1D QPMLHelmholtzElement elements:
991 /// Point elements
992 //=======================================================================
993 template<unsigned NNODE_1D>
995 : public virtual PointElement
996 {
997 public:
998 /// Constructor: Call the constructor for the
999 /// appropriate lower-dimensional QElement
1001 };
1002
1003
1004 /// /////////////////////////////////////////////////////////////////////
1005 /// /////////////////////////////////////////////////////////////////////
1006 /// /////////////////////////////////////////////////////////////////////
1007
1008
1009 //==========================================================
1010 /// PMLHelmholtz upgraded to become projectable
1011 //==========================================================
1012 template<class HELMHOLTZ_ELEMENT>
1014 : public virtual ProjectableElement<HELMHOLTZ_ELEMENT>
1015 {
1016 public:
1017 /// Constructor [this was only required explicitly
1018 /// from gcc 4.5.2 onwards...]
1020
1021 /// Specify the values associated with field fld.
1022 /// The information is returned in a vector of pairs which comprise
1023 /// the Data object and the value within it, that correspond to field fld.
1025 {
1026#ifdef PARANOID
1027 if (fld > 1)
1028 {
1029 std::stringstream error_stream;
1030 error_stream << "PMLHelmholtz elements only store 2 fields so fld = "
1031 << fld << " is illegal \n";
1032 throw OomphLibError(
1033 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1034 }
1035#endif
1036
1037 // Create the vector
1038 unsigned nnod = this->nnode();
1039 Vector<std::pair<Data*, unsigned>> data_values(nnod);
1040
1041 // Loop over all nodes
1042 for (unsigned j = 0; j < nnod; j++)
1043 {
1044 // Add the data value associated field: The node itself
1045 data_values[j] = std::make_pair(this->node_pt(j), fld);
1046 }
1047
1048 // Return the vector
1049 return data_values;
1050 }
1051
1052 /// Number of fields to be projected: 2 (real and imag part)
1054 {
1055 return 2;
1056 }
1057
1058 /// Number of history values to be stored for fld-th field.
1059 /// (Note: count includes current value!)
1060 unsigned nhistory_values_for_projection(const unsigned& fld)
1061 {
1062#ifdef PARANOID
1063 if (fld > 1)
1064 {
1065 std::stringstream error_stream;
1066 error_stream << "PMLHelmholtz elements only store two fields so fld = "
1067 << fld << " is illegal\n";
1068 throw OomphLibError(
1069 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1070 }
1071#endif
1072 return this->node_pt(0)->ntstorage();
1073 }
1074
1075 /// Number of positional history values
1076 /// (Note: count includes current value!)
1078 {
1079 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
1080 }
1081
1082 /// Return Jacobian of mapping and shape functions of field fld
1083 /// at local coordinate s
1084 double jacobian_and_shape_of_field(const unsigned& fld,
1085 const Vector<double>& s,
1086 Shape& psi)
1087 {
1088#ifdef PARANOID
1089 if (fld > 1)
1090 {
1091 std::stringstream error_stream;
1092 error_stream << "PMLHelmholtz elements only store two fields so fld = "
1093 << fld << " is illegal.\n";
1094 throw OomphLibError(
1095 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1096 }
1097#endif
1098 unsigned n_dim = this->dim();
1099 unsigned n_node = this->nnode();
1100 Shape test(n_node);
1101 DShape dpsidx(n_node, n_dim), dtestdx(n_node, n_dim);
1102 double J = this->dshape_and_dtest_eulerian_helmholtz(
1103 s, psi, dpsidx, test, dtestdx);
1104 return J;
1105 }
1106
1107
1108 /// Return interpolated field fld at local coordinate s, at time
1109 /// level t (t=0: present; t>0: history values)
1110 double get_field(const unsigned& t,
1111 const unsigned& fld,
1112 const Vector<double>& s)
1113 {
1114#ifdef PARANOID
1115 if (fld > 1)
1116 {
1117 std::stringstream error_stream;
1118 error_stream << "PMLHelmholtz elements only store two fields so fld = "
1119 << fld << " is illegal\n";
1120 throw OomphLibError(
1121 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1122 }
1123#endif
1124 // Find the index at which the variable is stored
1125 std::complex<unsigned> complex_u_nodal_index = this->u_index_helmholtz();
1126 unsigned u_nodal_index = 0;
1127 if (fld == 0)
1128 {
1129 u_nodal_index = complex_u_nodal_index.real();
1130 }
1131 else
1132 {
1133 u_nodal_index = complex_u_nodal_index.imag();
1134 }
1135
1136
1137 // Local shape function
1138 unsigned n_node = this->nnode();
1139 Shape psi(n_node);
1140
1141 // Find values of shape function
1142 this->shape(s, psi);
1143
1144 // Initialise value of u
1145 double interpolated_u = 0.0;
1146
1147 // Sum over the local nodes
1148 for (unsigned l = 0; l < n_node; l++)
1149 {
1150 interpolated_u += this->nodal_value(t, l, u_nodal_index) * psi[l];
1151 }
1152 return interpolated_u;
1153 }
1154
1155
1156 /// Return number of values in field fld: One per node
1157 unsigned nvalue_of_field(const unsigned& fld)
1158 {
1159#ifdef PARANOID
1160 if (fld > 1)
1161 {
1162 std::stringstream error_stream;
1163 error_stream << "PMLHelmholtz elements only store two fields so fld = "
1164 << fld << " is illegal\n";
1165 throw OomphLibError(
1166 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1167 }
1168#endif
1169 return this->nnode();
1170 }
1171
1172
1173 /// Return local equation number of value j in field fld.
1174 int local_equation(const unsigned& fld, const unsigned& j)
1175 {
1176#ifdef PARANOID
1177 if (fld > 1)
1178 {
1179 std::stringstream error_stream;
1180 error_stream << "PMLHelmholtz elements only store two fields so fld = "
1181 << fld << " is illegal\n";
1182 throw OomphLibError(
1183 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1184 }
1185#endif
1186 std::complex<unsigned> complex_u_nodal_index = this->u_index_helmholtz();
1187 unsigned u_nodal_index = 0;
1188 if (fld == 0)
1189 {
1190 u_nodal_index = complex_u_nodal_index.real();
1191 }
1192 else
1193 {
1194 u_nodal_index = complex_u_nodal_index.imag();
1195 }
1196 return this->nodal_local_eqn(j, u_nodal_index);
1197 }
1198
1199 /// Output FE representation of soln: x,y,u or x,y,z,u at
1200 /// n_plot^DIM plot points
1201 void output(std::ostream& outfile, const unsigned& nplot)
1202 {
1203 HELMHOLTZ_ELEMENT::output(outfile, nplot);
1204 }
1205 };
1206
1207
1208 //=======================================================================
1209 /// Face geometry for element is the same as that for the underlying
1210 /// wrapped element
1211 //=======================================================================
1212 template<class ELEMENT>
1214 : public virtual FaceGeometry<ELEMENT>
1215 {
1216 public:
1217 FaceGeometry() : FaceGeometry<ELEMENT>() {}
1218 };
1219
1220 /// /////////////////////////////////////////////////////////////////////
1221 /// /////////////////////////////////////////////////////////////////////
1222 /// /////////////////////////////////////////////////////////////////////
1223
1224 //=======================================================================
1225 /// Policy class defining the elements to be used in the actual
1226 /// PML layers. Same!
1227 //=======================================================================
1228 template<unsigned NNODE_1D>
1230 : public virtual QPMLHelmholtzElement<2, NNODE_1D>
1231 {
1232 public:
1233 /// Constructor: Call the constructor for the
1234 /// appropriate QElement
1236 };
1237
1238} // End of Namespace oomph
1239
1240#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A mapping function propsed by Bermudez et al, appears to be the best for the Helmholtz equations and ...
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:879
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2862
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number....
Definition: elements.h:726
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
An OomphLibError object which should be thrown when an run-time error is encountered....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: pml_meshes.h:60
bool Pml_is_enabled
Boolean indicating if element is used in pml mode.
Definition: pml_meshes.h:119
std::vector< bool > Pml_direction_active
Coordinate direction along which pml boundary is constant; alternatively: coordinate direction in whi...
Definition: pml_meshes.h:124
Vector< double > Pml_outer_boundary
Coordinate of outer pml boundary (Storage is provided for any coordinate direction; only the entries ...
Definition: pml_meshes.h:134
Vector< double > Pml_inner_boundary
Coordinate of inner pml boundary (Storage is provided for any coordinate direction; only the entries ...
Definition: pml_meshes.h:129
A class for all isoparametric elements that solve the Helmholtz equations with pml capabilities....
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and element Jacobian matrix (wrapper)
PMLHelmholtzSourceFctPt Source_fct_pt
Pointer to source function:
void values_to_be_pinned_on_outer_pml_boundary(Vector< unsigned > &values_to_pin)
Pure virtual function in which we specify the values to be pinned (and set to zero) on the outer edge...
virtual void fill_in_generic_residual_contribution_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
virtual double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
void output(std::ostream &outfile)
Output with default number of plot points.
PMLMapping * Pml_mapping_pt
Pointer to class which holds the pml mapping function, also known as gamma.
virtual void get_source_helmholtz(const unsigned &ipt, const Vector< double > &x, std::complex< double > &source) const
Get source term at (Eulerian) position x. This function is virtual to allow overloading in multi-phys...
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: (dummy time-dependent version to keep intel compiler happy)
virtual double dshape_and_dtest_eulerian_helmholtz(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
const double & alpha() const
Alpha, wavenumber complex shift.
double k_squared()
Get the square of wavenumber.
void compute_pml_coefficients(const unsigned &ipt, const Vector< double > &x, Vector< std::complex< double > > &pml_laplace_factor, std::complex< double > &pml_k_squared_factor)
Compute pml coefficients at position x and integration point ipt. pml_laplace_factor is used in the r...
PMLHelmholtzSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
PMLHelmholtzEquations(const PMLHelmholtzEquations &dummy)=delete
Broken copy constructor.
void scalar_value_fct_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
void compute_norm(double &norm)
Compute norm of fe solution.
void output_imag_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for imaginary part of full time-dependent fct u = Im( (u_r +i u_i) exp(-i omega t)) a...
void output_real(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for real part of full time-dependent solution u = Re( (u_r +i u_i) exp(-i omega t)) a...
void output(FILE *file_pt)
C_style output with default number of plot points.
double *& alpha_pt()
Pointer to Alpha, wavenumber complex shift.
static BermudezPMLMapping Default_pml_mapping
Static so that the class doesn't need to instantiate a new default everytime it uses it.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void(* PMLHelmholtzSourceFctPt)(const Vector< double > &x, std::complex< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: real and imag...
double * Alpha_pt
Pointer to wavenumber complex shift.
double * K_squared_pt
Pointer to wave number (must be set!)
PMLMapping *& pml_mapping_pt()
Return a pointer to the PML Mapping object.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
PMLMapping *const & pml_mapping_pt() const
Return a pointer to the PML Mapping object (const version)
double *& k_squared_pt()
Get pointer to k_squared.
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class.
void output_imag(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for imaginary part of full time-dependent solution u = Im( (u_r +i u_i) exp(-i omega ...
static double Default_Physical_Constant_Value
Static default value for the physical constants (initialised to zero)
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_re_exact,u_im_exact or x,y,z,u_re_exact,u_im_exact at n_plot^DIM plot points...
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
void get_flux(const Vector< double > &s, Vector< std::complex< double > > &flux) const
Get flux: flux[i] = du/dx_i for real and imag part.
unsigned self_test()
Self-test: Return 0 for OK.
PMLHelmholtzSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
void output_real_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for real part of full time-dependent fct u = Re( (u_r +i u_i) exp(-i omega t) at phas...
void output_total_real(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt incoming_wave_fct_pt, const double &phi, const unsigned &nplot)
Output function for real part of full time-dependent solution constructed by adding the scattered fie...
std::complex< double > interpolated_u_pml_helmholtz(const Vector< double > &s) const
Return FE representation of function value u_helmholtz(s) at local coordinate s.
PMLLayerElement()
Constructor: Call the constructor for the appropriate QElement.
General definition of policy class defining the elements to be used in the actual PML layers....
Definition: pml_meshes.h:48
Class to hold the mapping function (gamma) for the Pml which defines how the coordinates are transfor...
virtual std::complex< double > gamma(const double &nu_i, const double &pml_width_i, const double &wavenumber_squared, const double &alpha_shift=0.0)=0
Pure virtual to return Pml mapping gamma, where gamma is the as function of where where h is the v...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
ProjectablePMLHelmholtzElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!...
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
unsigned nfields_for_projection()
Number of fields to be projected: 2 (real and imag part)
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln: x,y,u or x,y,z,u at n_plot^DIM plot points.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
QPMLHelmholtzElement(const QPMLHelmholtzElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot ...
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
void output_real(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for real part of full time-dependent solution u = Re( (u_r +i u_i) exp(-i omega t) at...
void output_imag(std::ostream &outfile, const double &phi, const unsigned &n_plot)
Output function for imaginary part of full time-dependent solution u = Im( (u_r +i u_i) exp(-i omega ...
double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt....
QPMLHelmholtzElement()
Constructor: Call constructors for QElement and PMLHelmholtz equations.
double dshape_and_dtest_eulerian_helmholtz(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
void output_real_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for real part of full time-dependent fct u = Re( (u_r +i u_i) exp(-i omega t) at phas...
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
void output_imag_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for imaginary part of full time-dependent fct u = Im( (u_r +i u_i) exp(-i omega t)) a...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
void output()
Doc the command line arguments.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...