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_FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENTS_HEADER
28#define OOMPH_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
46namespace oomph
47{
48 //========================================================================
49 /// Helper namespace for functions required for Helmholtz computations
50 //========================================================================
51 namespace Legendre_functions_helper
52 {
53 /// Factorial
54 extern double factorial(const unsigned& l);
55
56 /// Legendre polynomials depending on one parameter
57 extern double plgndr1(const unsigned& n, const double& x);
58
59 /// Legendre polynomials depending on two parameters
60 extern double plgndr2(const unsigned& l,
61 const unsigned& m,
62 const double& x);
63
64 } // namespace Legendre_functions_helper
65
66
67 /// ////////////////////////////////////////////////////////////////////
68 /// ////////////////////////////////////////////////////////////////////
69 /// ////////////////////////////////////////////////////////////////////
70
71
72 //=============================================================
73 /// A class for all isoparametric elements that solve the
74 /// Helmholtz equations.
75 /// \f[ \nabla^2 U + k^2 U = f \f]
76 /// in Fourier decomposed form (cylindrical polars):
77 /// \f[ U(r,\varphi,z) = \Re( u^{(n)}(r,z) \exp(-i n \varphi)) \f]
78 /// We are solving for \f$ u^{(n)}(r,z)\f$ for given parameters
79 /// \f$ k^2 \f$ and \f$ n \f$ .
80 /// This contains the generic maths. Shape functions, geometric
81 /// mapping etc. must get implemented in derived class.
82 //=============================================================
84 {
85 public:
86 /// Function pointer to source function fct(x,f(x)) --
87 /// x is a Vector!
89 const Vector<double>& x, std::complex<double>& f);
90
91
92 /// Constructor
95 {
96 }
97
98 /// Broken copy constructor
100 const FourierDecomposedHelmholtzEquations& dummy) = delete;
101
102 /// Broken assignment operator
103 // Commented out broken assignment operator because this can lead to a
104 // conflict warning when used in the virtual inheritence hierarchy.
105 // Essentially the compiler doesn't realise that two separate
106 // implementations of the broken function are the same and so, quite
107 // rightly, it shouts.
108 /*void operator=(const FourierDecomposedHelmholtzEquations&) = delete;*/
109
110
111 /// Return the index at which the unknown value
112 /// is stored: Real/imag part of index contains (real) index of
113 /// real/imag part.
114 virtual inline std::complex<unsigned> u_index_fourier_decomposed_helmholtz()
115 const
116 {
117 return std::complex<unsigned>(0, 1);
118 }
119
120
121 /// Get pointer to square of wavenumber
122 double*& k_squared_pt()
123 {
124 return K_squared_pt;
125 }
126
127 /// Get the square of wavenumber
128 double k_squared()
129 {
130 if (K_squared_pt == 0)
131 {
132 return 0.0;
133 }
134 else
135 {
136 return *K_squared_pt;
137 }
138 }
139
140 /// Get pointer to Fourier wavenumber
142 {
143 return N_fourier_pt;
144 }
145
146 /// Get the Fourier wavenumber
148 {
149 if (N_fourier_pt == 0)
150 {
151 return 0;
152 }
153 else
154 {
155 return *N_fourier_pt;
156 }
157 }
158
159 /// Output with default number of plot points
160 void output(std::ostream& outfile)
161 {
162 const unsigned n_plot = 5;
163 output(outfile, n_plot);
164 }
165
166 /// Output FE representation of soln: x,y,u_re,u_im or
167 /// x,y,z,u_re,u_im at n_plot^2 plot points
168 void output(std::ostream& outfile, const unsigned& n_plot);
169
170 /// Output function for real part of full time-dependent solution
171 /// u = Re( (u_r +i u_i) exp(-i omega t)
172 /// at phase angle omega t = phi.
173 /// r,z,u at n_plot plot points in each coordinate
174 /// direction
175 void output_real(std::ostream& outfile,
176 const double& phi,
177 const unsigned& n_plot);
178
179 /// C_style output with default number of plot points
180 void output(FILE* file_pt)
181 {
182 const unsigned n_plot = 5;
183 output(file_pt, n_plot);
184 }
185
186 /// C-style output FE representation of soln: r,z,u_re,u_im or
187 /// at n_plot^2 plot points
188 void output(FILE* file_pt, const unsigned& n_plot);
189
190 /// Output exact soln: r,z,u_re_exact,u_im_exact
191 /// at n_plot^2 plot points
192 void output_fct(std::ostream& outfile,
193 const unsigned& n_plot,
195
196 /// Output exact soln: (dummy time-dependent version to
197 /// keep intel compiler happy)
198 virtual void output_fct(
199 std::ostream& outfile,
200 const unsigned& n_plot,
201 const double& time,
203 {
204 throw OomphLibError("There is no time-dependent output_fct() for "
205 "FourierDecomposedHelmholtz elements ",
206 OOMPH_CURRENT_FUNCTION,
207 OOMPH_EXCEPTION_LOCATION);
208 }
209
210
211 /// Output function for real part of full time-dependent fct
212 /// u = Re( (u_r +i u_i) exp(-i omega t)
213 /// at phase angle omega t = phi.
214 /// r,z,u at n_plot plot points in each coordinate
215 /// direction
216 void output_real_fct(std::ostream& outfile,
217 const double& phi,
218 const unsigned& n_plot,
220
221
222 /// Get error against and norm of exact solution
223 void compute_error(std::ostream& outfile,
225 double& error,
226 double& norm);
227
228
229 /// Dummy, time dependent error checker
230 void compute_error(std::ostream& outfile,
232 const double& time,
233 double& error,
234 double& norm)
235 {
236 throw OomphLibError("There is no time-dependent compute_error() for "
237 "FourierDecomposedHelmholtz elements",
238 OOMPH_CURRENT_FUNCTION,
239 OOMPH_EXCEPTION_LOCATION);
240 }
241
242 /// Compute norm of fe solution
243 void compute_norm(double& norm);
244
245 /// Access function: Pointer to source function
247 {
248 return Source_fct_pt;
249 }
250
251 /// Access function: Pointer to source function. Const version
253 {
254 return Source_fct_pt;
255 }
256
257 /// Get source term at (Eulerian) position x. This function is
258 /// virtual to allow overloading in multi-physics problems where
259 /// the strength of the source function might be determined by
260 /// another system of equations.
262 const unsigned& ipt,
263 const Vector<double>& x,
264 std::complex<double>& source) const
265 {
266 // If no source function has been set, return zero
267 if (Source_fct_pt == 0)
268 {
269 source = std::complex<double>(0.0, 0.0);
270 }
271 else
272 {
273 // Get source strength
274 (*Source_fct_pt)(x, source);
275 }
276 }
277
278
279 /// Get flux: flux[i] = du/dx_i for real and imag part
281 Vector<std::complex<double>>& flux) const
282 {
283 // Find out how many nodes there are in the element
284 const unsigned n_node = nnode();
285
286 // Set up memory for the shape and test functions
287 Shape psi(n_node);
288 DShape dpsidx(n_node, 2);
289
290 // Call the derivatives of the shape and test functions
291 dshape_eulerian(s, psi, dpsidx);
292
293 // Initialise to zero
294 const std::complex<double> zero(0.0, 0.0);
295 for (unsigned j = 0; j < 2; j++)
296 {
297 flux[j] = zero;
298 }
299
300 // Loop over nodes
301 for (unsigned l = 0; l < n_node; l++)
302 {
303 // Cache the complex value of the unknown
304 const std::complex<double> u_value(
307
308 // Loop over derivative directions
309 for (unsigned j = 0; j < 2; j++)
310 {
311 flux[j] += u_value * dpsidx(l, j);
312 }
313 }
314 }
315
316
317 /// Add the element's contribution to its residual vector (wrapper)
319 {
320 // Call the generic residuals function with flag set to 0
321 // using a dummy matrix argument
324 }
325
326
327 /// Add the element's contribution to its residual vector and
328 /// element Jacobian matrix (wrapper)
330 DenseMatrix<double>& jacobian)
331 {
332 // Call the generic routine with the flag set to 1
334 residuals, jacobian, 1);
335 }
336
337
338 /// Return FE representation of function value u(s)
339 /// at local coordinate s
341 const Vector<double>& s) const
342 {
343 // Find number of nodes
344 const unsigned n_node = nnode();
345
346 // Local shape function
347 Shape psi(n_node);
348
349 // Find values of shape function
350 shape(s, psi);
351
352 // Initialise value of u
353 std::complex<double> interpolated_u(0.0, 0.0);
354
355 // Get the index at which the helmholtz unknown is stored
356 const unsigned u_nodal_index_real =
358 const unsigned u_nodal_index_imag =
360
361 // Loop over the local nodes and sum
362 for (unsigned l = 0; l < n_node; l++)
363 {
364 // Make a temporary complex number from the stored data
365 const std::complex<double> u_value(
366 this->nodal_value(l, u_nodal_index_real),
367 this->nodal_value(l, u_nodal_index_imag));
368 // Add to the interpolated value
369 interpolated_u += u_value * psi[l];
370 }
371 return interpolated_u;
372 }
373
374
375 /// Self-test: Return 0 for OK
376 unsigned self_test();
377
378
379 protected:
380 /// Shape/test functions and derivs w.r.t. to global coords at
381 /// local coord. s; return Jacobian of mapping
383 const Vector<double>& s,
384 Shape& psi,
385 DShape& dpsidx,
386 Shape& test,
387 DShape& dtestdx) const = 0;
388
389
390 /// Shape/test functions and derivs w.r.t. to global coords at
391 /// integration point ipt; return Jacobian of mapping
393 const unsigned& ipt,
394 Shape& psi,
395 DShape& dpsidx,
396 Shape& test,
397 DShape& dtestdx) const = 0;
398
399 /// Compute element residual Vector only (if flag=and/or element
400 /// Jacobian matrix
402 Vector<double>& residuals,
403 DenseMatrix<double>& jacobian,
404 const unsigned& flag);
405
406 /// Pointer to source function:
408
409 /// Pointer to square of wavenumber
411
412 /// Pointer to Fourier wave number
414 };
415
416
417 /// ////////////////////////////////////////////////////////////////////////
418 /// ////////////////////////////////////////////////////////////////////////
419 /// ////////////////////////////////////////////////////////////////////////
420
421
422 //======================================================================
423 /// QFourierDecomposedHelmholtzElement elements are
424 /// linear/quadrilateral/brick-shaped FourierDecomposedHelmholtz
425 /// elements with isoparametric interpolation for the function.
426 //======================================================================
427 template<unsigned NNODE_1D>
429 : public virtual QElement<2, NNODE_1D>,
431 {
432 private:
433 /// Static int that holds the number of variables at
434 /// nodes: always the same
435 static const unsigned Initial_Nvalue;
436
437 public:
438 /// Constructor: Call constructors for QElement and
439 /// FourierDecomposedHelmholtz equations
442 {
443 }
444
445 /// Broken copy constructor
448
449 /// Broken assignment operator
450 /*void operator=(const QFourierDecomposedHelmholtzElement<NNODE_1D>&) =
451 * delete;*/
452
453
454 /// Required # of `values' (pinned or dofs)
455 /// at node n
456 inline unsigned required_nvalue(const unsigned& n) const
457 {
458 return Initial_Nvalue;
459 }
460
461 /// Output function: r,z,u
462 void output(std::ostream& outfile)
463 {
465 }
466
467 /// Output function:
468 /// r,z,u at n_plot^2 plot points
469 void output(std::ostream& outfile, const unsigned& n_plot)
470 {
472 }
473
474 /// Output function for real part of full time-dependent solution
475 /// u = Re( (u_r +i u_i) exp(-i omega t)
476 /// at phase angle omega t = phi.
477 /// r,z,u at n_plot plot points in each coordinate
478 /// direction
479 void output_real(std::ostream& outfile,
480 const double& phi,
481 const unsigned& n_plot)
482 {
484 }
485
486 /// C-style output function: r,z,u
487 void output(FILE* file_pt)
488 {
490 }
491
492 /// C-style output function:
493 /// r,z,u at n_plot^2 plot points
494 void output(FILE* file_pt, const unsigned& n_plot)
495 {
497 }
498
499 /// Output function for an exact solution:
500 /// r,z,u_exact at n_plot^2 plot points
501 void output_fct(std::ostream& outfile,
502 const unsigned& n_plot,
504 {
506 outfile, n_plot, exact_soln_pt);
507 }
508
509 /// Output function for real part of full time-dependent fct
510 /// u = Re( (u_r +i u_i) exp(-i omega t)
511 /// at phase angle omega t = phi.
512 /// r,z,u at n_plot plot points in each coordinate
513 /// direction
514 void output_real_fct(std::ostream& outfile,
515 const double& phi,
516 const unsigned& n_plot,
518 {
520 outfile, phi, n_plot, exact_soln_pt);
521 }
522
523
524 /// Output function for a time-dependent exact solution.
525 /// r,z,u_exact at n_plot^2 plot points
526 /// (Calls the steady version)
527 void output_fct(std::ostream& outfile,
528 const unsigned& n_plot,
529 const double& time,
531 {
533 outfile, n_plot, time, exact_soln_pt);
534 }
535
536 protected:
537 /// Shape, test functions & derivs. w.r.t. to global coords.
538 /// Return Jacobian.
540 const Vector<double>& s,
541 Shape& psi,
542 DShape& dpsidx,
543 Shape& test,
544 DShape& dtestdx) const;
545
546
547 /// Shape, test functions & derivs. w.r.t. to global coords. at
548 /// integration point ipt. Return Jacobian.
550 const unsigned& ipt,
551 Shape& psi,
552 DShape& dpsidx,
553 Shape& test,
554 DShape& dtestdx) const;
555 };
556
557
558 // Inline functions:
559
560
561 //======================================================================
562 /// Define the shape functions and test functions and derivatives
563 /// w.r.t. global coordinates and return Jacobian of mapping.
564 ///
565 /// Galerkin: Test functions = shape functions
566 //======================================================================
567 template<unsigned NNODE_1D>
570 const Vector<double>& s,
571 Shape& psi,
572 DShape& dpsidx,
573 Shape& test,
574 DShape& dtestdx) const
575 {
576 // Call the geometrical shape functions and derivatives
577 const double J = this->dshape_eulerian(s, psi, dpsidx);
578
579 // Set the test functions equal to the shape functions
580 test = psi;
581 dtestdx = dpsidx;
582
583 // Return the jacobian
584 return J;
585 }
586
587
588 //======================================================================
589 /// Define the shape functions and test functions and derivatives
590 /// w.r.t. global coordinates and return Jacobian of mapping.
591 ///
592 /// Galerkin: Test functions = shape functions
593 //======================================================================
594 template<unsigned NNODE_1D>
597 const unsigned& ipt,
598 Shape& psi,
599 DShape& dpsidx,
600 Shape& test,
601 DShape& dtestdx) const
602 {
603 // Call the geometrical shape functions and derivatives
604 const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
605
606 // Set the pointers of the test functions
607 test = psi;
608 dtestdx = dpsidx;
609
610 // Return the jacobian
611 return J;
612 }
613
614 /// /////////////////////////////////////////////////////////////////////
615 /// /////////////////////////////////////////////////////////////////////
616 /// /////////////////////////////////////////////////////////////////////
617
618
619 //=======================================================================
620 /// Face geometry for the QFourierDecomposedHelmholtzElement elements:
621 /// The spatial dimension of the face elements is one lower than that of the
622 /// bulk element but they have the same number of points
623 /// along their 1D edges.
624 //=======================================================================
625 template<unsigned NNODE_1D>
627 : public virtual QElement<1, NNODE_1D>
628 {
629 public:
630 /// Constructor: Call the constructor for the
631 /// appropriate lower-dimensional QElement
632 FaceGeometry() : QElement<1, NNODE_1D>() {}
633 };
634
635
636 /// /////////////////////////////////////////////////////////////////////
637 /// /////////////////////////////////////////////////////////////////////
638 /// /////////////////////////////////////////////////////////////////////
639
640
641 //==========================================================
642 /// Fourier decomposed Helmholtz upgraded to become projectable
643 //==========================================================
644 template<class FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENT>
646 : public virtual ProjectableElement<FOURIER_DECOMPOSED_HELMHOLTZ_ELEMENT>
647 {
648 public:
649 /// Constructor [this was only required explicitly
650 /// from gcc 4.5.2 onwards...]
652
653 /// Specify the values associated with field fld.
654 /// The information is returned in a vector of pairs which comprise
655 /// the Data object and the value within it, that correspond to field fld.
657 {
658#ifdef PARANOID
659 if (fld > 1)
660 {
661 std::stringstream error_stream;
662 error_stream << "Fourier decomposed Helmholtz elements only store 2 "
663 "fields so fld = "
664 << fld << " is illegal \n";
665 throw OomphLibError(
666 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
667 }
668#endif
669
670 // Create the vector
671 unsigned nnod = this->nnode();
672 Vector<std::pair<Data*, unsigned>> data_values(nnod);
673
674 // Loop over all nodes
675 for (unsigned j = 0; j < nnod; j++)
676 {
677 // Add the data value associated field: The node itself
678 data_values[j] = std::make_pair(this->node_pt(j), fld);
679 }
680
681 // Return the vector
682 return data_values;
683 }
684
685 /// Number of fields to be projected: 2 (real and imag part)
687 {
688 return 2;
689 }
690
691 /// Number of history values to be stored for fld-th field.
692 /// (Note: count includes current value!)
693 unsigned nhistory_values_for_projection(const unsigned& fld)
694 {
695#ifdef PARANOID
696 if (fld > 1)
697 {
698 std::stringstream error_stream;
699 error_stream << "Helmholtz elements only store two fields so fld = "
700 << fld << " is illegal\n";
701 throw OomphLibError(
702 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
703 }
704#endif
705 return this->node_pt(0)->ntstorage();
706 }
707
708 /// Number of positional history values
709 /// (Note: count includes current value!)
711 {
712 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
713 }
714
715 /// Return Jacobian of mapping and shape functions of field fld
716 /// at local coordinate s
717 double jacobian_and_shape_of_field(const unsigned& fld,
718 const Vector<double>& s,
719 Shape& psi)
720 {
721#ifdef PARANOID
722 if (fld > 1)
723 {
724 std::stringstream error_stream;
725 error_stream << "Helmholtz elements only store two fields so fld = "
726 << fld << " is illegal.\n";
727 throw OomphLibError(
728 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
729 }
730#endif
731 unsigned n_dim = this->dim();
732 unsigned n_node = this->nnode();
733 Shape test(n_node);
734 DShape dpsidx(n_node, n_dim), dtestdx(n_node, n_dim);
735 double J = this->dshape_and_dtest_eulerian_fourier_decomposed_helmholtz(
736 s, psi, dpsidx, test, dtestdx);
737 return J;
738 }
739
740
741 /// Return interpolated field fld at local coordinate s, at time
742 /// level t (t=0: present; t>0: history values)
743 double get_field(const unsigned& t,
744 const unsigned& fld,
745 const Vector<double>& s)
746 {
747#ifdef PARANOID
748 if (fld > 1)
749 {
750 std::stringstream error_stream;
751 error_stream << "Helmholtz elements only store two fields so fld = "
752 << fld << " is illegal\n";
753 throw OomphLibError(
754 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
755 }
756#endif
757 // Find the index at which the variable is stored
758 std::complex<unsigned> complex_u_nodal_index =
759 this->u_index_fourier_decomposed_helmholtz();
760 unsigned u_nodal_index = 0;
761 if (fld == 0)
762 {
763 u_nodal_index = complex_u_nodal_index.real();
764 }
765 else
766 {
767 u_nodal_index = complex_u_nodal_index.imag();
768 }
769
770
771 // Local shape function
772 unsigned n_node = this->nnode();
773 Shape psi(n_node);
774
775 // Find values of shape function
776 this->shape(s, psi);
777
778 // Initialise value of u
779 double interpolated_u = 0.0;
780
781 // Sum over the local nodes
782 for (unsigned l = 0; l < n_node; l++)
783 {
784 interpolated_u += this->nodal_value(t, l, u_nodal_index) * psi[l];
785 }
786 return interpolated_u;
787 }
788
789
790 /// Return number of values in field fld: One per node
791 unsigned nvalue_of_field(const unsigned& fld)
792 {
793#ifdef PARANOID
794 if (fld > 1)
795 {
796 std::stringstream error_stream;
797 error_stream << "Helmholtz elements only store two fields so fld = "
798 << fld << " is illegal\n";
799 throw OomphLibError(
800 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
801 }
802#endif
803 return this->nnode();
804 }
805
806
807 /// Return local equation number of value j in field fld.
808 int local_equation(const unsigned& fld, const unsigned& j)
809 {
810#ifdef PARANOID
811 if (fld > 1)
812 {
813 std::stringstream error_stream;
814 error_stream << "Helmholtz elements only store two fields so fld = "
815 << fld << " is illegal\n";
816 throw OomphLibError(
817 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
818 }
819#endif
820 std::complex<unsigned> complex_u_nodal_index =
821 this->u_index_fourier_decomposed_helmholtz();
822 unsigned u_nodal_index = 0;
823 if (fld == 0)
824 {
825 u_nodal_index = complex_u_nodal_index.real();
826 }
827 else
828 {
829 u_nodal_index = complex_u_nodal_index.imag();
830 }
831 return this->nodal_local_eqn(j, u_nodal_index);
832 }
833
834
835 /// Output FE representation of soln: x,y,u or x,y,z,u at
836 /// n_plot^DIM plot points
837 void output(std::ostream& outfile, const unsigned& nplot)
838 {
840 }
841 };
842
843
844 //=======================================================================
845 /// Face geometry for element is the same as that for the underlying
846 /// wrapped element
847 //=======================================================================
848 template<class ELEMENT>
850 : public virtual FaceGeometry<ELEMENT>
851 {
852 public:
853 FaceGeometry() : FaceGeometry<ELEMENT>() {}
854 };
855
856
857 //=======================================================================
858 /// Face geometry of the Face Geometry for element is the same as
859 /// that for the underlying wrapped element
860 //=======================================================================
861 template<class ELEMENT>
864 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
865 {
866 public:
868 };
869
870
871} // namespace oomph
872
873#endif
static char t char * s
Definition: cfortran.h:568
char t
Definition: cfortran.h:568
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
FourierDecomposedHelmholtzSourceFctPt Source_fct_pt
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...
virtual void fill_in_generic_residual_contribution_fourier_decomposed_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Broken assignment operator.
std::complex< double > interpolated_u_fourier_decomposed_helmholtz(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
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)
void(* FourierDecomposedHelmholtzSourceFctPt)(const Vector< double > &x, std::complex< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
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 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 compute_norm(double &norm)
Compute norm of fe solution.
FourierDecomposedHelmholtzSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
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.
virtual double dshape_and_dtest_eulerian_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...
virtual double dshape_and_dtest_eulerian_at_knot_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 ...
FourierDecomposedHelmholtzSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
FourierDecomposedHelmholtzEquations(const FourierDecomposedHelmholtzEquations &dummy)=delete
Broken copy constructor.
void output(std::ostream &outfile)
Output with default number of plot points.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
int *& fourier_wavenumber_pt()
Get pointer to Fourier wavenumber.
double *& k_squared_pt()
Get pointer to square of wavenumber.
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 get_source_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...
void output(FILE *file_pt)
C_style output with default number of plot points.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
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....
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
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...
unsigned nfields_for_projection()
Number of fields to be projected: 2 (real and imag part)
ProjectableFourierDecomposedHelmholtzElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
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 nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
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.
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_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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_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...
void output(FILE *file_pt)
C-style output function: r,z,u.
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(std::ostream &outfile)
Output function: r,z,u.
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...
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
QFourierDecomposedHelmholtzElement()
Constructor: Call constructors for QElement and FourierDecomposedHelmholtz equations.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,z,u at n_plot^2 plot points.
double dshape_and_dtest_eulerian_at_knot_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....
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,z,u at n_plot^2 plot points.
double dshape_and_dtest_eulerian_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.
QFourierDecomposedHelmholtzElement(const QFourierDecomposedHelmholtzElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...