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