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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 namespace oomph
35 {
36 
37 
38 //========================================================================
39 /// Helper namespace for functions required for Helmholtz computations
40 //========================================================================
41  namespace Legendre_functions_helper
42  {
43 
44 //========================================================================
45 // Factorial
46 //========================================================================
47  double factorial(const unsigned& l)
48  {
49  if(l==0) return 1.0;
50  return double(l*factorial(l-1));
51  }
52 
53 
54 //========================================================================
55 /// Legendre polynomials depending on one parameter
56 //========================================================================
57  double plgndr1(const unsigned& n, const double& x)
58  {
59  unsigned i;
60  double pmm,pmm1;
61  double pmm2=0;
62 
63 #ifdef PARANOID
64  // Shout if things went wrong
65  if (std::fabs(x) > 1.0)
66  {
67  std::ostringstream error_stream;
68  error_stream << "Bad arguments in routine plgndr1: x=" << x
69  << " but should be less than 1 in absolute value.\n";
70  throw OomphLibError(error_stream.str(),
71  OOMPH_CURRENT_FUNCTION,
72  OOMPH_EXCEPTION_LOCATION);
73  }
74 #endif
75 
76  // Compute pmm : if l=m it's finished
77  pmm=1.0;
78  if (n == 0)
79  {
80  return pmm;
81  }
82 
83  pmm1=x*1.0;
84  if (n == 1)
85  {
86  return pmm1;
87  }
88  else
89  {
90  for (i=2;i<=n;i++)
91  {
92  pmm2=(x*(2*i-1)*pmm1-(i-1)*pmm)/i;
93  pmm=pmm1;
94  pmm1=pmm2;
95  }
96  return pmm2;
97  }
98 
99  }//end of plgndr1
100 
101 
102 
103 //========================================================================
104 // Legendre polynomials depending on two parameters
105 //========================================================================
106  double plgndr2(const unsigned& l, const unsigned& m, const double& x)
107  {
108  unsigned i,ll;
109  double fact,pmm,pmmp1,somx2;
110  double pll=0.0;
111 
112 #ifdef PARANOID
113  // Shout if things went wrong
114  if (std::fabs(x) > 1.0)
115  {
116  std::ostringstream error_stream;
117  error_stream << "Bad arguments in routine plgndr2: x=" << x
118  << " but should be less than 1 in absolute value.\n";
119  throw OomphLibError(error_stream.str(),
120  OOMPH_CURRENT_FUNCTION,
121  OOMPH_EXCEPTION_LOCATION);
122  }
123 #endif
124 
125  // This one is easy...
126  if (m > l)
127  {
128  return 0.0;
129  }
130 
131  // Compute pmm : if l=m it's finished
132  pmm=1.0;
133  if (m > 0)
134  {
135  somx2=sqrt((1.0-x)*(1.0+x));
136  fact=1.0;
137  for (i=1;i<=m;i++)
138  {
139  pmm *= -fact*somx2;
140  fact += 2.0;
141  }
142  }
143  if (l == m) return pmm;
144 
145  // Compute pmmp1 : if l=m+1 it's finished
146  else
147  {
148  pmmp1=x*(2*m+1)*pmm;
149  if (l == (m+1))
150  {
151  return pmmp1;
152  }
153  // Compute pll : if l>m+1 it's finished
154  else
155  {
156  for (ll=m+2;ll<=l;ll++)
157  {
158  pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
159  pmm=pmmp1;
160  pmmp1=pll;
161  }
162  return pll;
163  }
164  }
165  }//end of plgndr2
166 
167  } // end namespace
168 
169 
170 
171 ///////////////////////////////////////////////////////////////////////
172 ///////////////////////////////////////////////////////////////////////
173 ///////////////////////////////////////////////////////////////////////
174 
175 
176 //======================================================================
177 /// Set the data for the number of Variables at each node, always two
178 /// (real and imag part) in every case
179 //======================================================================
180  template<unsigned NNODE_1D>
181  const unsigned QFourierDecomposedHelmholtzElement<NNODE_1D>::Initial_Nvalue=2;
182 
183 
184 //======================================================================
185 /// Compute element residual Vector and/or element Jacobian matrix
186 ///
187 /// flag=1: compute both
188 /// flag=0: compute only residual Vector
189 ///
190 /// Pure version without hanging nodes
191 //======================================================================
194  Vector<double> &residuals,
195  DenseMatrix<double> &jacobian,
196  const unsigned& flag)
197  {
198  //Find out how many nodes there are
199  const unsigned n_node = nnode();
200 
201  //Set up memory for the shape and test functions
202  Shape psi(n_node), test(n_node);
203  DShape dpsidx(n_node,2), dtestdx(n_node,2);
204 
205  //Set the value of n_intpt
206  const unsigned n_intpt = integral_pt()->nweight();
207 
208  //Integers to store the local equation and unknown numbers
209  int local_eqn_real=0, local_unknown_real=0;
210  int local_eqn_imag=0, local_unknown_imag=0;
211 
212  //Loop over the integration points
213  for(unsigned ipt=0;ipt<n_intpt;ipt++)
214  {
215  //Get the integral weight
216  double w = integral_pt()->weight(ipt);
217 
218  //Call the derivatives of the shape and test functions
220  (ipt,psi,dpsidx,test,dtestdx);
221 
222  //Premultiply the weights and the Jacobian
223  double W = w*J;
224 
225  //Calculate local values of unknown
226  //Allocate and initialise to zero
227  std::complex<double> interpolated_u(0.0,0.0);
229  Vector< std::complex<double> > interpolated_dudx(2);
230 
231  //Calculate function value and derivatives:
232  //-----------------------------------------
233  // Loop over nodes
234  for(unsigned l=0;l<n_node;l++)
235  {
236  // Loop over directions
237  for(unsigned j=0;j<2;j++)
238  {
239  interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
240  }
241 
242  //Get the nodal value of the helmholtz unknown
243  const std::complex<double> u_value(
246 
247  //Add to the interpolated value
248  interpolated_u += u_value*psi(l);
249 
250  // Loop over directions
251  for(unsigned j=0;j<2;j++)
252  {
253  interpolated_dudx[j] += u_value*dpsidx(l,j);
254  }
255  }
256 
257  //Get source function
258  //-------------------
259  std::complex<double> source(0.0,0.0);
260  get_source_fourier_decomposed_helmholtz(ipt,interpolated_x,source);
261 
262  double r = interpolated_x[0];
263  double rr = r*r;
264  double n = (double)fourier_wavenumber();
265  double n_squared = n*n;
266 
267  // Assemble residuals and Jacobian
268  //--------------------------------
269 
270  // Loop over the test functions
271  for(unsigned l=0;l<n_node;l++)
272  {
273 
274  // first, compute the real part contribution
275  //-------------------------------------------
276 
277  //Get the local equation
278  local_eqn_real =
280 
281  /*IF it's not a boundary condition*/
282  if(local_eqn_real >= 0)
283  {
284  // Add body force/source term and Helmholtz bit
285 
286  residuals[local_eqn_real] +=
287  (source.real()-((-n_squared/rr)+k_squared())*
288  interpolated_u.real())*test(l)*r*W;
289 
290  // The Helmholtz bit itself
291  for(unsigned k=0;k<2;k++)
292  {
293  residuals[local_eqn_real] +=
294  interpolated_dudx[k].real()*dtestdx(l,k)*r*W;
295  }
296 
297  // Calculate the jacobian
298  //-----------------------
299  if(flag)
300  {
301  //Loop over the velocity shape functions again
302  for(unsigned l2=0;l2<n_node;l2++)
303  {
304  local_unknown_real =
306 
307  //If at a non-zero degree of freedom add in the entry
308  if(local_unknown_real >= 0)
309  {
310  //Add contribution to elemental Matrix
311  for(unsigned i=0;i<2;i++)
312  {
313  jacobian(local_eqn_real,local_unknown_real)
314  += dpsidx(l2,i)*dtestdx(l,i)*r*W;
315  }
316 
317  // Add the helmholtz contribution
318  jacobian(local_eqn_real,local_unknown_real)
319  -= ((-n_squared/rr)+k_squared())*psi(l2)*test(l)*r*W;
320 
321  }//end of local_unknown
322  }
323  }
324  }
325 
326  // Second, compute the imaginary part contribution
327  //------------------------------------------------
328 
329  //Get the local equation
330  local_eqn_imag = nodal_local_eqn
332 
333  /*IF it's not a boundary condition*/
334  if(local_eqn_imag >= 0)
335  {
336  // Add body force/source term and Helmholtz bit
337  residuals[local_eqn_imag] +=
338  (source.imag() - ( (-n_squared/rr) + k_squared() ) *
339  interpolated_u.imag() )*test(l)*r*W;
340 
341  // The Helmholtz bit itself
342  for(unsigned k=0;k<2;k++)
343  {
344  residuals[local_eqn_imag] += interpolated_dudx[k].imag()*
345  dtestdx(l,k)*r*W;
346  }
347 
348  // Calculate the jacobian
349  //-----------------------
350  if(flag)
351  {
352  //Loop over the velocity shape functions again
353  for(unsigned l2=0;l2<n_node;l2++)
354  {
355  local_unknown_imag = nodal_local_eqn
357 
358  //If at a non-zero degree of freedom add in the entry
359  if(local_unknown_imag >= 0)
360  {
361  //Add contribution to Elemental Matrix
362  for(unsigned i=0;i<2;i++)
363  {
364  jacobian(local_eqn_imag,local_unknown_imag)
365  += dpsidx(l2,i)*dtestdx(l,i)*r*W;
366  }
367  // Add the helmholtz contribution
368  jacobian(local_eqn_imag,local_unknown_imag)
369  -= ((-n_squared/rr)+k_squared())*psi(l2)*test(l)*r*W;
370  }
371  }
372  }
373  }
374 
375  }
376  } // End of loop over integration points
377  }
378 
379 
380 //======================================================================
381 /// Self-test: Return 0 for OK
382 //======================================================================
384  {
385 
386  bool passed=true;
387 
388  // Check lower-level stuff
389  if (FiniteElement::self_test()!=0)
390  {
391  passed=false;
392  }
393 
394  // Return verdict
395  if (passed)
396  {
397  return 0;
398  }
399  else
400  {
401  return 1;
402  }
403  }
404 
405 
406 
407 //======================================================================
408 /// Output function:
409 ///
410 /// r,z,u_re,u_imag
411 ///
412 /// nplot points in each coordinate direction
413 //======================================================================
414  void FourierDecomposedHelmholtzEquations::output(std::ostream &outfile,
415  const unsigned &nplot)
416  {
417 
418  //Vector of local coordinates
419  Vector<double> s(2);
420 
421  // Tecplot header info
422  outfile << tecplot_zone_string(nplot);
423 
424  // Loop over plot points
425  unsigned num_plot_points=nplot_points(nplot);
426  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
427  {
428 
429  // Get local coordinates of plot point
430  get_s_plot(iplot,nplot,s);
431  std::complex<double> u(interpolated_u_fourier_decomposed_helmholtz(s));
432  for(unsigned i=0;i<2;i++)
433  {
434  outfile << interpolated_x(s,i) << " ";
435  }
436  outfile << u.real() << " " << u.imag() << std::endl;
437 
438  }
439 
440  // Write tecplot footer (e.g. FE connectivity lists)
441  write_tecplot_zone_footer(outfile,nplot);
442 
443  }
444 
445 
446 
447 
448 //======================================================================
449 /// Output function for real part of full time-dependent solution
450 ///
451 /// u = Re( (u_r +i u_i) exp(-i omega t)
452 ///
453 /// at phase angle omega t = phi.
454 ///
455 /// r,z,u
456 ///
457 /// Output at nplot points in each coordinate direction
458 //======================================================================
460  const double& phi,
461  const unsigned &nplot)
462  {
463 
464  //Vector of local coordinates
465  Vector<double> s(2);
466 
467  // Tecplot header info
468  outfile << tecplot_zone_string(nplot);
469 
470  // Loop over plot points
471  unsigned num_plot_points=nplot_points(nplot);
472  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
473  {
474 
475  // Get local coordinates of plot point
476  get_s_plot(iplot,nplot,s);
477  std::complex<double> u(interpolated_u_fourier_decomposed_helmholtz(s));
478  for(unsigned i=0;i<2;i++)
479  {
480  outfile << interpolated_x(s,i) << " ";
481  }
482  outfile << u.real()*cos(phi)+u.imag()*sin(phi) << std::endl;
483 
484  }
485 
486  // Write tecplot footer (e.g. FE connectivity lists)
487  write_tecplot_zone_footer(outfile,nplot);
488 
489  }
490 
491 
492 //======================================================================
493 /// C-style output function:
494 ///
495 /// r,z,u
496 ///
497 /// nplot points in each coordinate direction
498 //======================================================================
500  const unsigned &nplot)
501  {
502  //Vector of local coordinates
503  Vector<double> s(2);
504 
505  // Tecplot header info
506  fprintf(file_pt,"%s",tecplot_zone_string(nplot).c_str());
507 
508  // Loop over plot points
509  unsigned num_plot_points=nplot_points(nplot);
510  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
511  {
512  // Get local coordinates of plot point
513  get_s_plot(iplot,nplot,s);
514  std::complex<double> u(interpolated_u_fourier_decomposed_helmholtz(s));
515 
516  for(unsigned i=0;i<2;i++)
517  {
518  fprintf(file_pt,"%g ",interpolated_x(s,i));
519  }
520 
521  for(unsigned i=0;i<2;i++)
522  {
523  fprintf(file_pt,"%g ",interpolated_x(s,i));
524  }
525  fprintf(file_pt,"%g ",u.real());
526  fprintf(file_pt,"%g \n",u.imag());
527 
528  }
529 
530  // Write tecplot footer (e.g. FE connectivity lists)
531  write_tecplot_zone_footer(file_pt,nplot);
532  }
533 
534 
535 
536 //======================================================================
537  /// Output exact solution
538  ///
539  /// Solution is provided via function pointer.
540  /// Plot at a given number of plot points.
541  ///
542  /// r,z,u_exact
543 //======================================================================
545  std::ostream &outfile,
546  const unsigned &nplot,
548  {
549  //Vector of local coordinates
550  Vector<double> s(2);
551 
552  // Vector for coordintes
553  Vector<double> x(2);
554 
555  // Tecplot header info
556  outfile << tecplot_zone_string(nplot);
557 
558  // Exact solution Vector
559  Vector<double> exact_soln(2);
560 
561  // Loop over plot points
562  unsigned num_plot_points=nplot_points(nplot);
563  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
564  {
565  // Get local coordinates of plot point
566  get_s_plot(iplot,nplot,s);
567 
568  // Get x position as Vector
569  interpolated_x(s,x);
570 
571  // Get exact solution at this point
572  (*exact_soln_pt)(x,exact_soln);
573 
574  //Output r,z,u_exact
575  for(unsigned i=0;i<2;i++)
576  {
577  outfile << x[i] << " ";
578  }
579  outfile << exact_soln[0] << " " << exact_soln[1] << std::endl;
580  }
581 
582  // Write tecplot footer (e.g. FE connectivity lists)
583  write_tecplot_zone_footer(outfile,nplot);
584  }
585 
586 
587 
588 //======================================================================
589 /// Output function for real part of full time-dependent fct
590 ///
591 /// u = Re( (u_r +i u_i) exp(-i omega t)
592 ///
593 /// at phase angle omega t = phi.
594 ///
595 /// r,z,u
596 ///
597 /// Output at nplot points in each coordinate direction
598 //======================================================================
600  std::ostream &outfile,
601  const double& phi,
602  const unsigned &nplot,
604  {
605  //Vector of local coordinates
606  Vector<double> s(2);
607 
608  // Vector for coordintes
609  Vector<double> x(2);
610 
611  // Tecplot header info
612  outfile << tecplot_zone_string(nplot);
613 
614  // Exact solution Vector
615  Vector<double> exact_soln(2);
616 
617  // Loop over plot points
618  unsigned num_plot_points=nplot_points(nplot);
619  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
620  {
621  // Get local coordinates of plot point
622  get_s_plot(iplot,nplot,s);
623 
624  // Get x position as Vector
625  interpolated_x(s,x);
626 
627  // Get exact solution at this point
628  (*exact_soln_pt)(x,exact_soln);
629 
630  //Output x,y,...,u_exact
631  for(unsigned i=0;i<2;i++)
632  {
633  outfile << x[i] << " ";
634  }
635  outfile << exact_soln[0]*cos(phi)+ exact_soln[1]*sin(phi) << std::endl;
636  }
637 
638  // Write tecplot footer (e.g. FE connectivity lists)
639  write_tecplot_zone_footer(outfile,nplot);
640  }
641 
642 
643 
644 
645 //======================================================================
646  /// Validate against exact solution
647  ///
648  /// Solution is provided via function pointer.
649  /// Plot error at a given number of plot points.
650  ///
651 //======================================================================
653  std::ostream &outfile,
655  double& error, double& norm)
656  {
657 
658  // Initialise
659  error=0.0;
660  norm=0.0;
661 
662  //Vector of local coordinates
663  Vector<double> s(2);
664 
665  // Vector for coordintes
666  Vector<double> x(2);
667 
668  //Find out how many nodes there are in the element
669  unsigned n_node = nnode();
670 
671  Shape psi(n_node);
672 
673  //Set the value of n_intpt
674  unsigned n_intpt = integral_pt()->nweight();
675 
676  // Tecplot
677  outfile << "ZONE" << std::endl;
678 
679  // Exact solution Vector
680  Vector<double> exact_soln(2);
681 
682  //Loop over the integration points
683  for(unsigned ipt=0;ipt<n_intpt;ipt++)
684  {
685 
686  //Assign values of s
687  for(unsigned i=0;i<2;i++)
688  {
689  s[i] = integral_pt()->knot(ipt,i);
690  }
691 
692  //Get the integral weight
693  double w = integral_pt()->weight(ipt);
694 
695  // Get jacobian of mapping
696  double J=J_eulerian(s);
697 
698  //Premultiply the weights and the Jacobian
699  double W = w*J;
700 
701  // Get x position as Vector
702  interpolated_x(s,x);
703 
704  // Get FE function value
705  std::complex<double> u_fe=interpolated_u_fourier_decomposed_helmholtz(s);
706 
707  // Get exact solution at this point
708  (*exact_soln_pt)(x,exact_soln);
709 
710  //Output r,z,error
711  for(unsigned i=0;i<2;i++)
712  {
713  outfile << x[i] << " ";
714  }
715  outfile << exact_soln[0] << " " << exact_soln[1] << " "
716  << exact_soln[0]-u_fe.real() << " " << exact_soln[1]-u_fe.imag()
717  << std::endl;
718 
719  // Add to error and norm
720  norm+=(exact_soln[0]*exact_soln[0]+exact_soln[1]*exact_soln[1])*W;
721  error+=( (exact_soln[0]-u_fe.real())*(exact_soln[0]-u_fe.real())+
722  (exact_soln[1]-u_fe.imag())*(exact_soln[1]-u_fe.imag()) )*W;
723 
724  }
725  }
726 
727 
728 
729 //======================================================================
730  /// Compute norm of fe solution
731 //======================================================================
733  {
734 
735  // Initialise
736  norm=0.0;
737 
738  //Vector of local coordinates
739  Vector<double> s(2);
740 
741  // Vector for coordintes
742  Vector<double> x(2);
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  //Loop over the integration points
753  for(unsigned ipt=0;ipt<n_intpt;ipt++)
754  {
755 
756  //Assign values of s
757  for(unsigned i=0;i<2;i++)
758  {
759  s[i] = integral_pt()->knot(ipt,i);
760  }
761 
762  //Get the integral weight
763  double w = integral_pt()->weight(ipt);
764 
765  // Get jacobian of mapping
766  double J=J_eulerian(s);
767 
768  //Premultiply the weights and the Jacobian
769  double W = w*J;
770 
771  // Get FE function value
772  std::complex<double> u_fe=interpolated_u_fourier_decomposed_helmholtz(s);
773 
774  // Add to norm
775  norm+=(u_fe.real()*u_fe.real()+u_fe.imag()*u_fe.imag())*W;
776 
777  }
778  }
779 
780 
781 
782 
783 
784 //====================================================================
785 // Force build of templates
786 //====================================================================
790 
791 }
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...
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
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
void output(std::ostream &outfile)
Output with default number of plot points.
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
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
std::complex< double > interpolated_u_fourier_decomposed_helmholtz(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Broken assignment operator.
double plgndr1(const unsigned &n, const double &x)
Legendre polynomials depending on one parameter.
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
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
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
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
virtual void get_source_fourier_decomposed_helmholtz(const unsigned &ipt, const Vector< double > &x, std::complex< double > &source) const
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
void compute_norm(double &norm)
Compute norm of fe solution.
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_real_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for real part of full time-dependent fct u = Re( (u_r +i u_i) exp(-i omega t) at phas...
virtual double dshape_and_dtest_eulerian_at_knot_fourier_decomposed_helmholtz(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
double plgndr2(const unsigned &l, const unsigned &m, const double &x)
Legendre polynomials depending on two parameters.
virtual void fill_in_generic_residual_contribution_fourier_decomposed_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.