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