fourier_decomposed_helmholtz_flux_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 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 elements that are used to apply prescribed flux
31 // boundary conditions to the Fourier decomposed Helmholtz equations
32 #ifndef OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_FLUX_ELEMENTS_HEADER
33 #define OOMPH_FOURIER_DECOMPOSED_HELMHOLTZ_FLUX_ELEMENTS_HEADER
34 
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38 #include <oomph-lib-config.h>
39 #endif
40 
41 #include "math.h"
42 #include <complex>
43 
44 // oomph-lib includes
45 #include "../generic/Qelements.h"
46 
47 
48 namespace oomph
49 {
50 
51 //======================================================================
52 /// \short A class for elements that allow the imposition of an
53 /// applied flux on the boundaries of Fourier decomposed Helmholtz elements.
54 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
55 /// policy class.
56 //======================================================================
57  template <class ELEMENT>
59  public virtual FaceGeometry<ELEMENT>,
60  public virtual FaceElement
61  {
62 
63  public:
64 
65  /// \short Function pointer to the prescribed-flux function fct(x,f(x)) --
66  /// x is a Vector and the flux is a complex
67 
69  (const Vector<double>& x, std::complex<double>& flux);
70 
71  /// \short Constructor, takes the pointer to the "bulk" element and the
72  /// index of the face to which the element is attached.
74  const int& face_index);
75 
76  ///\short Broken empty constructor
78  {
79  throw OomphLibError(
80  "Don't call empty constructor for FourierDecomposedHelmholtzFluxElement",
81  OOMPH_CURRENT_FUNCTION,
82  OOMPH_EXCEPTION_LOCATION);
83  }
84 
85  /// Broken copy constructor
88  {
89  BrokenCopy::broken_copy("FourierDecomposedHelmholtzFluxElement");
90  }
91 
92  /// Broken assignment operator
93 //Commented out broken assignment operator because this can lead to a conflict warning
94 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
95 //realise that two separate implementations of the broken function are the same and so,
96 //quite rightly, it shouts.
97  /*void operator=(const FourierDecomposedHelmholtzFluxElement&)
98  {
99  BrokenCopy::broken_assign("FourierDecomposedHelmholtzFluxElement");
100  }*/
101 
102 
103  /// Access function for the prescribed-flux function pointer
105  {return Flux_fct_pt;}
106 
107 
108  /// Add the element's contribution to its residual vector
110  {
111  //Call the generic residuals function with flag set to 0
112  //using a dummy matrix argument
115  }
116 
117 
118  /// \short Add the element's contribution to its residual vector and its
119  /// Jacobian matrix
121  DenseMatrix<double> &jacobian)
122  {
123  //Call the generic routine with the flag set to 1
125  (residuals,jacobian,1);
126  }
127 
128  /// Output function -- forward to broken version in FiniteElement
129  /// until somebody decides what exactly they want to plot here...
130  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
131 
132  /// \short Output function -- forward to broken version in FiniteElement
133  /// until somebody decides what exactly they want to plot here...
134  void output(std::ostream &outfile, const unsigned &n_plot)
135  {FiniteElement::output(outfile,n_plot);}
136 
137 
138  /// C-style output function -- forward to broken version in FiniteElement
139  /// until somebody decides what exactly they want to plot here...
140  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
141 
142  /// \short C-style output function -- forward to broken version in
143  /// FiniteElement until somebody decides what exactly they want to plot
144  /// here...
145  void output(FILE* file_pt, const unsigned &n_plot)
146  {FiniteElement::output(file_pt,n_plot);}
147 
148 
149  /// \short Return the index at which the unknown value
150  /// is stored. (Real/imag part gives real index of real/imag part).
151  virtual inline std::complex<unsigned>
153  {
154  return std::complex<unsigned>(
157  }
158 
159 
160  protected:
161 
162  /// \short Function to compute the shape and test functions and to return
163  /// the Jacobian of mapping between local and global (Eulerian)
164  /// coordinates
165  inline double shape_and_test(const Vector<double> &s, Shape &psi,
166  Shape &test)
167  const
168  {
169  //Find number of nodes
170  unsigned n_node = nnode();
171 
172  //Get the shape functions
173  shape(s,psi);
174 
175  //Set the test functions to be the same as the shape functions
176  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
177 
178  //Return the value of the jacobian
179  return J_eulerian(s);
180  }
181 
182 
183  /// \short Function to compute the shape and test functions and to return
184  /// the Jacobian of mapping between local and global (Eulerian)
185  /// coordinates
186  inline double shape_and_test_at_knot(const unsigned &ipt,
187  Shape &psi, Shape &test)
188  const
189  {
190  //Find number of nodes
191  unsigned n_node = nnode();
192 
193  //Get the shape functions
194  shape_at_knot(ipt,psi);
195 
196  //Set the test functions to be the same as the shape functions
197  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
198 
199  //Return the value of the jacobian
200  return J_eulerian_at_knot(ipt);
201  }
202 
203 
204  /// Function to calculate the prescribed flux at a given spatial
205  /// position
206  void get_flux(const Vector<double>& x, std::complex<double>& flux)
207  {
208  //If the function pointer is zero return zero
209  if(Flux_fct_pt == 0)
210  {
211  flux = std::complex<double>(0.0,0.0);
212  }
213  //Otherwise call the function
214  else
215  {
216  (*Flux_fct_pt)(x,flux);
217  }
218  }
219 
220 
221  /// \short The index at which the real and imag part of the
222  /// unknown is stored at the nodes
223  std::complex<unsigned> U_index_fourier_decomposed_helmholtz;
224 
225 
226  /// \short Add the element's contribution to its residual vector.
227  /// flag=1(or 0): do (or don't) compute the contribution to the
228  /// Jacobian as well.
229  virtual void
231  Vector<double> &residuals, DenseMatrix<double> &jacobian,
232  const unsigned& flag);
233 
234 
235  /// Function pointer to the (global) prescribed-flux function
237 
238  };
239 
240 //////////////////////////////////////////////////////////////////////
241 //////////////////////////////////////////////////////////////////////
242 //////////////////////////////////////////////////////////////////////
243 
244 
245 
246 //===========================================================================
247 /// Constructor, takes the pointer to the "bulk" element, the
248 /// index of the fixed local coordinate and its value represented
249 /// by an integer (+/- 1), indicating that the face is located
250 /// at the max. or min. value of the "fixed" local coordinate
251 /// in the bulk element.
252 //===========================================================================
253  template<class ELEMENT>
256  const int &face_index) :
257  FaceGeometry<ELEMENT>(), FaceElement()
258  {
259  // Let the bulk element build the FaceElement, i.e. setup the pointers
260  // to its nodes (by referring to the appropriate nodes in the bulk
261  // element), etc.
262  bulk_el_pt->build_face_element(face_index,this);
263 
264  // Initialise the prescribed-flux function pointer to zero
265  Flux_fct_pt = 0;
266 
267 
268  // Initialise index at which real and imag unknowns are stored
269  U_index_fourier_decomposed_helmholtz = std::complex<unsigned>(0,1);
270 
271  // Now read out indices from bulk element
273  dynamic_cast<FourierDecomposedHelmholtzEquations*>(bulk_el_pt);
274  //If the cast has failed die
275  if(eqn_pt==0)
276  {
277  std::string error_string =
278  "Bulk element must inherit from FourierDecomposedHelmholtzEquations.";
279  throw OomphLibError(
280  error_string,
281  OOMPH_CURRENT_FUNCTION,
282  OOMPH_EXCEPTION_LOCATION);
283  }
284  else
285  {
286  //Read the index from the (cast) bulk element
288  eqn_pt->u_index_fourier_decomposed_helmholtz();
289  }
290 
291  }
292 
293 
294 //===========================================================================
295 /// Compute the element's residual vector and the (zero) Jacobian matrix.
296 //===========================================================================
297  template<class ELEMENT>
300  Vector<double> &residuals, DenseMatrix<double> &jacobian,
301  const unsigned& flag)
302  {
303  //Find out how many nodes there are
304  const unsigned n_node = nnode();
305 
306  //Set up memory for the shape and test functions
307  Shape psif(n_node), testf(n_node);
308 
309  //Set the value of Nintpt
310  const unsigned n_intpt = integral_pt()->nweight();
311 
312  //Set the Vector to hold local coordinates
313  Vector<double> s(1);
314 
315  //Integers to hold the local equation and unknown numbers
316  int local_eqn_real=0 ,local_eqn_imag=0;
317 
318  //Loop over the integration points
319  //--------------------------------
320  for(unsigned ipt=0;ipt<n_intpt;ipt++)
321  {
322 
323  //Assign values of s
324  for(unsigned i=0;i<1;i++) {s[i] = integral_pt()->knot(ipt,i);}
325 
326  //Get the integral weight
327  double w = integral_pt()->weight(ipt);
328 
329  //Find the shape and test functions and return the Jacobian
330  //of the mapping
331  double J = shape_and_test(s,psif,testf);
332 
333  //Premultiply the weights and the Jacobian
334  double W = w*J;
335 
336  //Need to find position to feed into flux function, initialise to zero
337  Vector<double> interpolated_x(2,0.0);
338 
339  //Calculate coordinates
340  for(unsigned l=0;l<n_node;l++)
341  {
342  //Loop over coordinates
343  for(unsigned i=0;i<2;i++)
344  {
345  interpolated_x[i] += nodal_position(l,i)*psif[l];
346  }
347  }
348 
349  //first component
350  double r = interpolated_x[0];
351 
352  //Get the imposed flux
353  std::complex<double> flux(0.0,0.0);
354  get_flux(interpolated_x,flux);
355 
356  //Now add to the appropriate equations
357  //Loop over the test functions
358  for(unsigned l=0;l<n_node;l++)
359  {
360  local_eqn_real =
361  nodal_local_eqn(l,U_index_fourier_decomposed_helmholtz.real());
362 
363  /*IF it's not a boundary condition*/
364  if(local_eqn_real >= 0)
365  {
366  //Add the prescribed flux terms
367  residuals[local_eqn_real] -= flux.real()*testf[l]*r*W;
368 
369  // Imposed traction doesn't depend upon the solution,
370  // --> the Jacobian is always zero, so no Jacobian
371  // terms are required
372  }
373  local_eqn_imag =
374  nodal_local_eqn(l,U_index_fourier_decomposed_helmholtz.imag());
375 
376  /*IF it's not a boundary condition*/
377  if(local_eqn_imag >= 0)
378  {
379  //Add the prescribed flux terms
380  residuals[local_eqn_imag] -= flux.imag()*testf[l]*r*W;
381 
382  // Imposed traction doesn't depend upon the solution,
383  // --> the Jacobian is always zero, so no Jacobian
384  // terms are required
385  }
386  }
387  }
388  }
389 
390 }
391 
392 #endif
393 
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
FourierDecomposedHelmholtzFluxElement(const FourierDecomposedHelmholtzFluxElement &dummy)
Broken copy constructor.
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
Definition: elements.cc:5121
FourierDecomposedHelmholtzPrescribedFluxFctPt & flux_fct_pt()
Broken assignment operator.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
A general Finite Element class.
Definition: elements.h:1271
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the element's contribution to its residual vector and its Jacobian matrix.
double shape_and_test(const Vector< double > &s, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping between local ...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
static char t char * s
Definition: cfortran.h:572
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3158
std::complex< unsigned > U_index_fourier_decomposed_helmholtz
The index at which the real and imag part of the unknown is stored at the nodes.
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:5033
FourierDecomposedHelmholtzPrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
A class for elements that allow the imposition of an applied flux on the boundaries of Fourier decomp...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* FourierDecomposedHelmholtzPrescribedFluxFctPt)(const Vector< double > &x, std::complex< double > &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector and the flux is a comple...
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
void get_flux(const Vector< double > &x, std::complex< double > &flux)
virtual std::complex< unsigned > u_index_fourier_decomposed_helmholtz() const
Return the index at which the unknown value is stored. (Real/imag part gives real index of real/imag ...
virtual void fill_in_generic_residual_contribution_fourier_decomposed_helmholtz_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element's contribution to its residual vector. flag=1(or 0): do (or don't) compute the contri...
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...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
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