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