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