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