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