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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for elements that are used to apply surface loads to
27// the Darcy equations
28
29#ifndef OOMPH_POROELASITICTY_FACE_ELEMENTS_HEADER
30#define OOMPH_POROELASITICTY_FACE_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37
38// OOMPH-LIB headers
39#include "generic/Qelements.h"
40
41namespace oomph
42{
43 //=======================================================================
44 /// Namespace containing the zero pressure function for Darcy pressure
45 /// elements
46 //=======================================================================
47 namespace PoroelasticityFaceElementHelper
48 {
49 //=======================================================================
50 /// Default load function (zero traction)
51 //=======================================================================
52 void Zero_traction_fct(const double& time,
53 const Vector<double>& x,
54 const Vector<double>& N,
55 Vector<double>& load)
56 {
57 unsigned n_dim = load.size();
58 for (unsigned i = 0; i < n_dim; i++)
59 {
60 load[i] = 0.0;
61 }
62 }
63
64 //=======================================================================
65 /// Default load function (zero pressure)
66 //=======================================================================
67 void Zero_pressure_fct(const double& time,
68 const Vector<double>& x,
69 const Vector<double>& N,
70 double& load)
71 {
72 load = 0.0;
73 }
74
75 } // namespace PoroelasticityFaceElementHelper
76
77
78 //======================================================================
79 /// A class for elements that allow the imposition of an applied pressure
80 /// in the Darcy equations.
81 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
82 /// class and and thus, we can be generic enough without the need to have
83 /// a separate equations class.
84 //======================================================================
85 template<class ELEMENT>
86 class PoroelasticityFaceElement : public virtual FaceGeometry<ELEMENT>,
87 public virtual FaceElement
88 {
89 protected:
90 /// pointer to the bulk element this face element is attached to
91 ELEMENT* Element_pt;
92
93 /// Pointer to an imposed traction function. Arguments:
94 /// Eulerian coordinate; outer unit normal; applied traction.
95 /// (Not all of the input arguments will be required for all specific load
96 /// functions but the list should cover all cases)
97 void (*Traction_fct_pt)(const double& time,
98 const Vector<double>& x,
99 const Vector<double>& n,
100 Vector<double>& result);
101
102 /// Pointer to an imposed pressure function. Arguments:
103 /// Eulerian coordinate; outer unit normal; applied pressure.
104 /// (Not all of the input arguments will be required for all specific load
105 /// functions but the list should cover all cases)
106 void (*Pressure_fct_pt)(const double& time,
107 const Vector<double>& x,
108 const Vector<double>& n,
109 double& result);
110
111
112 /// Get the traction vector: Pass number of integration point
113 /// (dummy), Eulerlian coordinate and normal vector and return the pressure
114 /// (not all of the input arguments will be required for all specific load
115 /// functions but the list should cover all cases). This function is virtual
116 /// so it can be overloaded for FSI.
117 virtual void get_traction(const double& time,
118 const unsigned& intpt,
119 const Vector<double>& x,
120 const Vector<double>& n,
122 {
123 Traction_fct_pt(time, x, n, traction);
124 }
125
126 /// Get the pressure value: Pass number of integration point (dummy),
127 /// Eulerlian coordinate and normal vector and return the pressure
128 /// (not all of the input arguments will be required for all specific load
129 /// functions but the list should cover all cases). This function is virtual
130 /// so it can be overloaded for FSI.
131 virtual void get_pressure(const double& time,
132 const unsigned& intpt,
133 const Vector<double>& x,
134 const Vector<double>& n,
135 double& pressure)
136 {
137 Pressure_fct_pt(time, x, n, pressure);
138 }
139
140
141 /// Helper function that actually calculates the residuals
142 // This small level of indirection is required to avoid calling
143 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
144 // which causes all kinds of pain if overloading later on
146 Vector<double>& residuals);
147
148
149 public:
150 /// Constructor, which takes a "bulk" element and the value of the
151 /// index and its limit
153 const int& face_index)
154 : FaceGeometry<ELEMENT>(), FaceElement()
155 {
156#ifdef PARANOID
157 {
158 // Check that the element is not a refineable 3d element
159 ELEMENT* elem_pt = new ELEMENT;
160 // If it's three-d
161 if (elem_pt->dim() == 3)
162 {
163 // Is it refineable
164 if (dynamic_cast<RefineableElement*>(elem_pt))
165 {
166 // Issue a warning
167 OomphLibWarning("This flux element will not work correctly if "
168 "nodes are hanging\n",
169 OOMPH_CURRENT_FUNCTION,
170 OOMPH_EXCEPTION_LOCATION);
171 }
172 }
173 }
174#endif
175
176 // Set the pointer to the bulk element
177 Element_pt = dynamic_cast<ELEMENT*>(element_pt);
178
179 // Attach the geometrical information to the element. N.B. This function
180 // also assigns nbulk_value from the required_nvalue of the bulk element
181 element_pt->build_face_element(face_index, this);
182
183 // Zero traction
185
186 // Zero pressure
188 }
189
190
191 /// Reference to the traction function pointer
192 void (*&traction_fct_pt())(const double& time,
193 const Vector<double>& x,
194 const Vector<double>& n,
196 {
197 return Traction_fct_pt;
198 }
199
200 /// Reference to the pressure function pointer
201 void (*&pressure_fct_pt())(const double& time,
202 const Vector<double>& x,
203 const Vector<double>& n,
204 double& pressure)
205 {
206 return Pressure_fct_pt;
207 }
208
209
210 /// Return the residuals
212 {
214 }
215
216
217 /// Fill in contribution from Jacobian
219 DenseMatrix<double>& jacobian)
220 {
221 // Call the residuals
223 }
224
225 /// Specify the value of nodal zeta from the face geometry
226 /// The "global" intrinsic coordinate of the element when
227 /// viewed as part of a geometric object should be given by
228 /// the FaceElement representation, by default (needed to break
229 /// indeterminacy if bulk element is SolidElement)
230 double zeta_nodal(const unsigned& n,
231 const unsigned& k,
232 const unsigned& i) const
233 {
234 return FaceElement::zeta_nodal(n, k, i);
235 }
236
237 /// Output function
238 void output(std::ostream& outfile)
239 {
240 FiniteElement::output(outfile);
241 }
242
243 /// Output function
244 void output(std::ostream& outfile, const unsigned& n_plot)
245 {
246 FiniteElement::output(outfile, n_plot);
247 }
248
249 /// C_style output function
250 void output(FILE* file_pt)
251 {
252 FiniteElement::output(file_pt);
253 }
254
255 /// C-style output function
256 void output(FILE* file_pt, const unsigned& n_plot)
257 {
258 FiniteElement::output(file_pt, n_plot);
259 }
260
261
262 /// Compute traction vector at specified local coordinate
263 /// Should only be used for post-processing; ignores dependence
264 /// on integration point!
265 void traction(const double& time,
266 const Vector<double>& s,
268
269 /// Compute pressure value at specified local coordinate
270 /// Should only be used for post-processing; ignores dependence
271 /// on integration point!
272 void pressure(const double& time,
273 const Vector<double>& s,
274 double& pressure);
275 };
276
277 /// ////////////////////////////////////////////////////////////////////
278 /// ////////////////////////////////////////////////////////////////////
279 /// ////////////////////////////////////////////////////////////////////
280
281 //=====================================================================
282 /// Compute traction vector at specified local coordinate
283 /// Should only be used for post-processing; ignores dependence
284 /// on integration point!
285 //=====================================================================
286 template<class ELEMENT>
288 const Vector<double>& s,
289 Vector<double>& traction)
290 {
291 unsigned n_dim = this->nodal_dimension();
292
293 // Position vector
294 Vector<double> x(n_dim);
295 interpolated_x(s, x);
296
297 // Outer unit normal
298 Vector<double> unit_normal(n_dim);
299 outer_unit_normal(s, unit_normal);
300
301 // Dummy
302 unsigned ipt = 0;
303
304 // Pressure value
305 get_traction(time, ipt, x, unit_normal, traction);
306 }
307 //=====================================================================
308 /// Compute pressure value at specified local coordinate
309 /// Should only be used for post-processing; ignores dependence
310 /// on integration point!
311 //=====================================================================
312 template<class ELEMENT>
314 const Vector<double>& s,
315 double& pressure)
316 {
317 unsigned n_dim = this->nodal_dimension();
318
319 // Position vector
320 Vector<double> x(n_dim);
321 interpolated_x(s, x);
322
323 // Outer unit normal
324 Vector<double> unit_normal(n_dim);
325 outer_unit_normal(s, unit_normal);
326
327 // Dummy
328 unsigned ipt = 0;
329
330 // Pressure value
331 get_pressure(time, ipt, x, unit_normal, pressure);
332 }
333
334
335 //=====================================================================
336 /// Return the residuals for the PoroelasticityFaceElement equations
337 //=====================================================================
338 template<class ELEMENT>
341 {
342 // Find out how many nodes there are
343 unsigned n_node = nnode();
344
345 // Get continuous time from timestepper of first node
346 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
347
348#ifdef PARANOID
349 // Find out how many positional dofs there are
350 unsigned n_position_type = this->nnodal_position_type();
351 if (n_position_type != 1)
352 {
353 throw OomphLibError("Poroelasticity equations are not yet implemented "
354 "for more than one position type",
355 OOMPH_CURRENT_FUNCTION,
356 OOMPH_EXCEPTION_LOCATION);
357 }
358#endif
359
360 // Find out the dimension of the node
361 unsigned n_dim = this->nodal_dimension();
362
363 unsigned n_q_basis = Element_pt->nq_basis();
364 unsigned n_q_basis_edge = Element_pt->nq_basis_edge();
365
366 // Integer to hold the local equation number
367 int local_eqn = 0;
368
369 // Set up memory for the shape functions
370 // Note that in this case, the number of lagrangian coordinates is always
371 // equal to the dimension of the nodes
372 Shape psi(n_node);
373 DShape dpsids(n_node, n_dim - 1);
374
375 Shape q_basis(n_q_basis, n_dim);
376
377 // Set the value of n_intpt
378 unsigned n_intpt = integral_pt()->nweight();
379
380 // Storage for the local coordinates
381 Vector<double> s_face(n_dim - 1, 0.0), s_bulk(n_dim, 0.0);
382
383 // mjr TODO enable if/when eta is added to Poroelasticity elements
384 /*double eta = Element_pt->eta();*/
385
386 // Loop over the integration points
387 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
388 {
389 // Get the integral weight
390 double w = integral_pt()->weight(ipt);
391
392 // Only need to call the local derivatives
393 dshape_local_at_knot(ipt, psi, dpsids);
394
395 // Assign values of s in FaceElement and local coordinates in bulk element
396 for (unsigned i = 0; i < n_dim - 1; i++)
397 {
398 s_face[i] = integral_pt()->knot(ipt, i);
399 }
400
401 s_bulk = local_coordinate_in_bulk(s_face);
402
403 // Get the q basis at bulk local coordinate s_bulk, corresponding to face
404 // local coordinate s_face
405 Element_pt->get_q_basis(s_bulk, q_basis);
406
407 // Calculate the Eulerian and Lagrangian coordinates
408 Vector<double> interpolated_x(n_dim, 0.0);
409
410 // Also calculate the surface Vectors (derivatives wrt local coordinates)
411 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
412
413 // Calculate displacements and derivatives
414 for (unsigned l = 0; l < n_node; l++)
415 {
416 // Loop over directions
417 for (unsigned i = 0; i < n_dim; i++)
418 {
419 // Calculate the Eulerian coords
420 const double x_local = nodal_position(l, i);
421 interpolated_x[i] += x_local * psi(l);
422
423 // Loop over LOCAL derivative directions, to calculate the tangent(s)
424 for (unsigned j = 0; j < n_dim - 1; j++)
425 {
426 interpolated_A(j, i) += x_local * dpsids(l, j);
427 }
428 }
429 }
430
431 // Now find the local metric tensor from the tangent vectors
432 DenseMatrix<double> A(n_dim - 1);
433 for (unsigned i = 0; i < n_dim - 1; i++)
434 {
435 for (unsigned j = 0; j < n_dim - 1; j++)
436 {
437 // Initialise surface metric tensor to zero
438 A(i, j) = 0.0;
439
440 // Take the dot product
441 for (unsigned k = 0; k < n_dim; k++)
442 {
443 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
444 }
445 }
446 }
447
448 // Get the outer unit normal
449 Vector<double> interpolated_normal(n_dim);
450 outer_unit_normal(ipt, interpolated_normal);
451
452 // Find the determinant of the metric tensor
453 double Adet = 0.0;
454 switch (n_dim)
455 {
456 case 2:
457 Adet = A(0, 0);
458 break;
459 case 3:
460 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
461 break;
462 default:
463 throw OomphLibError("Wrong dimension in PoroelasticityFaceElement",
464 OOMPH_CURRENT_FUNCTION,
465 OOMPH_EXCEPTION_LOCATION);
466 }
467
468 // Premultiply the weights and the square-root of the determinant of
469 // the metric tensor
470 double W = w * sqrt(Adet);
471
472 // Now calculate the traction load
473 Vector<double> traction(n_dim);
474 get_traction(time, ipt, interpolated_x, interpolated_normal, traction);
475
476 // Now calculate the load
477 double pressure;
478 get_pressure(time, ipt, interpolated_x, interpolated_normal, pressure);
479
480 // Loop over the test functions, nodes of the element
481 for (unsigned l = 0; l < n_node; l++)
482 {
483 // Loop over the displacement components
484 for (unsigned i = 0; i < n_dim; i++)
485 {
486 local_eqn = this->nodal_local_eqn(l, Element_pt->u_index(i));
487 /*IF it's not a boundary condition*/
488 if (local_eqn >= 0)
489 {
490 // Add the traction loading terms to the residuals
491 residuals[local_eqn] -= traction[i] * psi(l) * W;
492 } // End of if not boundary condition
493 }
494 } // End of loop over shape functions
495
496 // Loop over the q edge test functions only (internal basis functions
497 // have zero normal component on the boundary)
498 for (unsigned l = 0; l < n_q_basis_edge; l++)
499 {
500 local_eqn = this->nodal_local_eqn(1, Element_pt->q_edge_index(l));
501
502 /*IF it's not a boundary condition*/
503 if (local_eqn >= 0)
504 {
505 // Loop over the displacement components
506 for (unsigned i = 0; i < n_dim; i++)
507 {
508 // Add the loading terms to the residuals
509 residuals[local_eqn] +=
510 pressure * q_basis(l, i) * interpolated_normal[i] * W;
511 }
512 } // End of if not boundary condition
513 } // End of loop over shape functions
514 } // End of loop over integration points
515 }
516
517} // namespace oomph
518
519#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
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,...
Definition: elements.h:4497
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
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:3050
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:5132
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
A class for elements that allow the imposition of an applied pressure in the Darcy equations....
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, double &pressure) pressure_fct_pt()
Reference to the pressure function pointer.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
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...
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 output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
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...
void output(std::ostream &outfile)
Output function.
void(*&)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
ELEMENT * Element_pt
pointer to the bulk element this face element is attached to
void fill_in_contribution_to_residuals_darcy_face(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void output(FILE *file_pt)
C_style output function.
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...
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, const unsigned &n_plot)
Output function.
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...
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...
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
void Zero_traction_fct(const double &time, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
void Zero_pressure_fct(const double &time, const Vector< double > &x, const Vector< double > &N, double &load)
Default load function (zero pressure)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...