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