pml_fourier_decomposed_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 Fourier-decomposed Helmholtz elements
27#ifndef OOMPH_PML_FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENTS_HEADER
28#define OOMPH_PML_FOURIER_DECOMPOSED_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
40// OOMPH-LIB headers
41#include "../generic/projection.h"
42#include "../generic/nodes.h"
43#include "../generic/Qelements.h"
44#include "../generic/oomph_utilities.h"
45#include "../generic/pml_meshes.h"
46#include "../generic/projection.h"
47#include "../generic/oomph_definitions.h"
48
49
50namespace oomph
51{
52 //========================================================================
53 /// Helper namespace for functions required for Helmholtz computations
54 //========================================================================
55 namespace Legendre_functions_helper
56 {
57 /// Factorial
58 extern double factorial(const unsigned& l);
59
60 /// Legendre polynomials depending on one parameter
61 extern double plgndr1(const unsigned& n, const double& x);
62
63 /// Legendre polynomials depending on two parameters
64 extern double plgndr2(const unsigned& l,
65 const unsigned& m,
66 const double& x);
67
68 } // namespace Legendre_functions_helper
69
70
71 /// ////////////////////////////////////////////////////////////////////
72 /// ////////////////////////////////////////////////////////////////////
73 /// ////////////////////////////////////////////////////////////////////
74
75 //=======================================================================
76 /// Class to hold the mapping function for the PML
77 ///
78 //=======================================================================
80 {
81 public:
82 /// Default constructor (empty)
84
85 /// Pure virtual to return PML mapping gamma, where gamma is the
86 /// \f$d\tilde x / d x\f$ as function of \f$\nu\f$ where \f$\nu = x - h\f$
87 /// where h is the vector from the origin to the start of the PML
88 virtual std::complex<double> gamma(const double& nu_i,
89 const double& pml_width_i,
90 const double& k_squared) = 0;
91
92 /// Pure virtual to return PML transformed coordinate, also known as
93 /// \f$d\tilde x \f$ as function of \f$\nu\f$ where \f$\nu = x - h\f$ where
94 /// h is the vector from the origin to the start of the PML
95 virtual std::complex<double> transformed_coordinate(
96 const double& nu_i,
97 const double& pml_width_i,
98 const double& pml_inner_boundary,
99 const double& k_squared) = 0;
100 };
101
102
103 //=======================================================================
104 /// The mapping function propsed by Bermudez et al, appears to be the best
105 /// and so this will be the default mapping (see definition of
106 /// PMLHelmholtzEquations)
107 //=======================================================================
110 {
111 public:
112 /// Default constructor (empty)
114
115 /// Overwrite the pure PML mapping coefficient function to return the
116 /// mapping function proposed by Bermudez et al
117 std::complex<double> gamma(const double& nu_i,
118 const double& pml_width_i,
119 const double& k_squared)
120 {
121 /// return \f$\gamma=1 + (1/k)(i/|outer_boundary - x|)\f$ or more
122 /// abstractly \f$\gamma = 1 + \frac i {k\delta_{pml}}(1/|1-\bar\nu|)\f$
123 return 1.0 +
124 std::complex<double>(1.0 / sqrt(k_squared), 0) *
125 std::complex<double>(0.0, 1.0 / (std::fabs(pml_width_i - nu_i)));
126 }
127
128 /// Overwrite the pure PML mapping coefficient function to return the
129 /// transformed coordinate proposed by Bermudez et al
130 std::complex<double> transformed_coordinate(
131 const double& nu_i,
132 const double& pml_width_i,
133 const double& pml_inner_boundary,
134 const double& k_squared)
135 {
136 /// return \f$\tilde x = h + \nu + \log(1-|\nu / \delta|)\f$
137 double log_arg = 1.0 - std::fabs(nu_i / pml_width_i);
138 return std::complex<double>(pml_inner_boundary + nu_i,
139 -log(log_arg) / sqrt(k_squared));
140 }
141 };
142
143
144 /// /////////////////////////////////////////////////////////////////////
145 /// /////////////////////////////////////////////////////////////////////
146 /// /////////////////////////////////////////////////////////////////////
147
148
149 //=============================================================
150 /// A class for all isoparametric elements that solve the
151 /// Helmholtz equations with pml capabilities.
152 /// in Fourier decomposed form (cylindrical polars):
153 /// \f[ U(r,\varphi,z) = \Re( u^{(n)}(r,z) \exp(-i n \varphi)) \f]
154 /// We are solving for \f$ u^{(n)}(r,z)\f$ for given parameters
155 /// \f$ k^2 \f$ and \f$ n \f$ .
156 /// This contains the generic maths. Shape functions, geometric
157 /// mapping etc. must get implemented in derived class.
158 //=============================================================
160 : public virtual PMLElementBase<2>,
161 public virtual FiniteElement
162 {
163 public:
164 /// Function pointer to source function fct(x,f(x)) --
165 /// x is a Vector!
167 const Vector<double>& x, std::complex<double>& f);
168
169 /// Constructor
172 {
173 // Provide pointer to static method (Save memory)
177
179 }
180
181
182 /// Broken copy constructor
184 const PMLFourierDecomposedHelmholtzEquations& dummy) = delete;
185
186 /// Broken assignment operator
187 // Commented out broken assignment operator because this can lead to a
188 // conflict warning when used in the virtual inheritence hierarchy.
189 // Essentially the compiler doesn't realise that two separate
190 // implementations of the broken function are the same and so, quite
191 // rightly, it shouts.
192 /*void operator=(const PMLFourierDecomposedHelmholtzEquations&) = delete;*/
193
194 /// Return the index at which the unknown value
195 /// is stored: Real/imag part of index contains (real) index of
196 /// real/imag part.
197 virtual inline std::complex<unsigned> u_index_pml_fourier_decomposed_helmholtz()
198 const
199 {
200 return std::complex<unsigned>(0, 1);
201 }
202
203
204 /// Get pointer to frequency
205 double*& k_squared_pt()
206 {
207 return K_squared_pt;
208 }
209
210
211 /// Get k squared
212 double k_squared()
213 {
214#ifdef PARANOID
215 if (K_squared_pt == 0)
216 {
217 throw OomphLibError(
218 "Please set pointer to k_squared using access fct to pointer!",
219 OOMPH_CURRENT_FUNCTION,
220 OOMPH_EXCEPTION_LOCATION);
221 }
222#endif
223 return *K_squared_pt;
224 }
225
226 /// Get pointer to complex shift
227 double*& alpha_pt()
228 {
229 return Alpha_pt;
230 }
231
232
233 /// Get complex shift
234 double alpha()
235 {
236 return *Alpha_pt;
237 }
238
239 /// Get pointer to Fourier wavenumber
241 {
242 return N_pml_fourier_pt;
243 }
244
245 /// Get the Fourier wavenumber
247 {
248 if (N_pml_fourier_pt == 0)
249 {
250 return 0;
251 }
252 else
253 {
254 return *N_pml_fourier_pt;
255 }
256 }
257
258
259 /// Output with default number of plot points
260 void output(std::ostream& outfile)
261 {
262 const unsigned n_plot = 5;
263 output(outfile, n_plot);
264 }
265
266 /// Output FE representation of soln: x,y,u_re,u_im or
267 /// x,y,z,u_re,u_im at n_plot^2 plot points
268 void output(std::ostream& outfile, const unsigned& n_plot);
269
270 /// Output function for real part of full time-dependent solution
271 /// u = Re( (u_r +i u_i) exp(-i omega t)
272 /// at phase angle omega t = phi.
273 /// r,z,u at n_plot plot points in each coordinate
274 /// direction
275 void output_real(std::ostream& outfile,
276 const double& phi,
277 const unsigned& n_plot);
278
279 /// C_style output with default number of plot points
280 void output(FILE* file_pt)
281 {
282 const unsigned n_plot = 5;
283 output(file_pt, n_plot);
284 }
285
286 /// C-style output FE representation of soln: r,z,u_re,u_im or
287 /// at n_plot^2 plot points
288 void output(FILE* file_pt, const unsigned& n_plot);
289
290 /// Output exact soln: r,z,u_re_exact,u_im_exact
291 /// at n_plot^2 plot points
292 void output_fct(std::ostream& outfile,
293 const unsigned& n_plot,
295
296 /// Output exact soln: (dummy time-dependent version to
297 /// keep intel compiler happy)
298 virtual void output_fct(
299 std::ostream& outfile,
300 const unsigned& n_plot,
301 const double& time,
303 {
304 throw OomphLibError("There is no time-dependent output_fct() for "
305 "PMLFourierDecomposedHelmholtz elements ",
306 OOMPH_CURRENT_FUNCTION,
307 OOMPH_EXCEPTION_LOCATION);
308 }
309
310
311 /// Output function for real part of full time-dependent fct
312 /// u = Re( (u_r +i u_i) exp(-i omega t)
313 /// at phase angle omega t = phi.
314 /// r,z,u at n_plot plot points in each coordinate
315 /// direction
316 void output_real_fct(std::ostream& outfile,
317 const double& phi,
318 const unsigned& n_plot,
320
321
322 /// Get error against and norm of exact solution
323 void compute_error(std::ostream& outfile,
325 double& error,
326 double& norm);
327
328
329 /// Dummy, time dependent error checker
330 void compute_error(std::ostream& outfile,
332 const double& time,
333 double& error,
334 double& norm)
335 {
336 throw OomphLibError("There is no time-dependent compute_error() for "
337 "PMLFourierDecomposedHelmholtz elements",
338 OOMPH_CURRENT_FUNCTION,
339 OOMPH_EXCEPTION_LOCATION);
340 }
341
342 /// Compute norm of fe solution
343 void compute_norm(double& norm);
344
345
346 /// Access function: Pointer to source function
348 {
349 return Source_fct_pt;
350 }
351
352
353 /// Access function: Pointer to source function. Const version
355 {
356 return Source_fct_pt;
357 }
358
359
360 /// Get source term at (Eulerian) position x. This function is
361 /// virtual to allow overloading in multi-physics problems where
362 /// the strength of the source function might be determined by
363 /// another system of equations.
365 const unsigned& ipt,
366 const Vector<double>& x,
367 std::complex<double>& source) const
368 {
369 // If no source function has been set, return zero
370 if (Source_fct_pt == 0)
371 {
372 source = std::complex<double>(0.0, 0.0);
373 }
374 else
375 {
376 // Get source strength
377 (*Source_fct_pt)(x, source);
378 }
379 }
380
381 /// Pure virtual function in which we specify the
382 /// values to be pinned (and set to zero) on the outer edge of
383 /// the pml layer. All of them! Vector is resized internally.
385 Vector<unsigned>& values_to_pin)
386 {
387 values_to_pin.resize(2);
388 for (unsigned j = 0; j < 2; j++)
389 {
390 values_to_pin[j] = j;
391 }
392 }
393
394
395 /// Get flux: flux[i] = du/dx_i for real and imag part
397 Vector<std::complex<double>>& flux) const
398 {
399 // Find out how many nodes there are in the element
400 const unsigned n_node = nnode();
401
402 // Set up memory for the shape and test functions
403 Shape psi(n_node);
404 DShape dpsidx(n_node, 2);
405
406 // Call the derivatives of the shape and test functions
407 dshape_eulerian(s, psi, dpsidx);
408
409 // Initialise to zero
410 const std::complex<double> zero(0.0, 0.0);
411 for (unsigned j = 0; j < 2; j++)
412 {
413 flux[j] = zero;
414 }
415
416 // Loop over nodes
417 for (unsigned l = 0; l < n_node; l++)
418 {
419 // Cache the complex value of the unknown
420 const std::complex<double> u_value(
421 this->nodal_value(l,
423 this->nodal_value(l,
425
426 // Loop over derivative directions
427 for (unsigned j = 0; j < 2; j++)
428 {
429 flux[j] += u_value * dpsidx(l, j);
430 }
431 }
432 }
433
434
435 /// Add the element's contribution to its residual vector (wrapper)
437 {
438 // Call the generic residuals function with flag set to 0
439 // using a dummy matrix argument
442 }
443
444
445 /// Add the element's contribution to its residual vector and
446 /// element Jacobian matrix (wrapper)
448 DenseMatrix<double>& jacobian)
449 {
450 // Call the generic routine with the flag set to 1
452 residuals, jacobian, 1);
453 }
454
455
456 /// Return FE representation of function value u(s)
457 /// at local coordinate s
459 const Vector<double>& s) const
460 {
461 // Find number of nodes
462 const unsigned n_node = nnode();
463
464 // Local shape function
465 Shape psi(n_node);
466
467 // Find values of shape function
468 shape(s, psi);
469
470 // Initialise value of u
471 std::complex<double> interpolated_u(0.0, 0.0);
472
473 // Get the index at which the helmholtz unknown is stored
474 const unsigned u_nodal_index_real =
476 const unsigned u_nodal_index_imag =
478
479 // Loop over the local nodes and sum
480 for (unsigned l = 0; l < n_node; l++)
481 {
482 // Make a temporary complex number from the stored data
483 const std::complex<double> u_value(
484 this->nodal_value(l, u_nodal_index_real),
485 this->nodal_value(l, u_nodal_index_imag));
486 // Add to the interpolated value
487 interpolated_u += u_value * psi[l];
488 }
489 return interpolated_u;
490 }
491
492
493 /// Self-test: Return 0 for OK
494 unsigned self_test();
495
496
497 protected:
498 /// Compute pml coefficients at position x and integration point ipt.
499 /// pml_laplace_factor is used in the residual contribution from the laplace
500 /// operator, similarly pml_k_squared_factor is used in the contribution
501 /// from the k^2 of the Helmholtz operator.
503 const unsigned& ipt,
504 const Vector<double>& x,
505 Vector<std::complex<double>>& pml_laplace_factor,
506 std::complex<double>& pml_k_squared_factor)
507 {
508 /// Vector which points from the inner boundary to x
509 Vector<double> nu(2);
510 for (unsigned k = 0; k < 2; k++)
511 {
512 nu[k] = x[k] - this->Pml_inner_boundary[k];
513 }
514
515 /// Vector which points from the inner boundary to the edge of the
516 /// boundary
517 Vector<double> pml_width(2);
518 for (unsigned k = 0; k < 2; k++)
519 {
520 pml_width[k] =
521 this->Pml_outer_boundary[k] - this->Pml_inner_boundary[k];
522 }
523
524 // Declare gamma_i vectors of complex numbers for PML weights
525 Vector<std::complex<double>> pml_gamma(2);
526
527 if (this->Pml_is_enabled)
528 {
529 // Cache k_squared to pass into mapping function
530 double k_squared_local = k_squared();
531
532 for (unsigned k = 0; k < 2; k++)
533 {
534 // If PML is enabled in the respective direction
535 if (this->Pml_direction_active[k])
536 {
538 nu[k], pml_width[k], k_squared_local);
539 }
540 else
541 {
542 pml_gamma[k] = 1.0;
543 }
544 }
545
546 /// for 2D, in order:
547 /// g_y/g_x, g_x/g_y for Laplace bit and g_x*g_y for Helmholtz bit
548 /// for 3D, in order: g_y*g_x/g_x, g*x*g_z/g_y, g_x*g_y/g_z for Laplace
549 /// bit and g_x*g_y*g_z for Helmholtz factor
550 pml_laplace_factor[0] = pml_gamma[1] / pml_gamma[0];
551 pml_laplace_factor[1] = pml_gamma[0] / pml_gamma[1];
552 pml_k_squared_factor = pml_gamma[0] * pml_gamma[1];
553 }
554 else
555 {
556 /// The weights all default to 1.0 as if the propagation
557 /// medium is the physical domain
558 for (unsigned k = 0; k < 2; k++)
559 {
560 pml_laplace_factor[k] = std::complex<double>(1.0, 0.0);
561 }
562
563 pml_k_squared_factor = std::complex<double>(1.0, 0.0);
564 }
565 }
566
567
568 /// Compute complex variable r at position x[0] and
569 /// integration point ipt
570 void compute_complex_r(const unsigned& ipt,
571 const Vector<double>& x,
572 std::complex<double>& complex_r)
573 {
574 // Cache current position r
575 double r = x[0];
576
577 /// The complex r variable is only imaginary on two
578 /// conditions: First, the decaying nature of the
579 /// pml layers is active. Secondly, the
580 /// integration point is contained in the right pml
581 /// layer or the two corner pml layers.
582
583 // If the complex r variable is imaginary
584 if (this->Pml_is_enabled && (this->Pml_direction_active[0]))
585 {
586 double nu = x[0] - Pml_inner_boundary[0];
587 double pml_width = Pml_outer_boundary[0] - Pml_inner_boundary[0];
588 double k_squared_local = k_squared();
589
590 // Determine the complex r variable
591 complex_r =
593 nu, pml_width, Pml_inner_boundary[0], k_squared_local);
594 }
595 else
596 {
597 // The complex r variable is infact purely real, and
598 // is equal to x[0]
599 complex_r = std::complex<double>(r, 0.0);
600 }
601
602 } // end of compute_complex_r
603
604 /// Return a pointer to the PML Mapping object
606 {
608 }
609
610 /// Return a pointer to the PML Mapping object (const version)
612 const
613 {
615 }
616
617 /// Static so that the class doesn't need to instantiate a new default
618 /// everytime it uses it
621
622
623 /// Shape/test functions and derivs w.r.t. to global coords at
624 /// local coord. s; return Jacobian of mapping
626 const Vector<double>& s,
627 Shape& psi,
628 DShape& dpsidx,
629 Shape& test,
630 DShape& dtestdx) const = 0;
631
632
633 /// Shape/test functions and derivs w.r.t. to global coords at
634 /// integration point ipt; return Jacobian of mapping
636 const unsigned& ipt,
637 Shape& psi,
638 DShape& dpsidx,
639 Shape& test,
640 DShape& dtestdx) const = 0;
641
642 /// Compute element residual Vector only (if flag=and/or element
643 /// Jacobian matrix
645 Vector<double>& residuals,
646 DenseMatrix<double>& jacobian,
647 const unsigned& flag);
648
649 /// Pointer to source function:
651
652 /// Pointer to k^2 (wavenumber squared)
654
655 /// Pointer to class which holds the pml mapping function (also known
656 /// as gamma) and the associated transformed coordinate
659
660 /// Pointer to wavenumber complex shift
661 double* Alpha_pt;
662
663 /// Static default value for the physical constants (initialised to zero)
665
666 /// Pointer to Fourier wave number
668 };
669
670
671 /// ////////////////////////////////////////////////////////////////////////
672 /// ////////////////////////////////////////////////////////////////////////
673 /// ////////////////////////////////////////////////////////////////////////
674
675
676 //======================================================================
677 /// QPMLFourierDecomposedHelmholtzElement elements are
678 /// linear/quadrilateral/brick-shaped PMLFourierDecomposedHelmholtz
679 /// elements with isoparametric interpolation for the function.
680 //======================================================================
681 template<unsigned NNODE_1D>
683 : public virtual QElement<2, NNODE_1D>,
685 {
686 private:
687 /// Static int that holds the number of variables at
688 /// nodes: always the same
689 static const unsigned Initial_Nvalue;
690
691 public:
692 /// Constructor: Call constructors for QElement and
693 /// PMLFourierDecomposedHelmholtz equations
696 {
697 }
698
699 /// Broken copy constructor
702
703 /// Broken assignment operator
704 /*void operator=(const
705 QPMLFourierDecomposedHelmholtzElement<NNODE_1D>&) = delete;*/
706
707
708 /// Required # of `values' (pinned or dofs)
709 /// at node n
710 inline unsigned required_nvalue(const unsigned& n) const
711 {
712 return Initial_Nvalue;
713 }
714
715 /// Output function: r,z,u
716 void output(std::ostream& outfile)
717 {
719 }
720
721 /// Output function:
722 /// r,z,u at n_plot^2 plot points
723 void output(std::ostream& outfile, const unsigned& n_plot)
724 {
726 }
727
728 /// Output function for real part of full time-dependent solution
729 /// u = Re( (u_r +i u_i) exp(-i omega t)
730 /// at phase angle omega t = phi.
731 /// r,z,u at n_plot plot points in each coordinate
732 /// direction
733 void output_real(std::ostream& outfile,
734 const double& phi,
735 const unsigned& n_plot)
736 {
738 }
739
740 /// C-style output function: r,z,u
741 void output(FILE* file_pt)
742 {
744 }
745
746 /// C-style output function:
747 /// r,z,u at n_plot^2 plot points
748 void output(FILE* file_pt, const unsigned& n_plot)
749 {
751 }
752
753 /// Output function for an exact solution:
754 /// r,z,u_exact at n_plot^2 plot points
755 void output_fct(std::ostream& outfile,
756 const unsigned& n_plot,
758 {
760 outfile, n_plot, exact_soln_pt);
761 }
762
763 /// Output function for real part of full time-dependent fct
764 /// u = Re( (u_r +i u_i) exp(-i omega t)
765 /// at phase angle omega t = phi.
766 /// r,z,u at n_plot plot points in each coordinate
767 /// direction
768 void output_real_fct(std::ostream& outfile,
769 const double& phi,
770 const unsigned& n_plot,
772 {
774 outfile, phi, n_plot, exact_soln_pt);
775 }
776
777
778 /// Output function for a time-dependent exact solution.
779 /// r,z,u_exact at n_plot^2 plot points
780 /// (Calls the steady version)
781 void output_fct(std::ostream& outfile,
782 const unsigned& n_plot,
783 const double& time,
785 {
787 outfile, n_plot, time, exact_soln_pt);
788 }
789
790 protected:
791 /// Shape, test functions & derivs. w.r.t. to global coords.
792 /// Return Jacobian.
794 const Vector<double>& s,
795 Shape& psi,
796 DShape& dpsidx,
797 Shape& test,
798 DShape& dtestdx) const;
799
800
801 /// Shape, test functions & derivs. w.r.t. to global coords. at
802 /// integration point ipt. Return Jacobian.
804 const unsigned& ipt,
805 Shape& psi,
806 DShape& dpsidx,
807 Shape& test,
808 DShape& dtestdx) const;
809 };
810
811
812 // Inline functions:
813
814
815 //======================================================================
816 /// Define the shape functions and test functions and derivatives
817 /// w.r.t. global coordinates and return Jacobian of mapping.
818 ///
819 /// Galerkin: Test functions = shape functions
820 //======================================================================
821 template<unsigned NNODE_1D>
824 const Vector<double>& s,
825 Shape& psi,
826 DShape& dpsidx,
827 Shape& test,
828 DShape& dtestdx) const
829 {
830 // Call the geometrical shape functions and derivatives
831 const double J = this->dshape_eulerian(s, psi, dpsidx);
832
833 // Set the test functions equal to the shape functions
834 test = psi;
835 dtestdx = dpsidx;
836
837 // Return the jacobian
838 return J;
839 }
840
841
842 //======================================================================
843 /// Define the shape functions and test functions and derivatives
844 /// w.r.t. global coordinates and return Jacobian of mapping.
845 ///
846 /// Galerkin: Test functions = shape functions
847 //======================================================================
848 template<unsigned NNODE_1D>
851 const unsigned& ipt,
852 Shape& psi,
853 DShape& dpsidx,
854 Shape& test,
855 DShape& dtestdx) const
856 {
857 // Call the geometrical shape functions and derivatives
858 const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
859
860 // Set the pointers of the test functions
861 test = psi;
862 dtestdx = dpsidx;
863
864 // Return the jacobian
865 return J;
866 }
867
868 /// /////////////////////////////////////////////////////////////////////
869 /// /////////////////////////////////////////////////////////////////////
870 /// /////////////////////////////////////////////////////////////////////
871
872
873 //=======================================================================
874 /// Face geometry for the QPMLFourierDecomposedHelmholtzElement
875 /// elements:
876 /// The spatial dimension of the face elements is one lower than that of the
877 /// bulk element but they have the same number of points
878 /// along their 1D edges.
879 //=======================================================================
880 template<unsigned NNODE_1D>
882 : public virtual QElement<1, NNODE_1D>
883 {
884 public:
885 /// Constructor: Call the constructor for the
886 /// appropriate lower-dimensional QElement
887 FaceGeometry() : QElement<1, NNODE_1D>() {}
888 };
889
890
891 /// /////////////////////////////////////////////////////////////////////
892 /// /////////////////////////////////////////////////////////////////////
893 /// /////////////////////////////////////////////////////////////////////
894
895
896 //==========================================================
897 /// Fourier decomposed Helmholtz upgraded to become projectable
898 //==========================================================
899 template<class FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENT>
901 : public virtual ProjectableElement<FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENT>
902 {
903 public:
904 /// Constructor [this was only required explicitly
905 /// from gcc 4.5.2 onwards...]
907
908 /// Specify the values associated with field fld.
909 /// The information is returned in a vector of pairs which comprise
910 /// the Data object and the value within it, that correspond to field fld.
912 {
913#ifdef PARANOID
914 if (fld > 1)
915 {
916 std::stringstream error_stream;
917 error_stream << "Fourier decomposed Helmholtz elements only store 2 "
918 "fields so fld = "
919 << fld << " is illegal \n";
920 throw OomphLibError(
921 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
922 }
923#endif
924
925 // Create the vector
926 unsigned nnod = this->nnode();
927 Vector<std::pair<Data*, unsigned>> data_values(nnod);
928
929 // Loop over all nodes
930 for (unsigned j = 0; j < nnod; j++)
931 {
932 // Add the data value associated field: The node itself
933 data_values[j] = std::make_pair(this->node_pt(j), fld);
934 }
935
936 // Return the vector
937 return data_values;
938 }
939
940 /// Number of fields to be projected: 2 (real and imag part)
942 {
943 return 2;
944 }
945
946 /// Number of history values to be stored for fld-th field.
947 /// (Note: count includes current value!)
948 unsigned nhistory_values_for_projection(const unsigned& fld)
949 {
950#ifdef PARANOID
951 if (fld > 1)
952 {
953 std::stringstream error_stream;
954 error_stream << "Helmholtz elements only store two fields so fld = "
955 << fld << " is illegal\n";
956 throw OomphLibError(
957 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
958 }
959#endif
960 return this->node_pt(0)->ntstorage();
961 }
962
963 /// Number of positional history values
964 /// (Note: count includes current value!)
966 {
967 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
968 }
969
970 /// Return Jacobian of mapping and shape functions of field fld
971 /// at local coordinate s
972 double jacobian_and_shape_of_field(const unsigned& fld,
973 const Vector<double>& s,
974 Shape& psi)
975 {
976#ifdef PARANOID
977 if (fld > 1)
978 {
979 std::stringstream error_stream;
980 error_stream << "Helmholtz elements only store two fields so fld = "
981 << fld << " is illegal.\n";
982 throw OomphLibError(
983 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
984 }
985#endif
986 unsigned n_dim = this->dim();
987 unsigned n_node = this->nnode();
988 Shape test(n_node);
989 DShape dpsidx(n_node, n_dim), dtestdx(n_node, n_dim);
990 double J =
991 this->dshape_and_dtest_eulerian_pml_fourier_decomposed_helmholtz(
992 s, psi, dpsidx, test, dtestdx);
993 return J;
994 }
995
996
997 /// Return interpolated field fld at local coordinate s, at time
998 /// level t (t=0: present; t>0: history values)
999 double get_field(const unsigned& t,
1000 const unsigned& fld,
1001 const Vector<double>& s)
1002 {
1003#ifdef PARANOID
1004 if (fld > 1)
1005 {
1006 std::stringstream error_stream;
1007 error_stream << "Helmholtz elements only store two fields so fld = "
1008 << fld << " is illegal\n";
1009 throw OomphLibError(
1010 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1011 }
1012#endif
1013 // Find the index at which the variable is stored
1014 std::complex<unsigned> complex_u_nodal_index =
1015 this->u_index_pml_fourier_decomposed_helmholtz();
1016 unsigned u_nodal_index = 0;
1017 if (fld == 0)
1018 {
1019 u_nodal_index = complex_u_nodal_index.real();
1020 }
1021 else
1022 {
1023 u_nodal_index = complex_u_nodal_index.imag();
1024 }
1025
1026
1027 // Local shape function
1028 unsigned n_node = this->nnode();
1029 Shape psi(n_node);
1030
1031 // Find values of shape function
1032 this->shape(s, psi);
1033
1034 // Initialise value of u
1035 double interpolated_u = 0.0;
1036
1037 // Sum over the local nodes
1038 for (unsigned l = 0; l < n_node; l++)
1039 {
1040 interpolated_u += this->nodal_value(t, l, u_nodal_index) * psi[l];
1041 }
1042 return interpolated_u;
1043 }
1044
1045
1046 /// Return number of values in field fld: One per node
1047 unsigned nvalue_of_field(const unsigned& fld)
1048 {
1049#ifdef PARANOID
1050 if (fld > 1)
1051 {
1052 std::stringstream error_stream;
1053 error_stream << "Helmholtz elements only store two fields so fld = "
1054 << fld << " is illegal\n";
1055 throw OomphLibError(
1056 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1057 }
1058#endif
1059 return this->nnode();
1060 }
1061
1062
1063 /// Return local equation number of value j in field fld.
1064 int local_equation(const unsigned& fld, const unsigned& j)
1065 {
1066#ifdef PARANOID
1067 if (fld > 1)
1068 {
1069 std::stringstream error_stream;
1070 error_stream << "Helmholtz elements only store two fields so fld = "
1071 << fld << " is illegal\n";
1072 throw OomphLibError(
1073 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1074 }
1075#endif
1076 std::complex<unsigned> complex_u_nodal_index =
1077 this->u_index_pml_fourier_decomposed_helmholtz();
1078 unsigned u_nodal_index = 0;
1079 if (fld == 0)
1080 {
1081 u_nodal_index = complex_u_nodal_index.real();
1082 }
1083 else
1084 {
1085 u_nodal_index = complex_u_nodal_index.imag();
1086 }
1087 return this->nodal_local_eqn(j, u_nodal_index);
1088 }
1089
1090
1091 /// Output FE representation of soln: x,y,u or x,y,z,u at
1092 /// n_plot^DIM plot points
1093 void output(std::ostream& outfile, const unsigned& nplot)
1094 {
1096 }
1097 };
1098
1099
1100 //=======================================================================
1101 /// Face geometry for element is the same as that for the underlying
1102 /// wrapped element
1103 //=======================================================================
1104 template<class ELEMENT>
1106 : public virtual FaceGeometry<ELEMENT>
1107 {
1108 public:
1109 FaceGeometry() : FaceGeometry<ELEMENT>() {}
1110 };
1111
1112
1113 //=======================================================================
1114 /// Face geometry of the Face Geometry for element is the same as
1115 /// that for the underlying wrapped element
1116 //=======================================================================
1117 template<class ELEMENT>
1120 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
1121 {
1122 public:
1124 };
1125
1126
1127 /// /////////////////////////////////////////////////////////////////////
1128 /// /////////////////////////////////////////////////////////////////////
1129 /// /////////////////////////////////////////////////////////////////////
1130
1131 //=======================================================================
1132 /// Policy class defining the elements to be used in the actual
1133 /// PML layers. Same!
1134 //=======================================================================
1135 template<unsigned NNODE_1D>
1137 : public virtual QPMLFourierDecomposedHelmholtzElement<NNODE_1D>
1138 {
1139 public:
1140 /// Constructor: Call the constructor for the
1141 /// appropriate QElement
1143 };
1144
1145} // namespace oomph
1146
1147#endif
static char t char * s
Definition: cfortran.h:568
char t
Definition: cfortran.h:568
The mapping function propsed by Bermudez et al, appears to be the best and so this will be the defaul...
std::complex< double > transformed_coordinate(const double &nu_i, const double &pml_width_i, const double &pml_inner_boundary, const double &k_squared)
Overwrite the pure PML mapping coefficient function to return the transformed coordinate proposed by ...
std::complex< double > gamma(const double &nu_i, const double &pml_width_i, const double &k_squared)
Overwrite the pure PML mapping coefficient function to return the mapping function proposed by Bermud...
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
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...
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
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
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
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
virtual void fill_in_generic_residual_contribution_pml_fourier_decomposed_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::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
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...
PMLMappingAndTransformedCoordinate *const & pml_mapping_and_transformed_coordinate_pt() const
Return a pointer to the PML Mapping object (const version)
void(* PMLFourierDecomposedHelmholtzSourceFctPt)(const Vector< double > &x, std::complex< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
void output(std::ostream &outfile)
Output with default number of plot points.
virtual double dshape_and_dtest_eulerian_pml_fourier_decomposed_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...
void compute_complex_r(const unsigned &ipt, const Vector< double > &x, std::complex< double > &complex_r)
Compute complex variable r at position x[0] and integration point ipt.
static double Default_Physical_Constant_Value
Static default value for the physical constants (initialised to zero)
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
std::complex< double > interpolated_u_pml_fourier_decomposed_helmholtz(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
virtual std::complex< unsigned > u_index_pml_fourier_decomposed_helmholtz() const
Broken assignment operator.
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)
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.
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)
PMLFourierDecomposedHelmholtzEquations(const PMLFourierDecomposedHelmholtzEquations &dummy)=delete
Broken copy constructor.
virtual void get_source_pml_fourier_decomposed_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...
PMLFourierDecomposedHelmholtzSourceFctPt Source_fct_pt
Pointer to source function:
PMLFourierDecomposedHelmholtzSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,z,u_re_exact,u_im_exact at n_plot^2 plot points.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
int *& pml_fourier_wavenumber_pt()
Get pointer to Fourier 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...
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...
void output(FILE *file_pt)
C_style output with default number of 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...
virtual double dshape_and_dtest_eulerian_at_knot_pml_fourier_decomposed_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 ...
PMLFourierDecomposedHelmholtzSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
static BermudezPMLMappingAndTransformedCoordinate Default_pml_mapping_and_transformed_coordinate
Static so that the class doesn't need to instantiate a new default everytime it uses it.
PMLMappingAndTransformedCoordinate * Pml_mapping_and_transformed_coordinate_pt
Pointer to class which holds the pml mapping function (also known as gamma) and the associated transf...
PMLMappingAndTransformedCoordinate *& pml_mapping_and_transformed_coordinate_pt()
Return a pointer to the PML Mapping object.
General definition of policy class defining the elements to be used in the actual PML layers....
Definition: pml_meshes.h:48
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual std::complex< double > transformed_coordinate(const double &nu_i, const double &pml_width_i, const double &pml_inner_boundary, const double &k_squared)=0
Pure virtual to return PML transformed coordinate, also known as as function of where where h is t...
virtual std::complex< double > gamma(const double &nu_i, const double &pml_width_i, const double &k_squared)=0
Pure virtual to return PML mapping gamma, where gamma is the as function of where where h is the v...
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
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...
ProjectablePMLFourierDecomposedHelmholtzElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
unsigned nfields_for_projection()
Number of fields to be projected: 2 (real and imag part)
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
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 ...
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.
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.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
double dshape_and_dtest_eulerian_pml_fourier_decomposed_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: r,z,u at n_plot^2 plot points.
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. r,z,u_exact at n_plot^2 plot points (Calls the s...
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: r,z,u_exact at n_plot^2 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...
QPMLFourierDecomposedHelmholtzElement(const QPMLFourierDecomposedHelmholtzElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
void output(FILE *file_pt)
C-style output function: r,z,u.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
double dshape_and_dtest_eulerian_at_knot_pml_fourier_decomposed_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....
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
QPMLFourierDecomposedHelmholtzElement()
Constructor: Call constructors for QElement and PMLFourierDecomposedHelmholtz equations.
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(std::ostream &outfile)
Output function: r,z,u.
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.
double plgndr1(const unsigned &n, const double &x)
Legendre polynomials depending on one parameter.
double plgndr2(const unsigned &l, const unsigned &m, const double &x)
Legendre polynomials depending on two parameters.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...