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$
7 //LIC//
8 //LIC// $LastChangedDate$
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
124  void output(std::ostream &outfile)
125  {
126  const unsigned n_plot=5;
127  output(outfile,n_plot);
128  }
129 
130  /// \short Output function
131  void output(std::ostream &outfile, const unsigned &nplot)
132  {
133 
134  // Dimension of element
135  unsigned el_dim=dim();
136 
137  //Vector of local coordinates
138  Vector<double> s(el_dim);
139 
140  // Tecplot header info
141  outfile << tecplot_zone_string(nplot);
142 
143  // Loop over plot points
144  unsigned num_plot_points=nplot_points(nplot);
145  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
146  {
147 
148  // Get local coordinates of plot point
149  get_s_plot(iplot,nplot,s);
150 
151  Vector<double> x(el_dim+1);
152  for(unsigned i=0;i<el_dim+1;i++)
153  {
154  x[i]=interpolated_x(s,i);
155  outfile << x[i] << " ";
156  }
157  double flux=0.0;
158  get_flux(x,flux);
159  outfile << flux << std::endl;
160  }
161 
162  // Write tecplot footer (e.g. FE connectivity lists)
163  write_tecplot_zone_footer(outfile,nplot);
164 
165  }
166 
167 
168  /// C-style output function -- forward to broken version in FiniteElement
169  /// until somebody decides what exactly they want to plot here...
170  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
171 
172  /// \short C-style output function -- forward to broken version in
173  /// FiniteElement until somebody decides what exactly they want to plot
174  /// here...
175  void output(FILE* file_pt, const unsigned &n_plot)
176  {FiniteElement::output(file_pt,n_plot);}
177 
178 
179 protected:
180 
181  /// \short Function to compute the shape and test functions and to return
182  /// the Jacobian of mapping between local and global (Eulerian)
183  /// coordinates
184  inline double shape_and_test(const Vector<double> &s, Shape &psi, Shape &test)
185  const
186  {
187  //Find number of nodes
188  unsigned n_node = nnode();
189 
190  //Get the shape functions
191  shape(s,psi);
192 
193  //Set the test functions to be the same as the shape functions
194  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
195 
196  //Return the value of the jacobian
197  return J_eulerian(s);
198  }
199 
200 
201  /// \short Function to compute the shape and test functions and to return
202  /// the Jacobian of mapping between local and global (Eulerian)
203  /// coordinates
204  inline double shape_and_test_at_knot(const unsigned &ipt,
205  Shape &psi, Shape &test)
206  const
207  {
208  //Find number of nodes
209  unsigned n_node = nnode();
210 
211  //Get the shape functions
212  shape_at_knot(ipt,psi);
213 
214  //Set the test functions to be the same as the shape functions
215  for(unsigned i=0;i<n_node;i++) {test[i] = psi[i];}
216 
217  //Return the value of the jacobian
218  return J_eulerian_at_knot(ipt);
219  }
220 
221 
222  /// Function to calculate the prescribed flux at a given spatial
223  /// position
224  void get_flux(const Vector<double>& x, double& flux)
225  {
226  //If the function pointer is zero return zero
227  if(Flux_fct_pt == 0)
228  {
229  flux=0.0;
230  }
231  //Otherwise call the function
232  else
233  {
234  (*Flux_fct_pt)(x,flux);
235  }
236  }
237 
238 private:
239 
240 
241  /// \short Add the element's contribution to its residual vector.
242  /// flag=1(or 0): do (or don't) compute the contribution to the
243  /// Jacobian as well.
245  Vector<double> &residuals, DenseMatrix<double> &jacobian,
246  const unsigned& flag);
247 
248 
249  /// Function pointer to the (global) prescribed-flux function
251 
252  ///The spatial dimension of the problem
253  unsigned Dim;
254 
255  ///The index at which the unknown is stored at the nodes
256  unsigned U_index_poisson;
257 
258 
259 };
260 
261 //////////////////////////////////////////////////////////////////////
262 //////////////////////////////////////////////////////////////////////
263 //////////////////////////////////////////////////////////////////////
264 
265 
266 
267 //===========================================================================
268 /// Constructor, takes the pointer to the "bulk" element, the
269 /// index of the fixed local coordinate and its value represented
270 /// by an integer (+/- 1), indicating that the face is located
271 /// at the max. or min. value of the "fixed" local coordinate
272 /// in the bulk element.
273 //===========================================================================
274 template<class ELEMENT>
276 PoissonFluxElement(FiniteElement* const &bulk_el_pt,
277  const int &face_index) :
278  FaceGeometry<ELEMENT>(), FaceElement()
279  {
280  // Let the bulk element build the FaceElement, i.e. setup the pointers
281  // to its nodes (by referring to the appropriate nodes in the bulk
282  // element), etc.
283  bulk_el_pt->build_face_element(face_index,this);
284 
285 #ifdef PARANOID
286  {
287  //Check that the element is not a refineable 3d element
288  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(bulk_el_pt);
289  //If it's three-d
290  if(elem_pt->dim()==3)
291  {
292  //Is it refineable
293  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
294  if(ref_el_pt!=0)
295  {
296  if (this->has_hanging_nodes())
297  {
298  throw OomphLibError(
299  "This flux element will not work correctly if nodes are hanging\n",
300  OOMPH_CURRENT_FUNCTION,
301  OOMPH_EXCEPTION_LOCATION);
302  }
303  }
304  }
305  }
306 #endif
307 
308  // Initialise the prescribed-flux function pointer to zero
309  Flux_fct_pt = 0;
310 
311  // Extract the dimension of the problem from the dimension of
312  // the first node
313  Dim = this->node_pt(0)->ndim();
314 
315  //Set up U_index_poisson. Initialise to zero, which probably won't change
316  //in most cases, oh well, the price we pay for generality
317  U_index_poisson = 0;
318 
319  //Cast to the appropriate PoissonEquation so that we can
320  //find the index at which the variable is stored
321  //We assume that the dimension of the full problem is the same
322  //as the dimension of the node, if this is not the case you will have
323  //to write custom elements, sorry
324  switch(Dim)
325  {
326  //One dimensional problem
327  case 1:
328  {
329  PoissonEquations<1>* eqn_pt =
330  dynamic_cast<PoissonEquations<1>*>(bulk_el_pt);
331  //If the cast has failed die
332  if(eqn_pt==0)
333  {
334  std::string error_string =
335  "Bulk element must inherit from PoissonEquations.";
336  error_string +=
337  "Nodes are one dimensional, but cannot cast the bulk element to\n";
338  error_string += "PoissonEquations<1>\n.";
339  error_string +=
340  "If you desire this functionality, you must implement it yourself\n";
341 
342  throw OomphLibError(error_string,
343  OOMPH_CURRENT_FUNCTION,
344  OOMPH_EXCEPTION_LOCATION);
345  }
346  //Otherwise read out the value
347  else
348  {
349  //Read the index from the (cast) bulk element
350  U_index_poisson = eqn_pt->u_index_poisson();
351  }
352  }
353  break;
354 
355  //Two dimensional problem
356  case 2:
357  {
358  PoissonEquations<2>* eqn_pt =
359  dynamic_cast<PoissonEquations<2>*>(bulk_el_pt);
360  //If the cast has failed die
361  if(eqn_pt==0)
362  {
363  std::string error_string =
364  "Bulk element must inherit from PoissonEquations.";
365  error_string +=
366  "Nodes are two dimensional, but cannot cast the bulk element to\n";
367  error_string += "PoissonEquations<2>\n.";
368  error_string +=
369  "If you desire this functionality, you must implement it yourself\n";
370 
371  throw OomphLibError(error_string,
372  OOMPH_CURRENT_FUNCTION,
373  OOMPH_EXCEPTION_LOCATION);
374  }
375  else
376  {
377  //Read the index from the (cast) bulk element.
378  U_index_poisson = eqn_pt->u_index_poisson();
379  }
380  }
381  break;
382 
383  //Three dimensional problem
384  case 3:
385  {
386  PoissonEquations<3>* eqn_pt =
387  dynamic_cast<PoissonEquations<3>*>(bulk_el_pt);
388  //If the cast has failed die
389  if(eqn_pt==0)
390  {
391  std::string error_string =
392  "Bulk element must inherit from PoissonEquations.";
393  error_string +=
394  "Nodes are three dimensional, but cannot cast the bulk element to\n";
395  error_string += "PoissonEquations<3>\n.";
396  error_string +=
397  "If you desire this functionality, you must implement it yourself\n";
398 
399  throw OomphLibError(error_string,
400  OOMPH_CURRENT_FUNCTION,
401  OOMPH_EXCEPTION_LOCATION);
402 
403  }
404  else
405  {
406  //Read the index from the (cast) bulk element.
407  U_index_poisson = eqn_pt->u_index_poisson();
408  }
409  }
410  break;
411 
412  //Any other case is an error
413  default:
414  std::ostringstream error_stream;
415  error_stream << "Dimension of node is " << Dim
416  << ". It should be 1,2, or 3!" << std::endl;
417 
418  throw OomphLibError(error_stream.str(),
419  OOMPH_CURRENT_FUNCTION,
420  OOMPH_EXCEPTION_LOCATION);
421  break;
422  }
423  }
424 
425 
426 //===========================================================================
427 /// Compute the element's residual vector and the (zero) Jacobian matrix.
428 //===========================================================================
429 template<class ELEMENT>
432  Vector<double> &residuals, DenseMatrix<double> &jacobian,
433  const unsigned& flag)
434 {
435  //Find out how many nodes there are
436  const unsigned n_node = nnode();
437 
438  //Set up memory for the shape and test functions
439  Shape psif(n_node), testf(n_node);
440 
441  //Set the value of Nintpt
442  const unsigned n_intpt = integral_pt()->nweight();
443 
444  //Set the Vector to hold local coordinates
445  Vector<double> s(Dim-1);
446 
447  //Integers to hold the local equation and unknown numbers
448  int local_eqn=0;
449 
450  // Locally cache the index at which the variable is stored
451  const unsigned u_index_poisson = U_index_poisson;
452 
453  //Loop over the integration points
454  //--------------------------------
455  for(unsigned ipt=0;ipt<n_intpt;ipt++)
456  {
457 
458  //Assign values of s
459  for(unsigned i=0;i<(Dim-1);i++) {s[i] = integral_pt()->knot(ipt,i);}
460 
461  //Get the integral weight
462  double w = integral_pt()->weight(ipt);
463 
464  //Find the shape and test functions and return the Jacobian
465  //of the mapping
466  double J = shape_and_test(s,psif,testf);
467 
468  //Premultiply the weights and the Jacobian
469  double W = w*J;
470 
471  //Need to find position to feed into flux function, initialise to zero
473 
474  //Calculate velocities and derivatives
475  for(unsigned l=0;l<n_node;l++)
476  {
477  //Loop over velocity components
478  for(unsigned i=0;i<Dim;i++)
479  {
480  interpolated_x[i] += nodal_position(l,i)*psif[l];
481  }
482  }
483 
484  //Get the imposed flux
485  double flux;
486  get_flux(interpolated_x,flux);
487 
488  //Now add to the appropriate equations
489 
490  //Loop over the test functions
491  for(unsigned l=0;l<n_node;l++)
492  {
493  local_eqn = nodal_local_eqn(l,u_index_poisson);
494  /*IF it's not a boundary condition*/
495  if(local_eqn >= 0)
496  {
497  //Add the prescribed flux terms
498  residuals[local_eqn] -= flux*testf[l]*W;
499 
500  // Imposed traction doesn't depend upon the solution,
501  // --> the Jacobian is always zero, so no Jacobian
502  // terms are required
503  }
504  }
505  }
506 }
507 
508 
509 }
510 
511 #endif
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2990
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:5201
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exa...
void fill_in_generic_residual_contribution_poisson_flux(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add the element&#39;s contribution to its residual vector. flag=1(or 0): do (or don&#39;t) compute the contri...
bool has_hanging_nodes() const
Return boolean to indicate if any of the element&#39;s nodes are geometrically hanging.
Definition: elements.h:2356
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
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 get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2978
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:2884
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2227
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:4412
A general Finite Element class.
Definition: elements.h:1274
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:4292
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4322
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&#39;s contribution to its residual vector.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
void output(std::ostream &outfile)
Output function.
void(* PoissonPrescribedFluxFctPt)(const Vector< double > &x, double &flux)
Function pointer to the prescribed-flux function fct(x,f(x)) – x is a Vector!
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:3017
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:5113
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...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1908
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:3160
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3003
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 ...
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:2109
unsigned Dim
The spatial dimension of the problem.
PoissonPrescribedFluxFctPt & flux_fct_pt()
Access function for the prescribed-flux function pointer.
void operator=(const PoissonFluxElement &)
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:2482
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void output(std::ostream &outfile, const unsigned &nplot)
Output function.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
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...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn&#39;t been defined.
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:5002
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2146
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&#39;s contribution to its residual vector and its Jacobian matrix.
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:1386