pml_fourier_decomposed_helmholtz_elements.cc
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// Non-inline functions for Helmholtz elements
28
29/// ///////////////////////////////////////////////////////////////////
30/// ///////////////////////////////////////////////////////////////////
31/// ///////////////////////////////////////////////////////////////////
32
33
34namespace oomph
35{
36 //========================================================================
37 /// Helper namespace for functions required for Helmholtz computations
38 //========================================================================
39 namespace Legendre_functions_helper
40 {
41 //========================================================================
42 // Factorial
43 //========================================================================
44 double factorial(const unsigned& l)
45 {
46 if (l == 0) return 1.0;
47 return double(l * factorial(l - 1));
48 }
49
50
51 //========================================================================
52 /// Legendre polynomials depending on one parameter
53 //========================================================================
54 double plgndr1(const unsigned& n, const double& x)
55 {
56 unsigned i;
57 double pmm, pmm1;
58 double pmm2 = 0;
59
60#ifdef PARANOID
61 // Shout if things went wrong
62 if (std::fabs(x) > 1.0)
63 {
64 std::ostringstream error_stream;
65 error_stream << "Bad arguments in routine plgndr1: x=" << x
66 << " but should be less than 1 in absolute value.\n";
67 throw OomphLibError(
68 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
69 }
70#endif
71
72 // Compute pmm : if l=m it's finished
73 pmm = 1.0;
74 if (n == 0)
75 {
76 return pmm;
77 }
78
79 pmm1 = x * 1.0;
80 if (n == 1)
81 {
82 return pmm1;
83 }
84 else
85 {
86 for (i = 2; i <= n; i++)
87 {
88 pmm2 = (x * (2 * i - 1) * pmm1 - (i - 1) * pmm) / i;
89 pmm = pmm1;
90 pmm1 = pmm2;
91 }
92 return pmm2;
93 }
94
95 } // end of plgndr1
96
97
98 //========================================================================
99 // Legendre polynomials depending on two parameters
100 //========================================================================
101 double plgndr2(const unsigned& l, const unsigned& m, const double& x)
102 {
103 unsigned i, ll;
104 double fact, pmm, pmmp1, somx2;
105 double pll = 0.0;
106
107#ifdef PARANOID
108 // Shout if things went wrong
109 if (std::fabs(x) > 1.0)
110 {
111 std::ostringstream error_stream;
112 error_stream << "Bad arguments in routine plgndr2: x=" << x
113 << " but should be less than 1 in absolute value.\n";
114 throw OomphLibError(
115 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
116 }
117#endif
118
119 // This one is easy...
120 if (m > l)
121 {
122 return 0.0;
123 }
124
125 // Compute pmm : if l=m it's finished
126 pmm = 1.0;
127 if (m > 0)
128 {
129 somx2 = sqrt((1.0 - x) * (1.0 + x));
130 fact = 1.0;
131 for (i = 1; i <= m; i++)
132 {
133 pmm *= -fact * somx2;
134 fact += 2.0;
135 }
136 }
137 if (l == m) return pmm;
138
139 // Compute pmmp1 : if l=m+1 it's finished
140 else
141 {
142 pmmp1 = x * (2 * m + 1) * pmm;
143 if (l == (m + 1))
144 {
145 return pmmp1;
146 }
147 // Compute pll : if l>m+1 it's finished
148 else
149 {
150 for (ll = m + 2; ll <= l; ll++)
151 {
152 pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
153 pmm = pmmp1;
154 pmmp1 = pll;
155 }
156 return pll;
157 }
158 }
159 } // end of plgndr2
160
161 } // namespace Legendre_functions_helper
162
163
164 /// ////////////////////////////////////////////////////////////////////
165 /// ////////////////////////////////////////////////////////////////////
166 /// ////////////////////////////////////////////////////////////////////
167
168
169 //======================================================================
170 /// Set the data for the number of Variables at each node, always two
171 /// (real and imag part) in every case
172 //======================================================================
173 template<unsigned NNODE_1D>
174 const unsigned
176
177 /// PML Helmholtz equations static data, so that by default we can point to a
178 /// 0
179 double
181 0.0;
182
183 //======================================================================
184 /// Compute element residual Vector and/or element Jacobian matrix
185 ///
186 /// flag=1: compute both
187 /// flag=0: compute only residual Vector
188 ///
189 /// Pure version without hanging nodes
190 //======================================================================
193 Vector<double>& residuals,
194 DenseMatrix<double>& jacobian,
195 const unsigned& flag)
196 {
197 // Find out how many nodes there are
198 const unsigned n_node = nnode();
199
200 // Set up memory for the shape and test functions
201 Shape psi(n_node), test(n_node);
202 DShape dpsidx(n_node, 2), dtestdx(n_node, 2);
203
204 // Set the value of n_intpt
205 const unsigned n_intpt = integral_pt()->nweight();
206
207 // Integers to store the local equation and unknown numbers
208 int local_eqn_real = 0, local_unknown_real = 0;
209 int local_eqn_imag = 0, local_unknown_imag = 0;
210
211 // Loop over the integration points
212 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
213 {
214 // Get the integral weight
215 double w = integral_pt()->weight(ipt);
216
217 // Call the derivatives of the shape and test functions
218 double J =
220 ipt, psi, dpsidx, test, dtestdx);
221
222 // Premultiply the weights and the Jacobian
223 double W = w * J;
224
225 // Calculate local values of unknown
226 // Allocate and initialise to zero
227 std::complex<double> interpolated_u(0.0, 0.0);
229 Vector<std::complex<double>> interpolated_dudx(2);
230
231 // Calculate function value and derivatives:
232 //-----------------------------------------
233 // Loop over nodes
234 for (unsigned l = 0; l < n_node; l++)
235 {
236 // Loop over directions
237 for (unsigned j = 0; j < 2; j++)
238 {
239 interpolated_x[j] += raw_nodal_position(l, j) * psi(l);
240 }
241
242 // Get the nodal value of the helmholtz unknown
243 const std::complex<double> u_value(
247
248 // Add to the interpolated value
249 interpolated_u += u_value * psi(l);
250
251 // Loop over directions
252 for (unsigned j = 0; j < 2; j++)
253 {
254 interpolated_dudx[j] += u_value * dpsidx(l, j);
255 }
256 }
257
258 // Get source function
259 //-------------------
260 std::complex<double> source(0.0, 0.0);
262
263 double n = (double)pml_fourier_wavenumber();
264 double n_squared = n * n;
265
266
267 // Declare a vector of complex numbers for pml weights on the Laplace bit
268 Vector<std::complex<double>> pml_laplace_factor(2);
269 // Declare a complex number for pml weights on the mass matrix bit
270 std::complex<double> pml_k_squared_factor =
271 std::complex<double>(1.0, 0.0);
272
273 // All the PML weights that participate in the assemby process
274 // are computed here. pml_laplace_factor will contain the entries
275 // for the Laplace bit, while pml_k_squared_factor contains the
276 // contributions to the Helmholtz bit. Both default to 1.0, should the PML
277 // not be enabled via enable_pml.
279 ipt, interpolated_x, pml_laplace_factor, pml_k_squared_factor);
280
281
282 // Determine the complex r variable. The variable is
283 // only complex once it enters the right pml domain or either
284 // of the two corner pml domains, otherwise it acts like the
285 // variable r.
286 std::complex<double> complex_r = std::complex<double>(1.0, 0.0);
287 compute_complex_r(ipt, interpolated_x, complex_r);
288
289 // Calculate Jacobian
290 // ------------------
291 std::complex<double> pml_k_squared_jacobian =
292 -pml_k_squared_factor *
293 (k_squared() - n_squared / (complex_r * complex_r)) * complex_r * W;
294
295 // Term from the Laplace operator
296 Vector<std::complex<double>> pml_laplace_jacobian(2);
297 for (unsigned k = 0; k < 2; k++)
298 {
299 pml_laplace_jacobian[k] = pml_laplace_factor[k] * complex_r * W;
300 }
301
302 // Calculate residuals
303 // -------------------
304 // Note: it is a linear equation so residual = jacobian * u
305 std::complex<double> pml_k_squared_residual =
306 pml_k_squared_jacobian * interpolated_u;
307
308 // Term from the Laplace operator
309 Vector<std::complex<double>> pml_laplace_residual(2);
310 for (unsigned k = 0; k < 2; k++)
311 {
312 pml_laplace_residual[k] =
313 pml_laplace_jacobian[k] * interpolated_dudx[k];
314 }
315
316 // Assemble residuals and Jacobian
317 //--------------------------------
318 // Loop over the test functions
319 for (unsigned l = 0; l < n_node; l++)
320 {
321 // Get the equation numbers
322 local_eqn_real =
324 local_eqn_imag =
326
327 // first, add the real contributions
328 //-------------------------------------------
329
330 /*IF it's not a boundary condition*/
331 if (local_eqn_real >= 0)
332 {
333 // Add k squared term of equation
334 residuals[local_eqn_real] += pml_k_squared_residual.real() * test(l);
335
336 // Add the term from the Laplace operator
337 for (unsigned k = 0; k < 2; k++)
338 {
339 residuals[local_eqn_real] +=
340 pml_laplace_residual[k].real() * dtestdx(l, k);
341 }
342
343 // Add contributions to the jacobian
344 //----------------------------------
345 if (flag)
346 {
347 // Loop over the velocity shape functions again
348 for (unsigned l2 = 0; l2 < n_node; l2++)
349 {
350 // Get the unknown numbers
351 local_unknown_real = nodal_local_eqn(
353 local_unknown_imag = nodal_local_eqn(
355
356 // If at a non-zero degree of freedom add in the entry
357 if (local_unknown_real >= 0)
358 {
359 // Add the helmholtz contribution to the jacobian
360 jacobian(local_eqn_real, local_unknown_real) +=
361 pml_k_squared_jacobian.real() * psi(l2) * test(l);
362
363 // Add the term from the Laplace operator to the jacobian
364 for (unsigned k = 0; k < 2; k++)
365 {
366 jacobian(local_eqn_real, local_unknown_real) +=
367 pml_laplace_jacobian[k].real() * dpsidx(l2, k) *
368 dtestdx(l, k);
369 }
370 }
371 // If at a non-zero degree of freedom add in the entry
372 if (local_unknown_imag >= 0)
373 {
374 // Add k squared term to jacobian
375 jacobian(local_eqn_real, local_unknown_imag) +=
376 -pml_k_squared_jacobian.imag() * psi(l2) * test(l);
377
378 // Add contribution to elemental Matrix
379 for (unsigned k = 0; k < 2; k++)
380 {
381 jacobian(local_eqn_real, local_unknown_imag) +=
382 -pml_laplace_jacobian[k].imag() * dpsidx(l2, k) *
383 dtestdx(l, k);
384 }
385 }
386
387 } // end of loop over velocity shape functions
388 } // end of if(flag)
389 }
390
391 // Second, add the imaginary contribution
392 //------------------------------------------------
393
394 // IF it's not a boundary condition
395 if (local_eqn_imag >= 0)
396 {
397 // Add k squared term of equation
398 residuals[local_eqn_imag] += pml_k_squared_residual.imag() * test(l);
399
400 // Add the term from the Laplace operator
401 for (unsigned k = 0; k < 2; k++)
402 {
403 residuals[local_eqn_imag] +=
404 pml_laplace_residual[k].imag() * dtestdx(l, k);
405 }
406
407 // Add the contribution to the jacobian
408 //-----------------------
409 if (flag)
410 {
411 // Loop over the velocity shape functions again
412 for (unsigned l2 = 0; l2 < n_node; l2++)
413 {
414 local_unknown_imag = nodal_local_eqn(
416 local_unknown_real = nodal_local_eqn(
418
419 // If at a non-zero degree of freedom add in the entry
420 if (local_unknown_imag >= 0)
421 {
422 // Add the helmholtz contribution
423 jacobian(local_eqn_imag, local_unknown_imag) +=
424 pml_k_squared_jacobian.real() * psi(l2) * test(l);
425
426 // Add contribution to elemental Matrix
427 for (unsigned k = 0; k < 2; k++)
428 {
429 jacobian(local_eqn_imag, local_unknown_imag) +=
430 pml_laplace_jacobian[k].real() * dpsidx(l2, k) *
431 dtestdx(l, k);
432 }
433 }
434 // If at a non-zero degree of freedom add in the entry
435 if (local_unknown_real >= 0)
436 {
437 // Add the helmholtz contribution
438 jacobian(local_eqn_imag, local_unknown_real) +=
439 pml_k_squared_jacobian.imag() * psi(l2) * test(l);
440
441 // Add contribution to elemental Matrix
442 for (unsigned k = 0; k < 2; k++)
443 {
444 jacobian(local_eqn_imag, local_unknown_real) +=
445 pml_laplace_jacobian[k].imag() * dpsidx(l2, k) *
446 dtestdx(l, k);
447 }
448 }
449 }
450 }
451 }
452 }
453 } // End of loop over integration points
454 }
455
456
457 //======================================================================
458 /// Self-test: Return 0 for OK
459 //======================================================================
461 {
462 bool passed = true;
463
464 // Check lower-level stuff
465 if (FiniteElement::self_test() != 0)
466 {
467 passed = false;
468 }
469
470 // Return verdict
471 if (passed)
472 {
473 return 0;
474 }
475 else
476 {
477 return 1;
478 }
479 }
480
481
482 //======================================================================
483 /// Output function:
484 ///
485 /// r,z,u_re,u_imag
486 ///
487 /// nplot points in each coordinate direction
488 //======================================================================
490 const unsigned& nplot)
491 {
492 // Vector of local coordinates
493 Vector<double> s(2);
494
495 // Tecplot header info
496 outfile << tecplot_zone_string(nplot);
497
498 // Loop over plot points
499 unsigned num_plot_points = nplot_points(nplot);
500 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
501 {
502 // Get local coordinates of plot point
503 get_s_plot(iplot, nplot, s);
504 std::complex<double> u(
506 for (unsigned i = 0; i < 2; i++)
507 {
508 outfile << interpolated_x(s, i) << " ";
509 }
510 outfile << u.real() << " " << u.imag() << std::endl;
511 }
512
513 // Write tecplot footer (e.g. FE connectivity lists)
514 write_tecplot_zone_footer(outfile, nplot);
515 }
516
517
518 //======================================================================
519 /// Output function for real part of full time-dependent solution
520 ///
521 /// u = Re( (u_r +i u_i) exp(-i omega t)
522 ///
523 /// at phase angle omega t = phi.
524 ///
525 /// r,z,u
526 ///
527 /// Output at nplot points in each coordinate direction
528 //======================================================================
530 std::ostream& outfile, const double& phi, const unsigned& nplot)
531 {
532 // Vector of local coordinates
533 Vector<double> s(2);
534
535 // Tecplot header info
536 outfile << tecplot_zone_string(nplot);
537
538 // Loop over plot points
539 unsigned num_plot_points = nplot_points(nplot);
540 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
541 {
542 // Get local coordinates of plot point
543 get_s_plot(iplot, nplot, s);
544 std::complex<double> u(
546 for (unsigned i = 0; i < 2; i++)
547 {
548 outfile << interpolated_x(s, i) << " ";
549 }
550 outfile << u.real() * cos(phi) + u.imag() * sin(phi) << std::endl;
551 }
552
553 // Write tecplot footer (e.g. FE connectivity lists)
554 write_tecplot_zone_footer(outfile, nplot);
555 }
556
557
558 //======================================================================
559 /// C-style output function:
560 ///
561 /// r,z,u
562 ///
563 /// nplot points in each coordinate direction
564 //======================================================================
566 const unsigned& nplot)
567 {
568 // Vector of local coordinates
569 Vector<double> s(2);
570
571 // Tecplot header info
572 fprintf(file_pt, "%s", tecplot_zone_string(nplot).c_str());
573
574 // Loop over plot points
575 unsigned num_plot_points = nplot_points(nplot);
576 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
577 {
578 // Get local coordinates of plot point
579 get_s_plot(iplot, nplot, s);
580 std::complex<double> u(
582
583 for (unsigned i = 0; i < 2; i++)
584 {
585 fprintf(file_pt, "%g ", interpolated_x(s, i));
586 }
587
588 for (unsigned i = 0; i < 2; i++)
589 {
590 fprintf(file_pt, "%g ", interpolated_x(s, i));
591 }
592 fprintf(file_pt, "%g ", u.real());
593 fprintf(file_pt, "%g \n", u.imag());
594 }
595
596 // Write tecplot footer (e.g. FE connectivity lists)
597 write_tecplot_zone_footer(file_pt, nplot);
598 }
599
600
601 //======================================================================
602 /// Output exact solution
603 ///
604 /// Solution is provided via function pointer.
605 /// Plot at a given number of plot points.
606 ///
607 /// r,z,u_exact
608 //======================================================================
610 std::ostream& outfile,
611 const unsigned& nplot,
613 {
614 // Vector of local coordinates
615 Vector<double> s(2);
616
617 // Vector for coordintes
618 Vector<double> x(2);
619
620 // Tecplot header info
621 outfile << tecplot_zone_string(nplot);
622
623 // Exact solution Vector
624 Vector<double> exact_soln(2);
625
626 // Loop over plot points
627 unsigned num_plot_points = nplot_points(nplot);
628 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
629 {
630 // Get local coordinates of plot point
631 get_s_plot(iplot, nplot, s);
632
633 // Get x position as Vector
634 interpolated_x(s, x);
635
636 // Get exact solution at this point
637 (*exact_soln_pt)(x, exact_soln);
638
639 // Output r,z,u_exact
640 for (unsigned i = 0; i < 2; i++)
641 {
642 outfile << x[i] << " ";
643 }
644 outfile << exact_soln[0] << " " << exact_soln[1] << std::endl;
645 }
646
647 // Write tecplot footer (e.g. FE connectivity lists)
648 write_tecplot_zone_footer(outfile, nplot);
649 }
650
651
652 //======================================================================
653 /// Output function for real part of full time-dependent fct
654 ///
655 /// u = Re( (u_r +i u_i) exp(-i omega t)
656 ///
657 /// at phase angle omega t = phi.
658 ///
659 /// r,z,u
660 ///
661 /// Output at nplot points in each coordinate direction
662 //======================================================================
664 std::ostream& outfile,
665 const double& phi,
666 const unsigned& nplot,
668 {
669 // Vector of local coordinates
670 Vector<double> s(2);
671
672 // Vector for coordintes
673 Vector<double> x(2);
674
675 // Tecplot header info
676 outfile << tecplot_zone_string(nplot);
677
678 // Exact solution Vector
679 Vector<double> exact_soln(2);
680
681 // Loop over plot points
682 unsigned num_plot_points = nplot_points(nplot);
683 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
684 {
685 // Get local coordinates of plot point
686 get_s_plot(iplot, nplot, s);
687
688 // Get x position as Vector
689 interpolated_x(s, x);
690
691 // Get exact solution at this point
692 (*exact_soln_pt)(x, exact_soln);
693
694 // Output x,y,...,u_exact
695 for (unsigned i = 0; i < 2; i++)
696 {
697 outfile << x[i] << " ";
698 }
699 outfile << exact_soln[0] * cos(phi) + exact_soln[1] * sin(phi)
700 << std::endl;
701 }
702
703 // Write tecplot footer (e.g. FE connectivity lists)
704 write_tecplot_zone_footer(outfile, nplot);
705 }
706
707
708 //======================================================================
709 /// Validate against exact solution
710 ///
711 /// Solution is provided via function pointer.
712 /// Plot error at a given number of plot points.
713 ///
714 //======================================================================
716 std::ostream& outfile,
718 double& error,
719 double& norm)
720 {
721 // Initialise
722 error = 0.0;
723 norm = 0.0;
724
725 // Vector of local coordinates
726 Vector<double> s(2);
727
728 // Vector for coordintes
729 Vector<double> x(2);
730
731 // Find out how many nodes there are in the element
732 unsigned n_node = nnode();
733
734 Shape psi(n_node);
735
736 // Set the value of n_intpt
737 unsigned n_intpt = integral_pt()->nweight();
738
739 // Tecplot
740 outfile << "ZONE" << std::endl;
741
742 // Exact solution Vector
743 Vector<double> exact_soln(2);
744
745 // Loop over the integration points
746 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
747 {
748 // Assign values of s
749 for (unsigned i = 0; i < 2; i++)
750 {
751 s[i] = integral_pt()->knot(ipt, i);
752 }
753
754 // Get the integral weight
755 double w = integral_pt()->weight(ipt);
756
757 // Get jacobian of mapping
758 double J = J_eulerian(s);
759
760 // Premultiply the weights and the Jacobian
761 double W = w * J;
762
763 // Get x position as Vector
764 interpolated_x(s, x);
765
766 // Get FE function value
767 std::complex<double> u_fe =
769
770 // Get exact solution at this point
771 (*exact_soln_pt)(x, exact_soln);
772
773 // Output r,z,error
774 for (unsigned i = 0; i < 2; i++)
775 {
776 outfile << x[i] << " ";
777 }
778 outfile << exact_soln[0] << " " << exact_soln[1] << " "
779 << exact_soln[0] - u_fe.real() << " "
780 << exact_soln[1] - u_fe.imag() << std::endl;
781
782 // Add to error and norm
783 norm +=
784 (exact_soln[0] * exact_soln[0] + exact_soln[1] * exact_soln[1]) * W;
785 error += ((exact_soln[0] - u_fe.real()) * (exact_soln[0] - u_fe.real()) +
786 (exact_soln[1] - u_fe.imag()) * (exact_soln[1] - u_fe.imag())) *
787 W;
788 }
789 }
790
791
792 //======================================================================
793 /// Compute norm of fe solution
794 //======================================================================
796 {
797 // Initialise
798 norm = 0.0;
799
800 // Vector of local coordinates
801 Vector<double> s(2);
802
803 // Vector for coordintes
804 Vector<double> x(2);
805
806 // Find out how many nodes there are in the element
807 unsigned n_node = nnode();
808
809 Shape psi(n_node);
810
811 // Set the value of n_intpt
812 unsigned n_intpt = integral_pt()->nweight();
813
814 // Loop over the integration points
815 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
816 {
817 // Assign values of s
818 for (unsigned i = 0; i < 2; i++)
819 {
820 s[i] = integral_pt()->knot(ipt, i);
821 }
822
823 // Get the integral weight
824 double w = integral_pt()->weight(ipt);
825
826 // Get jacobian of mapping
827 double J = J_eulerian(s);
828
829 // Premultiply the weights and the Jacobian
830 double W = w * J;
831
832 // Get FE function value
833 std::complex<double> u_fe =
835
836 // Add to norm
837 norm += (u_fe.real() * u_fe.real() + u_fe.imag() * u_fe.imag()) * W;
838 }
839 }
840
841
842 //====================================================================
843 // Force build of templates
844 //====================================================================
848
852
853} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
The mapping function propsed by Bermudez et al, appears to be the best and so this will be the defaul...
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition: elements.cc:4103
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned 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
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2576
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation....
Definition: elements.cc:1686
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
Definition: elements.cc:4440
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
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...
void output(std::ostream &outfile)
Output with default number of plot points.
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)
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.
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...
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 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 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 ...
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.
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...