poroelasticity_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_POROELASITICTY_FACE_ELEMENTS_HEADER
34 #define OOMPH_POROELASITICTY_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 "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 PoroelasticityFaceElementHelper
55  {
56 
57  //=======================================================================
58  /// Default load function (zero traction)
59  //=======================================================================
60  void Zero_traction_fct(const double& time,
61  const Vector<double> &x,
62  const Vector<double>& N,
63  Vector<double>& load)
64  {
65  unsigned n_dim=load.size();
66  for (unsigned i=0;i<n_dim;i++) {load[i]=0.0;}
67  }
68 
69  //=======================================================================
70  /// Default load function (zero pressure)
71  //=======================================================================
72  void Zero_pressure_fct(const double& time,
73  const Vector<double> &x,
74  const Vector<double>& N,
75  double &load)
76  {
77  load=0.0;
78  }
79 
80  }
81 
82 
83 //======================================================================
84 /// A class for elements that allow the imposition of an applied pressure
85 /// in the Darcy equations.
86 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
87 /// class and and thus, we can be generic enough without the need to have
88 /// a separate equations class.
89 //======================================================================
90 template <class ELEMENT>
91 class PoroelasticityFaceElement : public virtual FaceGeometry<ELEMENT>,
92 public virtual FaceElement
93 {
94 protected:
95 
96  /// pointer to the bulk element this face element is attached to
97  ELEMENT *Element_pt;
98 
99  /// \short Pointer to an imposed traction function. Arguments:
100  /// Eulerian coordinate; outer unit normal; applied traction.
101  /// (Not all of the input arguments will be required for all specific load
102  /// functions but the list should cover all cases)
103  void (*Traction_fct_pt)(const double& time,
104  const Vector<double> &x,
105  const Vector<double> &n,
106  Vector<double> &result);
107 
108  /// \short Pointer to an imposed pressure function. Arguments:
109  /// Eulerian coordinate; outer unit normal; applied pressure.
110  /// (Not all of the input arguments will be required for all specific load
111  /// functions but the list should cover all cases)
112  void (*Pressure_fct_pt)(const double& time,
113  const Vector<double> &x,
114  const Vector<double> &n,
115  double &result);
116 
117 
118  /// \short Get the traction vector: Pass number of integration point (dummy),
119  /// Eulerlian coordinate and normal vector and return the pressure
120  /// (not all of the input arguments will be required for all specific load
121  /// functions but the list should cover all cases). This function is virtual
122  /// so it can be overloaded for FSI.
123  virtual void get_traction(const double& time,
124  const unsigned& intpt,
125  const Vector<double>& x,
126  const Vector<double>& n,
128  {
129  Traction_fct_pt(time,x,n,traction);
130  }
131 
132  /// \short Get the pressure value: Pass number of integration point (dummy),
133  /// Eulerlian coordinate and normal vector and return the pressure
134  /// (not all of the input arguments will be required for all specific load
135  /// functions but the list should cover all cases). This function is virtual
136  /// so it can be overloaded for FSI.
137  virtual void get_pressure(const double& time,
138  const unsigned& intpt,
139  const Vector<double>& x,
140  const Vector<double>& n,
141  double &pressure)
142  {
143  Pressure_fct_pt(time,x,n,pressure);
144  }
145 
146 
147  /// \short Helper function that actually calculates the residuals
148  // This small level of indirection is required to avoid calling
149  // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
150  // which causes all kinds of pain if overloading later on
152  Vector<double> &residuals);
153 
154 
155 public:
156 
157  /// \short Constructor, which takes a "bulk" element and the value of the
158  /// index and its limit
160  const int &face_index) :
161  FaceGeometry<ELEMENT>(), FaceElement()
162  {
163 #ifdef PARANOID
164  {
165  // Check that the element is not a refineable 3d element
166  ELEMENT* elem_pt = new ELEMENT;
167  // If it's three-d
168  if(elem_pt->dim()==3)
169  {
170  // Is it refineable
171  if(dynamic_cast<RefineableElement*>(elem_pt))
172  {
173  // Issue a warning
175  "This flux element will not work correctly if nodes are hanging\n",
176  OOMPH_CURRENT_FUNCTION,
177  OOMPH_EXCEPTION_LOCATION);
178  }
179  }
180  }
181 #endif
182 
183  // Set the pointer to the bulk element
184  Element_pt=dynamic_cast<ELEMENT*>(element_pt);
185 
186  // Attach the geometrical information to the element. N.B. This function
187  // also assigns nbulk_value from the required_nvalue of the bulk element
188  element_pt->build_face_element(face_index,this);
189 
190  // Zero traction
192 
193  // Zero pressure
195  }
196 
197 
198  /// Reference to the traction function pointer
199  void (* &traction_fct_pt())(const double& time,
200  const Vector<double>& x,
201  const Vector<double>& n,
203  {return Traction_fct_pt;}
204 
205  /// Reference to the pressure function pointer
206  void (* &pressure_fct_pt())(const double& time,
207  const Vector<double>& x,
208  const Vector<double>& n,
209  double &pressure)
210  {return Pressure_fct_pt;}
211 
212 
213  /// Return the residuals
215  {
217  }
218 
219 
220 
221  /// Fill in contribution from Jacobian
223  DenseMatrix<double> &jacobian)
224  {
225  // Call the residuals
227  }
228 
229  /// Specify the value of nodal zeta from the face geometry
230  /// \short The "global" intrinsic coordinate of the element when
231  /// viewed as part of a geometric object should be given by
232  /// the FaceElement representation, by default (needed to break
233  /// indeterminacy if bulk element is SolidElement)
234  double zeta_nodal(const unsigned &n, const unsigned &k,
235  const unsigned &i) const
236  {return FaceElement::zeta_nodal(n,k,i);}
237 
238  /// \short Output function
239  void output(std::ostream &outfile)
240  {FiniteElement::output(outfile);}
241 
242  /// \short Output function
243  void output(std::ostream &outfile, const unsigned &n_plot)
244  {FiniteElement::output(outfile,n_plot);}
245 
246  /// \short C_style output function
247  void output(FILE* file_pt)
248  {FiniteElement::output(file_pt);}
249 
250  /// \short C-style output function
251  void output(FILE* file_pt, const unsigned &n_plot)
252  {FiniteElement::output(file_pt,n_plot);}
253 
254 
255  /// \short Compute traction vector at specified local coordinate
256  /// Should only be used for post-processing; ignores dependence
257  /// on integration point!
258  void traction(const double& time,
259  const Vector<double>& s,
261 
262  /// \short Compute pressure value at specified local coordinate
263  /// Should only be used for post-processing; ignores dependence
264  /// on integration point!
265  void pressure(const double& time,
266  const Vector<double>& s,
267  double &pressure);
268 
269 };
270 
271 ///////////////////////////////////////////////////////////////////////
272 ///////////////////////////////////////////////////////////////////////
273 ///////////////////////////////////////////////////////////////////////
274 
275 //=====================================================================
276 /// Compute traction vector at specified local coordinate
277 /// Should only be used for post-processing; ignores dependence
278 /// on integration point!
279 //=====================================================================
280 template<class ELEMENT>
282  const double& time, const Vector<double>& s, Vector<double> &traction)
283  {
284  unsigned n_dim = this->nodal_dimension();
285 
286  // Position vector
287  Vector<double> x(n_dim);
288  interpolated_x(s,x);
289 
290  // Outer unit normal
291  Vector<double> unit_normal(n_dim);
292  outer_unit_normal(s,unit_normal);
293 
294  // Dummy
295  unsigned ipt=0;
296 
297  // Pressure value
298  get_traction(time,ipt,x,unit_normal,traction);
299 
300  }
301 //=====================================================================
302 /// Compute pressure value at specified local coordinate
303 /// Should only be used for post-processing; ignores dependence
304 /// on integration point!
305 //=====================================================================
306 template<class ELEMENT>
308  const double& time, const Vector<double>& s, double &pressure)
309  {
310  unsigned n_dim = this->nodal_dimension();
311 
312  // Position vector
313  Vector<double> x(n_dim);
314  interpolated_x(s,x);
315 
316  // Outer unit normal
317  Vector<double> unit_normal(n_dim);
318  outer_unit_normal(s,unit_normal);
319 
320  // Dummy
321  unsigned ipt=0;
322 
323  // Pressure value
324  get_pressure(time,ipt,x,unit_normal,pressure);
325 
326  }
327 
328 
329 //=====================================================================
330 /// Return the residuals for the PoroelasticityFaceElement equations
331 //=====================================================================
332 template<class ELEMENT>
335  Vector<double> &residuals)
336  {
337 
338  // Find out how many nodes there are
339  unsigned n_node = nnode();
340 
341  // Get continuous time from timestepper of first node
342  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
343 
344  #ifdef PARANOID
345  // Find out how many positional dofs there are
346  unsigned n_position_type = this->nnodal_position_type();
347  if(n_position_type != 1)
348  {
349  throw OomphLibError(
350  "Poroelasticity equations are not yet implemented for more than one position type",
351  OOMPH_CURRENT_FUNCTION,
352  OOMPH_EXCEPTION_LOCATION);
353  }
354 #endif
355 
356  // Find out the dimension of the node
357  unsigned n_dim = this->nodal_dimension();
358 
359  unsigned n_q_basis = Element_pt->nq_basis();
360  unsigned n_q_basis_edge = Element_pt->nq_basis_edge();
361 
362  // Integer to hold the local equation number
363  int local_eqn=0;
364 
365  // Set up memory for the shape functions
366  // Note that in this case, the number of lagrangian coordinates is always
367  // equal to the dimension of the nodes
368  Shape psi(n_node);
369  DShape dpsids(n_node,n_dim-1);
370 
371  Shape q_basis(n_q_basis,n_dim);
372 
373  // Set the value of n_intpt
374  unsigned n_intpt = integral_pt()->nweight();
375 
376  // Storage for the local coordinates
377  Vector<double> s_face(n_dim-1,0.0), s_bulk(n_dim,0.0);
378 
379  //mjr TODO enable if/when eta is added to Poroelasticity elements
380  /*double eta = Element_pt->eta();*/
381 
382  // Loop over the integration points
383  for(unsigned ipt=0;ipt<n_intpt;ipt++)
384  {
385  // Get the integral weight
386  double w = integral_pt()->weight(ipt);
387 
388  // Only need to call the local derivatives
389  dshape_local_at_knot(ipt,psi,dpsids);
390 
391  // Assign values of s in FaceElement and local coordinates in bulk element
392  for(unsigned i=0;i<n_dim-1;i++)
393  {
394  s_face[i] = integral_pt()->knot(ipt,i);
395  }
396 
397  s_bulk=local_coordinate_in_bulk(s_face);
398 
399  // Get the q basis at bulk local coordinate s_bulk, corresponding to face local
400  // coordinate s_face
401  Element_pt->get_q_basis(s_bulk,q_basis);
402 
403  // Calculate the Eulerian and Lagrangian coordinates
404  Vector<double> interpolated_x(n_dim,0.0);
405 
406  // Also calculate the surface Vectors (derivatives wrt local coordinates)
407  DenseMatrix<double> interpolated_A(n_dim-1,n_dim,0.0);
408 
409  // Calculate displacements and derivatives
410  for(unsigned l=0;l<n_node;l++)
411  {
412  // Loop over directions
413  for(unsigned i=0;i<n_dim;i++)
414  {
415  // Calculate the Eulerian coords
416  const double x_local = nodal_position(l,i);
417  interpolated_x[i] += x_local*psi(l);
418 
419  // Loop over LOCAL derivative directions, to calculate the tangent(s)
420  for(unsigned j=0;j<n_dim-1;j++)
421  {
422  interpolated_A(j,i) += x_local*dpsids(l,j);
423  }
424  }
425  }
426 
427  // Now find the local metric tensor from the tangent vectors
428  DenseMatrix<double> A(n_dim-1);
429  for(unsigned i=0;i<n_dim-1;i++)
430  {
431  for(unsigned j=0;j<n_dim-1;j++)
432  {
433  // Initialise surface metric tensor to zero
434  A(i,j) = 0.0;
435 
436  // Take the dot product
437  for(unsigned k=0;k<n_dim;k++)
438  {
439  A(i,j) += interpolated_A(i,k)*interpolated_A(j,k);
440  }
441  }
442  }
443 
444  // Get the outer unit normal
445  Vector<double> interpolated_normal(n_dim);
446  outer_unit_normal(ipt,interpolated_normal);
447 
448  // Find the determinant of the metric tensor
449  double Adet =0.0;
450  switch(n_dim)
451  {
452  case 2:
453  Adet = A(0,0);
454  break;
455  case 3:
456  Adet = A(0,0)*A(1,1) - A(0,1)*A(1,0);
457  break;
458  default:
459  throw
461  "Wrong dimension in PoroelasticityFaceElement",
462  OOMPH_CURRENT_FUNCTION,
463  OOMPH_EXCEPTION_LOCATION);
464  }
465 
466  // Premultiply the weights and the square-root of the determinant of
467  // the metric tensor
468  double W = w*sqrt(Adet);
469 
470  //Now calculate the traction load
471  Vector<double> traction(n_dim);
472  get_traction(time,
473  ipt,
474  interpolated_x,
475  interpolated_normal,
476  traction);
477 
478  // Now calculate the load
479  double pressure;
480  get_pressure(time,
481  ipt,
482  interpolated_x,
483  interpolated_normal,
484  pressure);
485 
486  //Loop over the test functions, nodes of the element
487  for(unsigned l=0;l<n_node;l++)
488  {
489  //Loop over the displacement components
490  for(unsigned i=0;i<n_dim;i++)
491  {
492  local_eqn = this->nodal_local_eqn(l,Element_pt->u_index(i));
493  /*IF it's not a boundary condition*/
494  if(local_eqn >= 0)
495  {
496  //Add the traction loading terms to the residuals
497  residuals[local_eqn] -= traction[i]*psi(l)*W;
498  } //End of if not boundary condition
499  }
500  } //End of loop over shape functions
501 
502  // Loop over the q edge test functions only (internal basis functions
503  // have zero normal component on the boundary)
504  for(unsigned l=0;l<n_q_basis_edge;l++)
505  {
506  local_eqn = this->nodal_local_eqn(1,Element_pt->q_edge_index(l));
507 
508  /*IF it's not a boundary condition*/
509  if(local_eqn >= 0)
510  {
511  // Loop over the displacement components
512  for(unsigned i=0;i<n_dim;i++)
513  {
514  // Add the loading terms to the residuals
515  residuals[local_eqn] +=
516  pressure*q_basis(l,i)*interpolated_normal[i]*W;
517  }
518  } // End of if not boundary condition
519  } //End of loop over shape functions
520  } //End of loop over integration points
521  }
522 
523 }
524 
525 #endif
void Zero_traction_fct(const double &time, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
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(FILE *file_pt)
C_style output function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
cstr elem_len * i
Definition: cfortran.h:607
ELEMENT * Element_pt
pointer to the bulk element this face element is attached to
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
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...
A general Finite Element class.
Definition: elements.h:1271
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(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function. Arguments: Eulerian coordinate; outer unit normal; applied t...
virtual void get_traction(const double &time, const unsigned &intpt, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Get the traction vector: Pass number of integration point (dummy), Eulerlian coordinate and normal ve...
static char t char * s
Definition: cfortran.h:572
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void output(std::ostream &outfile)
Output function.
void Zero_pressure_fct(const double &time, const Vector< double > &x, const Vector< double > &N, double &load)
Default load function (zero pressure)
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, double &pressure) pressure_fct_pt()
Reference to the pressure function pointer.
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_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void fill_in_contribution_to_residuals_darcy_face(Vector< double > &residuals)
Helper function that actually calculates the residuals.
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...
PoroelasticityFaceElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction 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 ...
void traction(const double &time, const Vector< double > &s, Vector< double > &traction)
Compute traction vector at specified local coordinate Should only be used for post-processing; ignore...