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-2023 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
27 #include "pml_helmholtz_elements.h"
28 
29 
30 namespace 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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...