linear_wave_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: 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 // Header file for elements that are used to apply prescribed flux
31 // boundary conditions to the LinearWave equations
32 #ifndef OOMPH_WAVE_FLUX_ELEMENTS_HEADER
33 #define OOMPH_WAVE_FLUX_ELEMENTS_HEADER
34 
35 #ifdef HAVE_CONFIG_H
36 #include <oomph-lib-config.h>
37 #endif
38 
39 // oomph-lib ncludes
40 #include "../generic/Qelements.h"
41 
42 namespace oomph
43 {
44 
45 ////////////////////////////////////////////////////////////////////////
46 ////////////////////////////////////////////////////////////////////////
47 ////////////////////////////////////////////////////////////////////////
48 
49 //======================================================================
50 /// \short A class for elements that allow the imposition of an
51 /// applied flux on the boundaries of LinearWave elements.
52 /// The element geometry is obtained from the FaceGeometry<ELEMENT>
53 /// policy class.
54 //======================================================================
55 template <class ELEMENT>
56 class LinearWaveFluxElement : public virtual FaceGeometry<ELEMENT>,
57 public virtual FaceElement
58 {
59 
60 public:
61 
62  /// \short Function pointer to the prescribed-flux function fct(x,f(x)) --
63  /// x is a Vector!
64  typedef void (*LinearWavePrescribedFluxFctPt)
65  (const double& time, const Vector<double>& x, double& flux);
66 
67 
68  /// \short Constructor, takes the pointer to the "bulk" element and the
69  /// index of the face to be created
70  LinearWaveFluxElement(FiniteElement *bulk_el_pt, const int &face_index);
71 
72 
73  ///\short Broken empty constructor
75  {
76  throw OomphLibError(
77  "Don't call empty constructor for LinearWaveFluxElement",
78  OOMPH_CURRENT_FUNCTION,
79  OOMPH_EXCEPTION_LOCATION);
80  }
81 
82  /// Broken copy constructor
84  {
85  BrokenCopy::broken_copy("LinearWaveFluxElement");
86  }
87 
88  /// Broken assignment operator
90  {
91  BrokenCopy::broken_assign("LinearWaveFluxElement");
92  }
93 
94 
95  /// Access function for the prescribed-flux function pointer
97 
98 
99  /// Compute the element residual vector
101  {
102  //Call the generic residuals function with flag set to 0
103  //using the dummy matrix argument
106  }
107 
108 
109  /// Compute the element's residual vector and its Jacobian matrix
111  DenseMatrix<double> &jacobian)
112  {
113  //Call the generic routine with the flag set to 1
115  }
116 
117  /// Specify the value of nodal zeta from the face geometry
118  /// \short The "global" intrinsic coordinate of the element when
119  /// viewed as part of a geometric object should be given by
120  /// the FaceElement representation, by default (needed to break
121  /// indeterminacy if bulk element is SolidElement)
122  double zeta_nodal(const unsigned &n, const unsigned &k,
123  const unsigned &i) const
124  {return FaceElement::zeta_nodal(n,k,i);}
125 
126  /// Output function -- forward to broken version in FiniteElement
127  /// until somebody decides what exactly they want to plot here...
128  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
129 
130  /// \short Output function -- forward to broken version in FiniteElement
131  /// until somebody decides what exactly they want to plot here...
132  void output(std::ostream &outfile, const unsigned &n_plot)
133  {FiniteElement::output(outfile,n_plot);}
134 
135 
136  /// Output function -- forward to broken version in FiniteElement
137  /// until somebody decides what exactly they want to plot here...
138  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
139 
140  /// \short Output function -- forward to broken version in FiniteElement
141  /// until somebody decides what exactly they want to plot here...
142  void output(FILE* file_pt, const unsigned &n_plot)
143  {FiniteElement::output(file_pt,n_plot);}
144 
145 
146 protected:
147 
148  /// \short Function to compute the shape and test functions and to return
149  /// the Jacobian of mapping between local and global (Eulerian)
150  /// coordinates
151  inline double shape_and_test(const Vector<double> &s, Shape &psi, Shape &test)
152  const
153  {
154  //Find number of nodes
155  unsigned n_node = nnode();
156 
157  //Get the shape functions
158  shape(s,psi);
159 
160  //Set the test functions to be the same as the shape functions
161  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
162 
163  //Return the value of the jacobian
164  return J_eulerian(s);
165  }
166 
167 
168 
169  /// Function to calculate the prescribed flux at a given spatial
170  /// position and at a gien time
171  void get_flux(const double& time, const Vector<double>& x, double& flux)
172  {
173  //If the function pointer is zero return zero
174  if(Flux_fct_pt == 0)
175  {
176  flux=0.0;
177  }
178  //Otherwise call the function
179  else
180  {
181  (*Flux_fct_pt)(time,x,flux);
182  }
183  }
184 
185 
186 
187 private:
188 
189 
190  /// \short Compute the element residual vector.
191  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
193  Vector<double> &residuals,
194  DenseMatrix<double> &jacobian,
195  unsigned flag);
196 
197  /// Function pointer to the (global) prescribed-flux function
199 
200  ///The spatial dimension of the problem
201  unsigned Dim;
202 
203  /// The index at which the unknown is stored at the nodes
205 
206 };
207 
208 
209 ////////////////////////////////////////////////////////////////////////
210 ////////////////////////////////////////////////////////////////////////
211 ////////////////////////////////////////////////////////////////////////
212 
213 
214 
215 //===========================================================================
216 /// Constructor, takes the pointer to the "bulk" element and the
217 /// index of the face to be created.
218 //===========================================================================
219 template<class ELEMENT>
221  const int &face_index) :
222  FaceGeometry<ELEMENT>(), FaceElement()
223  {
224 
225  // Let the bulk element build the FaceElement, i.e. setup the pointers
226  // to its nodes (by referring to the appropriate nodes in the bulk
227  // element), etc.
228  bulk_el_pt->build_face_element(face_index,this);
229 
230 #ifdef PARANOID
231  {
232  //Check that the element is not a refineable 3d element
233  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
234 
235  //If it's three-d
236  if(elem_pt->dim()==3)
237  {
238  //Is it refineable
239  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
240  if(ref_el_pt!=0)
241  {
242  if (this->has_hanging_nodes())
243  {
244  throw OomphLibError(
245  "This flux element will not work correctly if nodes are hanging\n",
246  OOMPH_CURRENT_FUNCTION,
247  OOMPH_EXCEPTION_LOCATION);
248  }
249  }
250  }
251  }
252 #endif
253 
254  // Initialise the prescribed-flux function pointer to zero
255  Flux_fct_pt = 0;
256 
257  // Extract the dimension of the problem from the dimension of
258  // the first node
259  Dim = node_pt(0)->ndim();
260 
261  //Set up U_index_lin_wave. Initialise to zero, which probably won't change
262  //in most cases, oh well, the price we pay for generality
263  U_index_lin_wave = 0;
264 
265  //Cast to the appropriate LinearWaveEquation so that we can
266  //find the index at which the variable is stored
267  //We assume that the dimension of the full problem is the same
268  //as the dimension of the node, if this is not the case you will have
269  //to write custom elements, sorry
270  switch(Dim)
271  {
272  //One dimensional problem
273  case 1:
274  {
275  LinearWaveEquations<1>* eqn_pt =
276  dynamic_cast<LinearWaveEquations<1>*>(bulk_el_pt);
277  //If the cast has failed die
278  if(eqn_pt==0)
279  {
280  std::string error_string =
281  "Bulk element must inherit from LinearWaveEquations.";
282  error_string +=
283  "Nodes are one dimensional, but cannot cast the bulk element to\n";
284  error_string += "LinearWaveEquations<1>\n.";
285  error_string +=
286  "If you desire this functionality, you must implement it yourself\n";
287 
288  throw OomphLibError(
289  error_string,
290  OOMPH_CURRENT_FUNCTION,
291  OOMPH_EXCEPTION_LOCATION);
292  }
293  //Otherwise read out the value
294  else
295  {
296  //Read the index from the (cast) bulk element
298  }
299  }
300  break;
301 
302  //Two dimensional problem
303  case 2:
304  {
305  LinearWaveEquations<2>* eqn_pt =
306  dynamic_cast<LinearWaveEquations<2>*>(bulk_el_pt);
307  //If the cast has failed die
308  if(eqn_pt==0)
309  {
310  std::string error_string =
311  "Bulk element must inherit from LinearWaveEquations.";
312  error_string +=
313  "Nodes are two dimensional, but cannot cast the bulk element to\n";
314  error_string += "LinearWaveEquations<2>\n.";
315  error_string +=
316  "If you desire this functionality, you must implement it yourself\n";
317 
318  throw OomphLibError(
319  error_string,
320  OOMPH_CURRENT_FUNCTION,
321  OOMPH_EXCEPTION_LOCATION);
322  }
323  else
324  {
325  //Read the index from the (cast) bulk element.
327  }
328  }
329  break;
330 
331  //Three dimensional problem
332  case 3:
333  {
334  LinearWaveEquations<3>* eqn_pt =
335  dynamic_cast<LinearWaveEquations<3>*>(bulk_el_pt);
336  //If the cast has failed die
337  if(eqn_pt==0)
338  {
339  std::string error_string =
340  "Bulk element must inherit from LinearWaveEquations.";
341  error_string +=
342  "Nodes are three dimensional, but cannot cast the bulk element to\n";
343  error_string += "LinearWaveEquations<3>\n.";
344  error_string +=
345  "If you desire this functionality, you must implement it yourself\n";
346 
347  throw OomphLibError(
348  error_string,
349  OOMPH_CURRENT_FUNCTION,
350  OOMPH_EXCEPTION_LOCATION);
351 
352  }
353  else
354  {
355  //Read the index from the (cast) bulk element.
357  }
358  }
359  break;
360 
361  //Any other case is an error
362  default:
363  std::ostringstream error_stream;
364  error_stream << "Dimension of node is " << Dim
365  << ". It should be 1,2, or 3!" << std::endl;
366 
367  throw OomphLibError(
368  error_stream.str(),
369  OOMPH_CURRENT_FUNCTION,
370  OOMPH_EXCEPTION_LOCATION);
371  break;
372  }
373  }
374 
375 
376 
377 //===========================================================================
378 /// Compute the element's residual vector and the (zero) Jacobian matrix.
379 //===========================================================================
380 template<class ELEMENT>
383  Vector<double> &residuals, DenseMatrix<double> &jacobian,
384  unsigned flag)
385 {
386  //Find out how many nodes there are
387  const unsigned n_node = nnode();
388 
389  // Get continuous time from timestepper of first node
390  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
391 
392  //Set up memory for the shape and test functions
393  Shape psif(n_node), testf(n_node);
394 
395  //Set the value of n_intpt
396  const unsigned n_intpt = integral_pt()->nweight();
397 
398  //Set the Vector to hold local coordinates
399  Vector<double> s(Dim-1);
400 
401  //Integer to store the local equation and unknowns
402  int local_eqn=0;
403 
404  // Locally cache the index at which the variable is stored
405  const unsigned u_index_lin_wave = U_index_lin_wave;
406 
407  //Loop over the integration points
408  //--------------------------------
409  for(unsigned ipt=0;ipt<n_intpt;ipt++)
410  {
411 
412  //Assign values of s
413  for(unsigned i=0;i<(Dim-1);i++)
414  {
415  s[i] = integral_pt()->knot(ipt,i);
416  }
417 
418  //Get the integral weight
419  double w = integral_pt()->weight(ipt);
420 
421  //Find the shape and test functions and return the Jacobian
422  //of the mapping
423  double J = shape_and_test(s,psif,testf);
424 
425  //Premultiply the weights and the Jacobian
426  double W = w*J;
427 
428  //Need to find position to feed into flux function
429  Vector<double> interpolated_x(Dim);
430 
431  //Initialise to zero
432  for(unsigned i=0;i<Dim;i++)
433  {
434  interpolated_x[i] = 0.0;
435  }
436 
437  //Calculate velocities and derivatives
438  for(unsigned l=0;l<n_node;l++)
439  {
440  //Loop over velocity components
441  for(unsigned i=0;i<Dim;i++)
442  {
443  interpolated_x[i] += nodal_position(l,i)*psif[l];
444  }
445  }
446 
447  //Get the imposed flux
448  double flux;
449  get_flux(time,interpolated_x,flux);
450 
451  //Now add to the appropriate equations
452 
453  //Loop over the test functions
454  for(unsigned l=0;l<n_node;l++)
455  {
456  local_eqn = nodal_local_eqn(l,u_index_lin_wave);
457  /*IF it's not a boundary condition*/
458  if(local_eqn >= 0)
459  {
460  //Add the prescribed flux terms
461  residuals[local_eqn] -= flux*testf[l]*W;
462 
463  // Imposed traction doesn't depend upon the solution,
464  // --> the Jacobian is always zero, so no Jacobian
465  // terms are required
466  }
467  }
468  }
469 }
470 
471 
472 }
473 
474 #endif
LinearWavePrescribedFluxFctPt Flux_fct_pt
Function pointer to the (global) prescribed-flux function.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element residual vector.
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
virtual unsigned u_index_lin_wave() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
A general Finite Element class.
Definition: elements.h:1271
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: .
A class for elements that allow the imposition of an applied flux on the boundaries of LinearWave ele...
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual vector and its Jacobian matrix.
void operator=(const LinearWaveFluxElement &)
Broken assignment operator.
LinearWaveFluxElement(const LinearWaveFluxElement &dummy)
Broken copy constructor.
unsigned U_index_lin_wave
The index at which the unknown is stored at the nodes.
void get_flux(const double &time, const Vector< double > &x, double &flux)
static char t char * s
Definition: cfortran.h:572
LinearWaveFluxElement()
Broken empty constructor.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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
void output(std::ostream &outfile)
void fill_in_generic_residual_contribution_lin_wave_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute the element residual vector. flag=1(or 0): do (or don't) compute the Jacobian as well...
void output(FILE *file_pt, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
unsigned Dim
The spatial dimension of the problem.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
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 ...
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(* LinearWavePrescribedFluxFctPt)(const double &time, const Vector< double > &x, double &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector!
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
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
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
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 output(std::ostream &outfile, const unsigned &n_plot)
Output function – forward to broken version in FiniteElement until somebody decides what exactly they...
LinearWavePrescribedFluxFctPt & flux_fct_pt()
Access function for the prescribed-flux function pointer.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...