helmholtz_elements.h
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: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
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 //Header file for Helmholtz elements
31 #ifndef OOMPH_HELMHOLTZ_ELEMENTS_HEADER
32 #define OOMPH_HELMHOLTZ_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37  #include <oomph-lib-config.h>
38 #endif
39 
40 
41 //OOMPH-LIB headers
42 #include "../generic/projection.h"
43 #include "../generic/nodes.h"
44 #include "../generic/Qelements.h"
45 #include "../generic/oomph_utilities.h"
46 #include "math.h"
47 #include <complex>
48 
49 namespace oomph
50 {
51 
52 //=============================================================
53 /// A class for all isoparametric elements that solve the
54 /// Helmholtz equations.
55 /// \f[
56 /// \frac{\partial^2 u}{\partial x_i^2} + k^2 u = f(x_j)
57 /// \f]
58 /// This contains the generic maths. Shape functions, geometric
59 /// mapping etc. must get implemented in derived class.
60 //=============================================================
61 template <unsigned DIM>
62 class HelmholtzEquations : public virtual FiniteElement
63 {
64 
65 public:
66 
67  /// \short Function pointer to source function fct(x,f(x)) --
68  /// x is a Vector!
69  typedef void (*HelmholtzSourceFctPt)(const Vector<double>& x,
70  std::complex<double>& f);
71 
72 
73  /// Constructor (must initialise the Source_fct_pt to null)
75  K_squared_pt(0)
76  {}
77 
78  /// Broken copy constructor
80  {
81  BrokenCopy::broken_copy("HelmholtzEquations");
82  }
83 
84  /// Broken assignment operator
85 //Commented out broken assignment operator because this can lead to a conflict warning
86 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
87 //realise that two separate implementations of the broken function are the same and so,
88 //quite rightly, it shouts.
89  /*void operator=(const HelmholtzEquations&)
90  {
91  BrokenCopy::broken_assign("HelmholtzEquations");
92  }*/
93 
94  /// \short Return the index at which the unknown value
95  /// is stored.
96  virtual inline std::complex<unsigned> u_index_helmholtz() const
97  {return std::complex<unsigned>(0,1);}
98 
99 
100  /// Get pointer to square of wavenumber
101  double*& k_squared_pt()
102  {
103  return K_squared_pt;
104  }
105 
106 
107  /// Get the square of wavenumber
108  double k_squared()
109  {
110 #ifdef PARANOID
111  if (K_squared_pt==0)
112  {
113  throw OomphLibError(
114  "Please set pointer to k_squared using access fct to pointer!",
115  OOMPH_CURRENT_FUNCTION,
116  OOMPH_EXCEPTION_LOCATION);
117  }
118 #endif
119  return *K_squared_pt;
120 
121  }
122 
123 
124  /// \short Number of scalars/fields output by this element. Reimplements
125  /// broken virtual function in base class.
126  unsigned nscalar_paraview() const
127  {
128  return 2;
129  }
130 
131  /// \short Write values of the i-th scalar field at the plot points. Needs
132  /// to be implemented for each new specific element type.
133  void scalar_value_paraview(std::ofstream& file_out,
134  const unsigned& i,
135  const unsigned& nplot) const
136  {
137  //Vector of local coordinates
138  Vector<double> s(DIM);
139 
140  // Loop over plot points
141  unsigned num_plot_points=nplot_points_paraview(nplot);
142  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
143  {
144  // Get local coordinates of plot point
145  get_s_plot(iplot,nplot,s);
146  std::complex<double> u(interpolated_u_helmholtz(s));
147 
148  // Paraview need to ouput the fileds seperatly so it loops through all
149  // the elements twice
150  switch(i)
151  {
152  // Real part first
153  case 0:
154  file_out << u.real() << std::endl;
155  break;
156 
157  // Imaginary part second
158  case 1:
159  file_out << u.imag() << std::endl;
160  break;
161 
162  // Never get here
163  default:
164  std::stringstream error_stream;
165  error_stream
166  << "Helmholtz elements only store 2 fields so i must be 0 or 1"
167  << std::endl;
168  throw OomphLibError(
169  error_stream.str(),
170  OOMPH_CURRENT_FUNCTION,
171  OOMPH_EXCEPTION_LOCATION);
172  break;
173  }
174  } // end of plotpoint loop
175  } // end scalar_value_paraview
176 
177  /// \short Name of the i-th scalar field. Default implementation
178  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
179  /// overloaded with more meaningful names in specific elements.
180  std::string scalar_name_paraview(const unsigned& i) const
181  {
182  switch(i)
183  {
184  case 0:
185  return "Real part";
186  break;
187 
188  case 1:
189  return "Imaginary part";
190  break;
191 
192  // Never get here
193  default:
194  std::stringstream error_stream;
195  error_stream
196  << "Helmholtz elements only store 2 fields so i must be 0 or 1"
197  << std::endl;
198  throw OomphLibError(
199  error_stream.str(),
200  OOMPH_CURRENT_FUNCTION,
201  OOMPH_EXCEPTION_LOCATION);
202 
203  // Dummy return for the default statement
204  return " ";
205  break;
206  }
207  }
208 
209  /// Output with default number of plot points
210  void output(std::ostream &outfile)
211  {
212  const unsigned n_plot=5;
213  output(outfile,n_plot);
214  }
215 
216  /// \short Output FE representation of soln: x,y,u_re,u_im or
217  /// x,y,z,u_re,u_im at n_plot^DIM plot points
218  void output(std::ostream &outfile, const unsigned &n_plot);
219 
220  /// \short Output function for real part of full time-dependent solution
221  /// u = Re( (u_r +i u_i) exp(-i omega t)
222  /// at phase angle omega t = phi.
223  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
224  /// direction
225  void output_real(std::ostream &outfile, const double& phi,
226  const unsigned &n_plot);
227 
228  /// C_style output with default number of plot points
229  void output(FILE* file_pt)
230  {
231  const unsigned n_plot=5;
232  output(file_pt,n_plot);
233  }
234 
235  /// \short C-style output FE representation of soln: x,y,u_re,u_im or
236  /// x,y,z,u_re,u_im at n_plot^DIM plot points
237  void output(FILE* file_pt, const unsigned &n_plot);
238 
239  /// Output exact soln: x,y,u_re_exact,u_im_exact
240  /// or x,y,z,u_re_exact,u_im_exact at n_plot^DIM plot points
241  void output_fct(std::ostream &outfile, const unsigned &n_plot,
243 
244  /// \short Output exact soln: (dummy time-dependent version to
245  /// keep intel compiler happy)
246  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
247  const double& time,
249  exact_soln_pt)
250  {
251  throw OomphLibError(
252  "There is no time-dependent output_fct() for Helmholtz elements ",
253  OOMPH_CURRENT_FUNCTION,
254  OOMPH_EXCEPTION_LOCATION);
255  }
256 
257 
258 
259  /// \short Output function for real part of full time-dependent fct
260  /// u = Re( (u_r +i u_i) exp(-i omega t)
261  /// at phase angle omega t = phi.
262  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
263  /// direction
264  void output_real_fct(std::ostream &outfile,
265  const double& phi,
266  const unsigned &n_plot,
268 
269 
270  /// Get error against and norm of exact solution
271  void compute_error(std::ostream &outfile,
273  double& error, double& norm);
274 
275 
276  /// Dummy, time dependent error checker
277  void compute_error(std::ostream &outfile,
279  const double& time, double& error, double& norm)
280  {
281  throw OomphLibError(
282  "There is no time-dependent compute_error() for Helmholtz elements",
283  OOMPH_CURRENT_FUNCTION,
284  OOMPH_EXCEPTION_LOCATION);
285  }
286 
287  /// Access function: Pointer to source function
289 
290  /// Access function: Pointer to source function. Const version
292 
293 
294  /// Get source term at (Eulerian) position x. This function is
295  /// virtual to allow overloading in multi-physics problems where
296  /// the strength of the source function might be determined by
297  /// another system of equations.
298  inline virtual void get_source_helmholtz(const unsigned& ipt,
299  const Vector<double>& x,
300  std::complex<double>& source) const
301  {
302  //If no source function has been set, return zero
303  if(Source_fct_pt==0)
304  {
305  source = std::complex<double>(0.0,0.0);
306  }
307  else
308  {
309  // Get source strength
310  (*Source_fct_pt)(x,source);
311  }
312  }
313 
314 
315  /// Get flux: flux[i] = du/dx_i for real and imag part
316  void get_flux(const Vector<double>& s,
317  Vector<std::complex <double> >& flux) const
318  {
319  //Find out how many nodes there are in the element
320  const unsigned n_node = nnode();
321 
322 
323  //Set up memory for the shape and test functions
324  Shape psi(n_node);
325  DShape dpsidx(n_node,DIM);
326 
327  //Call the derivatives of the shape and test functions
328  dshape_eulerian(s,psi,dpsidx);
329 
330  //Initialise to zero
331  const std::complex<double> zero(0.0,0.0);
332  for(unsigned j=0;j<DIM;j++)
333  {
334  flux[j] = zero;
335  }
336 
337  // Loop over nodes
338  for(unsigned l=0;l<n_node;l++)
339  {
340  //Cache the complex value of the unknown
341  const std::complex<double> u_value(
342  this->nodal_value(l,u_index_helmholtz().real()),
343  this->nodal_value(l,u_index_helmholtz().imag()));
344  //Loop over derivative directions
345  for(unsigned j=0;j<DIM;j++)
346  {
347  flux[j] += u_value*dpsidx(l,j);
348  }
349  }
350  }
351 
352 
353  /// Add the element's contribution to its residual vector (wrapper)
355  {
356  //Call the generic residuals function with flag set to 0
357  //using a dummy matrix argument
360  }
361 
362 
363 
364  /// \short Add the element's contribution to its residual vector and
365  /// element Jacobian matrix (wrapper)
367  DenseMatrix<double> &jacobian)
368  {
369  //Call the generic routine with the flag set to 1
371  }
372 
373 
374 
375  /// \short Return FE representation of function value u_helmholtz(s)
376  /// at local coordinate s
377  inline std::complex<double> interpolated_u_helmholtz(const Vector<double> &s) const
378  {
379  //Find number of nodes
380  const unsigned n_node = nnode();
381 
382  //Local shape function
383  Shape psi(n_node);
384 
385  //Find values of shape function
386  shape(s,psi);
387 
388  //Initialise value of u
389  std::complex<double> interpolated_u(0.0,0.0);
390 
391  //Get the index at which the helmholtz unknown is stored
392  const unsigned u_nodal_index_real = u_index_helmholtz().real();
393  const unsigned u_nodal_index_imag = u_index_helmholtz().imag();
394 
395  //Loop over the local nodes and sum
396  for(unsigned l=0;l<n_node;l++)
397  {
398  //Make a temporary complex number from the stored data
399  const std::complex<double> u_value(
400  this->nodal_value(l,u_nodal_index_real),
401  this->nodal_value(l,u_nodal_index_imag));
402  //Add to the interpolated value
403  interpolated_u += u_value*psi[l];
404  }
405  return interpolated_u;
406  }
407 
408 
409  /// \short Self-test: Return 0 for OK
410  unsigned self_test();
411 
412 
413 protected:
414 
415  /// \short Shape/test functions and derivs w.r.t. to global coords at
416  /// local coord. s; return Jacobian of mapping
418  Shape &psi,
419  DShape &dpsidx, Shape &test,
420  DShape &dtestdx) const=0;
421 
422 
423  /// \short Shape/test functions and derivs w.r.t. to global coords at
424  /// integration point ipt; return Jacobian of mapping
425  virtual double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned &ipt,
426  Shape &psi,
427  DShape &dpsidx,
428  Shape &test,
429  DShape &dtestdx)
430  const=0;
431 
432  /// \short Compute element residual Vector only (if flag=and/or element
433  /// Jacobian matrix
435  Vector<double> &residuals, DenseMatrix<double> &jacobian,
436  const unsigned& flag);
437 
438  /// Pointer to source function:
440 
441  /// Pointer to square of wavenumber
442  double* K_squared_pt;
443 
444 };
445 
446 
447 
448 
449 
450 
451 ///////////////////////////////////////////////////////////////////////////
452 ///////////////////////////////////////////////////////////////////////////
453 ///////////////////////////////////////////////////////////////////////////
454 
455 
456 
457 //======================================================================
458 /// QHelmholtzElement elements are linear/quadrilateral/brick-shaped
459 /// Helmholtz elements with isoparametric interpolation for the function.
460 //======================================================================
461 template <unsigned DIM, unsigned NNODE_1D>
462  class QHelmholtzElement : public virtual QElement<DIM,NNODE_1D>,
463  public virtual HelmholtzEquations<DIM>
464 {
465 
466 private:
467 
468  /// \short Static int that holds the number of variables at
469  /// nodes: always the same
470  static const unsigned Initial_Nvalue;
471 
472  public:
473 
474 
475  ///\short Constructor: Call constructors for QElement and
476  /// Helmholtz equations
477  QHelmholtzElement() : QElement<DIM,NNODE_1D>(), HelmholtzEquations<DIM>()
478  {}
479 
480  /// Broken copy constructor
482  {
483  BrokenCopy::broken_copy("QHelmholtzElement");
484  }
485 
486  /// Broken assignment operator
487  /*void operator=(const QHelmholtzElement<DIM,NNODE_1D>&)
488  {
489  BrokenCopy::broken_assign("QHelmholtzElement");
490  }*/
491 
492 
493  /// \short Required # of `values' (pinned or dofs)
494  /// at node n
495  inline unsigned required_nvalue(const unsigned &n) const
496  {return Initial_Nvalue;}
497 
498  /// \short Output function:
499  /// x,y,u or x,y,z,u
500  void output(std::ostream &outfile)
502 
503 
504  /// \short Output function:
505  /// x,y,u or x,y,z,u at n_plot^DIM plot points
506  void output(std::ostream &outfile, const unsigned &n_plot)
507  {HelmholtzEquations<DIM>::output(outfile,n_plot);}
508 
509  /// \short Output function for real part of full time-dependent solution
510  /// u = Re( (u_r +i u_i) exp(-i omega t)
511  /// at phase angle omega t = phi.
512  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
513  /// direction
514  void output_real(std::ostream &outfile, const double& phi,
515  const unsigned &n_plot)
516  {HelmholtzEquations<DIM>::output_real(outfile,phi,n_plot);}
517 
518 
519  /// \short C-style output function:
520  /// x,y,u or x,y,z,u
521  void output(FILE* file_pt)
523 
524 
525  /// \short C-style output function:
526  /// x,y,u or x,y,z,u at n_plot^DIM plot points
527  void output(FILE* file_pt, const unsigned &n_plot)
528  {HelmholtzEquations<DIM>::output(file_pt,n_plot);}
529 
530 
531  /// \short Output function for an exact solution:
532  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
533  void output_fct(std::ostream &outfile, const unsigned &n_plot,
535  {HelmholtzEquations<DIM>::output_fct(outfile,n_plot,exact_soln_pt);}
536 
537 
538  /// \short Output function for real part of full time-dependent fct
539  /// u = Re( (u_r +i u_i) exp(-i omega t)
540  /// at phase angle omega t = phi.
541  /// x,y,u or x,y,z,u at n_plot plot points in each coordinate
542  /// direction
543  void output_real_fct(std::ostream &outfile,
544  const double& phi,
545  const unsigned &n_plot,
547  {
548  HelmholtzEquations<DIM>::output_real_fct(outfile,phi,n_plot,exact_soln_pt);
549  }
550 
551 
552  /// \short Output function for a time-dependent exact solution.
553  /// x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points
554  /// (Calls the steady version)
555  void output_fct(std::ostream &outfile, const unsigned &n_plot,
556  const double& time,
558  {HelmholtzEquations<DIM>::output_fct(outfile,n_plot,time,exact_soln_pt);}
559 
560 
561 protected:
562 
563 /// Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
565  const Vector<double> &s, Shape &psi, DShape &dpsidx,
566  Shape &test, DShape &dtestdx) const;
567 
568 
569  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
570  /// integration point ipt. Return Jacobian.
571  inline double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned& ipt,
572  Shape &psi,
573  DShape &dpsidx,
574  Shape &test,
575  DShape &dtestdx)
576  const;
577 
578 };
579 
580 
581 
582 
583 //Inline functions:
584 
585 
586 //======================================================================
587 /// Define the shape functions and test functions and derivatives
588 /// w.r.t. global coordinates and return Jacobian of mapping.
589 ///
590 /// Galerkin: Test functions = shape functions
591 //======================================================================
592 template<unsigned DIM, unsigned NNODE_1D>
594  const Vector<double> &s,
595  Shape &psi,
596  DShape &dpsidx,
597  Shape &test,
598  DShape &dtestdx) const
599 {
600  //Call the geometrical shape functions and derivatives
601  const double J = this->dshape_eulerian(s,psi,dpsidx);
602 
603  //Set the test functions equal to the shape functions
604  test = psi;
605  dtestdx= dpsidx;
606 
607  //Return the jacobian
608  return J;
609 }
610 
611 
612 
613 
614 //======================================================================
615 /// Define the shape functions and test functions and derivatives
616 /// w.r.t. global coordinates and return Jacobian of mapping.
617 ///
618 /// Galerkin: Test functions = shape functions
619 //======================================================================
620 template<unsigned DIM, unsigned NNODE_1D>
623  const unsigned &ipt,
624  Shape &psi,
625  DShape &dpsidx,
626  Shape &test,
627  DShape &dtestdx) const
628 {
629  //Call the geometrical shape functions and derivatives
630  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
631 
632  //Set the pointers of the test functions
633  test = psi;
634  dtestdx = dpsidx;
635 
636  //Return the jacobian
637  return J;
638 }
639 
640 ////////////////////////////////////////////////////////////////////////
641 ////////////////////////////////////////////////////////////////////////
642 ////////////////////////////////////////////////////////////////////////
643 
644 
645 
646 //=======================================================================
647 /// Face geometry for the QHelmholtzElement elements: The spatial
648 /// dimension of the face elements is one lower than that of the
649 /// bulk element but they have the same number of points
650 /// along their 1D edges.
651 //=======================================================================
652 template<unsigned DIM, unsigned NNODE_1D>
653 class FaceGeometry<QHelmholtzElement<DIM,NNODE_1D> >:
654  public virtual QElement<DIM-1,NNODE_1D>
655 {
656 
657  public:
658 
659  /// \short Constructor: Call the constructor for the
660  /// appropriate lower-dimensional QElement
661  FaceGeometry() : QElement<DIM-1,NNODE_1D>() {}
662 
663 };
664 
665 ////////////////////////////////////////////////////////////////////////
666 ////////////////////////////////////////////////////////////////////////
667 ////////////////////////////////////////////////////////////////////////
668 
669 
670 //=======================================================================
671 /// Face geometry for the 1D QHelmholtzElement elements: Point elements
672 //=======================================================================
673 template<unsigned NNODE_1D>
674 class FaceGeometry<QHelmholtzElement<1,NNODE_1D> >:
675  public virtual PointElement
676 {
677 
678  public:
679 
680  /// \short Constructor: Call the constructor for the
681  /// appropriate lower-dimensional QElement
683 
684 };
685 
686 
687 
688 ////////////////////////////////////////////////////////////////////////
689 ////////////////////////////////////////////////////////////////////////
690 ////////////////////////////////////////////////////////////////////////
691 
692 
693 
694 //==========================================================
695 /// Helmholtz upgraded to become projectable
696 //==========================================================
697  template<class HELMHOLTZ_ELEMENT>
699  public virtual ProjectableElement<HELMHOLTZ_ELEMENT>
700  {
701 
702  public:
703 
704 
705  /// \short Constructor [this was only required explicitly
706  /// from gcc 4.5.2 onwards...]
708 
709  /// \short Specify the values associated with field fld.
710  /// The information is returned in a vector of pairs which comprise
711  /// the Data object and the value within it, that correspond to field fld.
713  {
714 
715 #ifdef PARANOID
716  if (fld>1)
717  {
718  std::stringstream error_stream;
719  error_stream
720  << "Helmholtz elements only store 2 fields so fld = "
721  << fld << " is illegal \n";
722  throw OomphLibError(
723  error_stream.str(),
724  OOMPH_CURRENT_FUNCTION,
725  OOMPH_EXCEPTION_LOCATION);
726  }
727 #endif
728 
729  // Create the vector
730  unsigned nnod=this->nnode();
731  Vector<std::pair<Data*,unsigned> > data_values(nnod);
732 
733  // Loop over all nodes
734  for (unsigned j=0;j<nnod;j++)
735  {
736  // Add the data value associated field: The node itself
737  data_values[j]=std::make_pair(this->node_pt(j),fld);
738  }
739 
740  // Return the vector
741  return data_values;
742  }
743 
744  /// \short Number of fields to be projected: 2 (real and imag part)
746  {
747  return 2;
748  }
749 
750  /// \short Number of history values to be stored for fld-th field.
751  /// (includes current value!)
752  unsigned nhistory_values_for_projection(const unsigned &fld)
753  {
754 #ifdef PARANOID
755  if (fld>1)
756  {
757  std::stringstream error_stream;
758  error_stream
759  << "Helmholtz elements only store two fields so fld = "
760  << fld << " is illegal\n";
761  throw OomphLibError(
762  error_stream.str(),
763  OOMPH_CURRENT_FUNCTION,
764  OOMPH_EXCEPTION_LOCATION);
765  }
766 #endif
767  return this->node_pt(0)->ntstorage();
768  }
769 
770  ///\short Number of positional history values
771  /// (includes current value!)
773  {
774  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
775  }
776 
777  /// \short Return Jacobian of mapping and shape functions of field fld
778  /// at local coordinate s
779  double jacobian_and_shape_of_field(const unsigned &fld,
780  const Vector<double> &s,
781  Shape &psi)
782  {
783 #ifdef PARANOID
784  if (fld>1)
785  {
786  std::stringstream error_stream;
787  error_stream
788  << "Helmholtz elements only store two fields so fld = "
789  << fld << " is illegal.\n";
790  throw OomphLibError(
791  error_stream.str(),
792  OOMPH_CURRENT_FUNCTION,
793  OOMPH_EXCEPTION_LOCATION);
794  }
795 #endif
796  unsigned n_dim=this->dim();
797  unsigned n_node=this->nnode();
798  Shape test(n_node);
799  DShape dpsidx(n_node,n_dim), dtestdx(n_node,n_dim);
800  double J=this->dshape_and_dtest_eulerian_helmholtz(s,psi,dpsidx,
801  test,dtestdx);
802  return J;
803  }
804 
805 
806 
807  /// \short Return interpolated field fld at local coordinate s, at time level
808  /// t (t=0: present; t>0: history values)
809  double get_field(const unsigned &t,
810  const unsigned &fld,
811  const Vector<double>& s)
812  {
813 #ifdef PARANOID
814  if (fld>1)
815  {
816  std::stringstream error_stream;
817  error_stream
818  << "Helmholtz elements only store two fields so fld = "
819  << fld << " is illegal\n";
820  throw OomphLibError(
821  error_stream.str(),
822  OOMPH_CURRENT_FUNCTION,
823  OOMPH_EXCEPTION_LOCATION);
824  }
825 #endif
826  //Find the index at which the variable is stored
827  std::complex<unsigned> complex_u_nodal_index = this->u_index_helmholtz();
828  unsigned u_nodal_index = 0;
829  if (fld==0)
830  {
831  u_nodal_index = complex_u_nodal_index.real();
832  }
833  else
834  {
835  u_nodal_index = complex_u_nodal_index.imag();
836  }
837 
838 
839  //Local shape function
840  unsigned n_node=this->nnode();
841  Shape psi(n_node);
842 
843  //Find values of shape function
844  this->shape(s,psi);
845 
846  //Initialise value of u
847  double interpolated_u = 0.0;
848 
849  //Sum over the local nodes
850  for(unsigned l=0;l<n_node;l++)
851  {
852  interpolated_u += this->nodal_value(t,l,u_nodal_index)*psi[l];
853  }
854  return interpolated_u;
855  }
856 
857 
858 
859 
860  ///Return number of values in field fld: One per node
861  unsigned nvalue_of_field(const unsigned &fld)
862  {
863 #ifdef PARANOID
864  if (fld>1)
865  {
866  std::stringstream error_stream;
867  error_stream
868  << "Helmholtz elements only store two fields so fld = "
869  << fld << " is illegal\n";
870  throw OomphLibError(
871  error_stream.str(),
872  OOMPH_CURRENT_FUNCTION,
873  OOMPH_EXCEPTION_LOCATION);
874  }
875 #endif
876  return this->nnode();
877  }
878 
879 
880  ///Return local equation number of value j in field fld.
881  int local_equation(const unsigned &fld,
882  const unsigned &j)
883  {
884 #ifdef PARANOID
885  if (fld>1)
886  {
887  std::stringstream error_stream;
888  error_stream
889  << "Helmholtz elements only store two fields so fld = "
890  << fld << " is illegal\n";
891  throw OomphLibError(
892  error_stream.str(),
893  OOMPH_CURRENT_FUNCTION,
894  OOMPH_EXCEPTION_LOCATION);
895  }
896 #endif
897  std::complex<unsigned> complex_u_nodal_index = this->u_index_helmholtz();
898  unsigned u_nodal_index = 0;
899  if (fld==0)
900  {
901  u_nodal_index = complex_u_nodal_index.real();
902  }
903  else
904  {
905  u_nodal_index = complex_u_nodal_index.imag();
906  }
907  return this->nodal_local_eqn(j,u_nodal_index);
908  }
909 
910 
911 
912 
913  /// \short Output FE representation of soln: x,y,u or x,y,z,u at
914  /// n_plot^DIM plot points
915  void output(std::ostream &outfile, const unsigned &nplot)
916  {
917  HELMHOLTZ_ELEMENT::output(outfile,nplot);
918  }
919 
920 
921  };
922 
923 
924 //=======================================================================
925 /// Face geometry for element is the same as that for the underlying
926 /// wrapped element
927 //=======================================================================
928  template<class ELEMENT>
930  : public virtual FaceGeometry<ELEMENT>
931  {
932  public:
933  FaceGeometry() : FaceGeometry<ELEMENT>() {}
934  };
935 
936 
937 //=======================================================================
938 /// Face geometry of the Face Geometry for element is the same as
939 /// that for the underlying wrapped element
940 //=======================================================================
941  template<class ELEMENT>
943  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
944  {
945  public:
947  };
948 
949 
950 
951 
952 }
953 
954 #endif
HelmholtzSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector (wrapper)
double k_squared()
Get the square of wavenumber.
double dshape_and_dtest_eulerian_at_knot_helmholtz(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt. Return Jacobian.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
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
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln: x,y,u or x,y,z,u at n_plot^DIM plot points.
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
HelmholtzSourceFctPt Source_fct_pt
Pointer to source function:
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
void output(FILE *file_pt)
C_style output with default number of plot points.
QHelmholtzElement()
Constructor: Call constructors for QElement and Helmholtz equations.
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (includes current value!) ...
virtual std::complex< unsigned > u_index_helmholtz() const
Broken assignment operator.
A general Finite Element class.
Definition: elements.h:1271
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...
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
Helmholtz upgraded to become projectable.
virtual double dshape_and_dtest_eulerian_helmholtz(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
virtual void get_source_helmholtz(const unsigned &ipt, const Vector< double > &x, std::complex< double > &source) const
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
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.
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
unsigned self_test()
Self-test: Return 0 for OK.
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: (dummy time-dependent version to keep intel compiler happy)
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
void output(std::ostream &outfile)
Output function: x,y,u or x,y,z,u.
unsigned nfields_for_projection()
Number of fields to be projected: 2 (real and imag part)
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
void(* HelmholtzSourceFctPt)(const Vector< double > &x, std::complex< double > &f)
Function pointer to source function fct(x,f(x)) – x is a Vector!
double dshape_and_dtest_eulerian_helmholtz(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points...
ProjectableHelmholtzElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (includes current value!)
static char t char * s
Definition: cfortran.h:572
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot ...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
void get_flux(const Vector< double > &s, Vector< std::complex< double > > &flux) const
Get flux: flux[i] = du/dx_i for real and imag part.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double *& k_squared_pt()
Get pointer to square of wavenumber.
void output_real_fct(std::ostream &outfile, const double &phi, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for real part of full time-dependent fct u = Re( (u_r +i u_i) exp(-i omega t) at phas...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and element Jacobian matrix (wrapper) ...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
QHelmholtzElement(const QHelmholtzElement< DIM, NNODE_1D > &dummy)
Broken copy constructor.
double * K_squared_pt
Pointer to square of wavenumber.
void output(std::ostream &outfile)
Output with default number of plot points.
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
HelmholtzSourceFctPt source_fct_pt() const
Access function: Pointer to source function. Const version.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld: One per node.
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 double dshape_and_dtest_eulerian_at_knot_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 ...
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement. ...
HelmholtzEquations(const HelmholtzEquations &dummy)
Broken copy constructor.
void output(std::ostream &outfile)
Output with default number of plot points.
void output(FILE *file_pt)
C-style output function: x,y,u or x,y,z,u.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: x,y,u or x,y,z,u at n_plot^DIM plot points.
HelmholtzEquations()
Constructor (must initialise the Source_fct_pt to null)
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
std::complex< double > interpolated_u_helmholtz(const Vector< double > &s) const
Return FE representation of function value u_helmholtz(s) at local coordinate s.
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.