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