darcy_face_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 surface loads to
31 // the Darcy equations
32 
33 #ifndef OOMPH_DARCY_FACE_ELEMENTS_HEADER
34 #define OOMPH_DARCY_FACE_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 
42 // OOMPH-LIB headers
43 #include "../../../src/generic/Qelements.h"
44 
45 namespace oomph
46 {
47 
48 
49 
50 //=======================================================================
51 /// Namespace containing the zero pressure function for Darcy pressure
52 /// elements
53 //=======================================================================
54 namespace DarcyFaceElementHelper
55  {
56 
57  //=======================================================================
58  /// Default load function (zero pressure)
59  //=======================================================================
60  void Zero_pressure_fct(const double& time,
61  const Vector<double> &x,
62  const Vector<double>& N,
63  double &load)
64  {
65  load=0.0;
66  }
67 
68  }
69 
70 
71 //======================================================================
72 /// A class for elements that allow the imposition of an applied pressure
73 /// in the Darcy equations.
74 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
75 /// class and and thus, we can be generic enough without the need to have
76 /// a separate equations class.
77 //======================================================================
78 template <class ELEMENT>
79 class DarcyFaceElement : public virtual FaceGeometry<ELEMENT>,
80 public virtual FaceElement
81 {
82 protected:
83 
84  /// \short Pointer to an imposed pressure function. Arguments:
85  /// Eulerian coordinate; outer unit normal; applied pressure.
86  /// (Not all of the input arguments will be required for all specific load
87  /// functions but the list should cover all cases)
88  void (*Pressure_fct_pt)(const double &time,
89  const Vector<double> &x,
90  const Vector<double> &n,
91  double &result);
92 
93 
94  /// \short Get the pressure value: Pass number of integration point (dummy),
95  /// Eulerlian coordinate and normal vector and return the pressure
96  /// (not all of the input arguments will be required for all specific load
97  /// functions but the list should cover all cases). This function is virtual
98  /// so it can be overloaded for FSI.
99  virtual void get_pressure(const double &time,
100  const unsigned &intpt,
101  const Vector<double> &x,
102  const Vector<double> &n,
103  double &pressure)
104  {
105  Pressure_fct_pt(time,x,n,pressure);
106  }
107 
108 
109  /// \short Helper function that actually calculates the residuals
110  // This small level of indirection is required to avoid calling
111  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
112  // which causes all kinds of pain if overloading later on
114  Vector<double> &residuals);
115 
116 
117 public:
118 
119  /// \short Constructor, which takes a "bulk" element and the value of the
120  /// index and its limit
121  DarcyFaceElement(FiniteElement* const &element_pt,
122  const int &face_index) :
123  FaceGeometry<ELEMENT>(), FaceElement()
124  {
125 
126 #ifdef PARANOID
127  {
128  //Check that the element is not a refineable 3d element
129  ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
130  //If it's three-d
131  if(elem_pt->dim()==3)
132  {
133  //Is it refineable
134  RefineableElement* ref_el_pt=dynamic_cast<RefineableElement*>(elem_pt);
135  if(ref_el_pt!=0)
136  {
137  if (this->has_hanging_nodes())
138  {
139  throw OomphLibError(
140  "This flux element will not work correctly if nodes are hanging\n",
141  OOMPH_CURRENT_FUNCTION,
142  OOMPH_EXCEPTION_LOCATION);
143  }
144  }
145  }
146  }
147 #endif
148 
149  // Attach the geometrical information to the element. N.B. This function
150  // also assigns nbulk_value from the required_nvalue of the bulk element
151  element_pt->build_face_element(face_index,this);
152 
153  // Zero pressure
155  }
156 
157 
158  /// Reference to the pressure function pointer
159  void (* &pressure_fct_pt())(const double& time,
160  const Vector<double>& x,
161  const Vector<double>& n,
162  double &pressure)
163  {return Pressure_fct_pt;}
164 
165 
166  /// Return the residuals
168  {
170  }
171 
172 
173 
174  /// Fill in contribution from Jacobian
176  DenseMatrix<double> &jacobian)
177  {
178  // Call the residuals
180  }
181 
182  /// Specify the value of nodal zeta from the face geometry
183  /// \short The "global" intrinsic coordinate of the element when
184  /// viewed as part of a geometric object should be given by
185  /// the FaceElement representation, by default (needed to break
186  /// indeterminacy if bulk element is SolidElement)
187  double zeta_nodal(const unsigned &n, const unsigned &k,
188  const unsigned &i) const
189  {return FaceElement::zeta_nodal(n,k,i);}
190 
191  /// \short Output function
192  void output(std::ostream &outfile)
194 
195  /// \short Output function
196  void output(std::ostream &outfile, const unsigned &n_plot)
197  {
198  FaceGeometry<ELEMENT>::output(outfile,n_plot);
199  }
200 
201  /// \short C_style output function
202  void output(FILE* file_pt)
204 
205  /// \short C-style output function
206  void output(FILE* file_pt, const unsigned &n_plot)
207  {FaceGeometry<ELEMENT>::output(file_pt,n_plot);}
208 
209 
210  /// \short Compute pressure value at specified local coordinate
211  /// Should only be used for post-processing; ignores dependence
212  /// on integration point!
213  void pressure(const double& time,
214  const Vector<double>& s,
215  double &pressure);
216 
217 };
218 
219 ///////////////////////////////////////////////////////////////////////
220 ///////////////////////////////////////////////////////////////////////
221 ///////////////////////////////////////////////////////////////////////
222 
223 //=====================================================================
224 /// Compute pressure value at specified local coordinate
225 /// Should only be used for post-processing; ignores dependence
226 /// on integration point!
227 //=====================================================================
228 template<class ELEMENT>
230  const double& time, const Vector<double>& s, double &pressure)
231  {
232  unsigned n_dim = this->nodal_dimension();
233 
234  // Position vector
235  Vector<double> x(n_dim);
236  interpolated_x(s,x);
237 
238  // Outer unit normal
239  Vector<double> unit_normal(n_dim);
240  outer_unit_normal(s,unit_normal);
241 
242  // Dummy
243  unsigned ipt=0;
244 
245  // Pressure value
246  get_pressure(time,ipt,x,unit_normal,pressure);
247 
248  }
249 
250 
251 //=====================================================================
252 /// Return the residuals for the DarcyFaceElement equations
253 //=====================================================================
254 template<class ELEMENT>
257  Vector<double> &residuals)
258  {
259 
260  // Find out how many nodes there are
261  unsigned n_node = nnode();
262 
263  // Get continuous time from timestepper of first node
264  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
265 
266  #ifdef PARANOID
267  // Find out how many positional dofs there are
268  unsigned n_position_type = this->nnodal_position_type();
269  if(n_position_type != 1)
270  {
271  throw OomphLibError(
272  "Darcy equations are not yet implemented for more than one position type",
273  OOMPH_CURRENT_FUNCTION,
274  OOMPH_EXCEPTION_LOCATION);
275  }
276 #endif
277 
278  // Find out the dimension of the node
279  unsigned n_dim = this->nodal_dimension();
280 
281  // Set the pointer to the bulk element
282  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(bulk_element_pt());
283 
284  unsigned n_q_basis = bulk_el_pt->nq_basis();
285  unsigned n_q_basis_edge = bulk_el_pt->nq_basis_edge();
286 
287  // Integer to hold the local equation number
288  int local_eqn=0;
289 
290  // Set up memory for the shape functions
291  // Note that in this case, the number of lagrangian coordinates is always
292  // equal to the dimension of the nodes
293  Shape psi(n_node);
294  DShape dpsids(n_node,n_dim-1);
295 
296  Shape q_basis(n_q_basis,n_dim);
297 
298  // Set the value of n_intpt
299  unsigned n_intpt = integral_pt()->nweight();
300 
301  // Storage for the local coordinates
302  Vector<double> s_face(n_dim-1,0.0), s_bulk(n_dim,0.0);
303 
304  // Loop over the integration points
305  for(unsigned ipt=0;ipt<n_intpt;ipt++)
306  {
307  // Get the integral weight
308  double w = integral_pt()->weight(ipt);
309 
310  // Only need to call the local derivatives
311  dshape_local_at_knot(ipt,psi,dpsids);
312 
313  // Assign values of s in FaceElement and local coordinates in bulk element
314  for(unsigned i=0;i<n_dim-1;i++)
315  {
316  s_face[i] = integral_pt()->knot(ipt,i);
317  }
318 
319  s_bulk=local_coordinate_in_bulk(s_face);
320 
321  // Get the q basis at bulk local coordinate s_bulk, corresponding to
322  // face local coordinate s_face
323  bulk_el_pt->get_q_basis(s_bulk,q_basis);
324 
325  // Calculate the Eulerian and Lagrangian coordinates
326  Vector<double> interpolated_x(n_dim,0.0);
327 
328  // Also calculate the surface Vectors (derivatives wrt local coordinates)
329  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
330 
331  // Calculate displacements and derivatives
332  for(unsigned l=0;l<n_node;l++)
333  {
334  // Loop over directions
335  for(unsigned i=0;i<n_dim;i++)
336  {
337  // Calculate the Eulerian coords
338  const double x_local = nodal_position(l,i);
339  interpolated_x[i] += x_local*psi(l);
340 
341  // Loop over LOCAL derivative directions, to calculate the tangent(s)
342  for(unsigned j=0;j<n_dim-1;j++)
343  {
344  interpolated_A(j,i) += x_local*dpsids(l,j);
345  }
346  }
347  }
348 
349  // Now find the local metric tensor from the tangent vectors
350  DenseMatrix<double> A(n_dim-1);
351  for(unsigned i=0;i<n_dim-1;i++)
352  {
353  for(unsigned j=0;j<n_dim-1;j++)
354  {
355  // Initialise surface metric tensor to zero
356  A(i,j) = 0.0;
357 
358  // Take the dot product
359  for(unsigned k=0;k<n_dim;k++)
360  {
361  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
362  }
363  }
364  }
365 
366  // Get the outer unit normal
367  Vector<double> interpolated_normal(n_dim);
368  outer_unit_normal(ipt,interpolated_normal);
369 
370  // Find the determinant of the metric tensor
371  double Adet =0.0;
372  switch(n_dim)
373  {
374  case 2:
375  Adet = A(0,0);
376  break;
377  case 3:
378  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
379  break;
380  default:
381  throw
383  "Wrong dimension in DarcyFaceElement",
384  OOMPH_CURRENT_FUNCTION,
385  OOMPH_EXCEPTION_LOCATION);
386  }
387 
388  // Premultiply the weights and the square-root of the determinant of
389  // the metric tensor
390  double W = w*sqrt(Adet);
391 
392  // Now calculate the load
393  double pressure;
394  get_pressure(time,
395  ipt,
396  interpolated_x,
397  interpolated_normal,
398  pressure);
399 
400  // Loop over the q edge test functions only (internal basis functions
401  // have zero normal component on the boundary)
402  for(unsigned l=0;l<n_q_basis_edge;l++)
403  {
404  local_eqn = this->nodal_local_eqn(1,bulk_el_pt->q_edge_index(l));
405 
406  /*IF it's not a boundary condition*/
407  if(local_eqn >= 0)
408  {
409  // Loop over the displacement components
410  for(unsigned i=0;i<n_dim;i++)
411  {
412  // Add the loading terms to the residuals
413  residuals[local_eqn] +=
414  pressure*q_basis(l,i)*interpolated_normal[i]*W;
415  }
416  } // End of if not boundary condition
417  } //End of loop over shape functions
418  } //End of loop over integration points
419  }
420 
421 }
422 
423 #endif
void(* Pressure_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, double &result)
Pointer to an imposed pressure function. Arguments: Eulerian coordinate; outer unit normal; applied p...
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 ...
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 Zero_pressure_fct(const double &time, const Vector< double > &x, const Vector< double > &N, double &load)
Default load function (zero pressure)
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, double &pressure) pressure_fct_pt()
Reference to the pressure function pointer.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
static char t char * s
Definition: cfortran.h:572
virtual void get_pressure(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, double &pressure)
Get the pressure value: Pass number of integration point (dummy), Eulerlian coordinate and normal vec...
void pressure(const double &time, const Vector< double > &s, double &pressure)
Compute pressure value at specified local coordinate Should only be used for post-processing; ignores...
void output(std::ostream &outfile)
Output with default number of plot points.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
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 output(FILE *file_pt)
C_style output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
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
void fill_in_contribution_to_residuals_darcy_face(Vector< double > &residuals)
Helper function that actually calculates the residuals.
DarcyFaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void output(std::ostream &outfile)
Output function.