polar_stress_integral_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 integrate the shear stress
27// along either side wall
28
29#ifndef OOMPH_POLAR_STRESS_INTEGRAL_ELEMENTS_HEADER
30#define OOMPH_POLAR_STRESS_INTEGRAL_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <config.h>
35#endif
36
37
38// OOMPH-LIB headers
39#include "../generic/Qelements.h"
40
41namespace oomph
42{
43 //======================================================================
44 /// A class for elements that allow the imposition of an applied traction
45 /// to the Navier--Stokes equations
46 /// The geometrical information can be read from the FaceGeometery<ELEMENT>
47 /// class and and thus, we can be generic enough without the need to have
48 /// a separate equations class
49 //======================================================================
50 template<class ELEMENT>
51 class PolarStressIntegralElement : public virtual FaceGeometry<ELEMENT>,
52 public virtual FaceElement
53 {
54 private:
55 /// Pointer to an imposed traction function
56 void (*Traction_fct_pt)(const double& time,
57 const Vector<double>& x,
58 Vector<double>& result);
59
60 /// The highest dimension of the problem
61 unsigned Dim;
62
63 protected:
64 /// Access function that returns the local equation numbers
65 /// for velocity components.
66 /// u_local_eqn(n,i) = local equation number or < 0 if pinned.
67 /// The default is to asssume that n is the local node number
68 /// and the i-th velocity component is the i-th unknown stored at the node.
69 virtual inline int u_local_eqn(const unsigned& n, const unsigned& i)
70 {
71 return nodal_local_eqn(n, i);
72 }
73
74 /// Function to compute the shape and test functions and to return
75 /// the Jacobian of mapping
76 inline double shape_and_test_at_knot(const unsigned& ipt,
77 Shape& psi,
78 Shape& test) const
79 {
80 // Find number of nodes
81 unsigned n_node = nnode();
82 // Calculate the shape functions
83 shape_at_knot(ipt, psi);
84 // Set the test functions to be the same as the shape functions
85 for (unsigned i = 0; i < n_node; i++)
86 {
87 test[i] = psi[i];
88 }
89 // Return the value of the jacobian
90 return J_eulerian_at_knot(ipt);
91 }
92
93 /// Pointer to the angle alpha
94 double* Alpha_pt;
95
96 // Traction elements need to know whether they're at the inlet or outlet
97 // as the unit outward normal has a differing sign dependent on
98 // the boundary
99 // phi=-1, phi=1
101
102 public:
103 /// Alpha
104 const double& alpha() const
105 {
106 return *Alpha_pt;
107 }
108
109 /// Pointer to Alpha
110 double*& alpha_pt()
111 {
112 return Alpha_pt;
113 }
114
115 /// Boundary
116 const int boundary() const
117 {
118 return Boundary;
119 }
120
121 /// Function to set boundary
122 void set_boundary(int bound)
123 {
124 Boundary = bound;
125 }
126
127 /// Function to calculate the shear stress along boundary
128 double get_shear_stress();
129
130 /// Constructor, which takes a "bulk" element and the value of the index
131 /// and its limit
133 const int& face_index)
134 : FaceGeometry<ELEMENT>(), FaceElement()
135 {
136 // Attach the geometrical information to the element. N.B. This function
137 // also assigns nbulk_value from the required_nvalue of the bulk element
138 element_pt->build_face_element(face_index, this);
139
140#ifdef PARANOID
141 {
142 // Check that the element is not a refineable 3d element
143 ELEMENT* elem_pt = dynamic_cast<ELEMENT*>(element_pt);
144
145 // If it's three-d
146 if (elem_pt->dim() == 3)
147 {
148 // Is it refineable
149 RefineableElement* ref_el_pt =
150 dynamic_cast<RefineableElement*>(elem_pt);
151 if (ref_el_pt != 0)
152 {
153 if (this->has_hanging_nodes())
154 {
155 throw OomphLibError("This flux element will not work correctly "
156 "if nodes are hanging\n",
157 OOMPH_CURRENT_FUNCTION,
158 OOMPH_EXCEPTION_LOCATION);
159 }
160 }
161 }
162 }
163#endif
164
165 // Set the dimension from the dimension of the first node
166 Dim = this->node_pt(0)->ndim();
167 }
168
169 /// Destructor should not delete anything
171
172 /// This function returns just the residuals
174 {
175 // Do nothing
176 }
177
178 /// This function returns the residuals and the jacobian
180 DenseMatrix<double>& jacobian)
181 {
182 // Do nothing
183 }
184
185 /// Compute the element's residual Vector and the jacobian matrix
186 /// Plus the mass matrix especially for eigenvalue problems
188 Vector<double>& residuals,
189 DenseMatrix<double>& jacobian,
190 DenseMatrix<double>& mass_matrix)
191 {
192 // Do nothing
193 }
194
195 /// Overload the output function
196 void output(std::ostream& outfile)
197 {
198 FiniteElement::output(outfile);
199 }
200
201 /// Output function: x,y,[z],u,v,[w],p in tecplot format
202 void output(std::ostream& outfile, const unsigned& nplot)
203 {
204 FiniteElement::output(outfile, nplot);
205 }
206
207 /// local velocities
208 double u(const unsigned& l, const unsigned& i)
209 {
210 return nodal_value(l, i);
211 }
212
213 /// local position
214 double x(const unsigned& l, const unsigned& i)
215 {
216 return nodal_position(l, i);
217 }
218 };
219
220 //============================================================================
221 /// Function that returns the shear stress
222 //============================================================================
223 template<class ELEMENT>
225 {
226 // Storage for shear stress
227 double dudphi, shear_contribution = 0.0;
228
229 // Set the value of n_intpt
230 unsigned n_intpt = integral_pt()->nweight();
231
232 // Storage for local coordinate
233 Vector<double> s_local(1);
234 // Storage for local coordinate in bulk
235 Vector<double> s_bulk(2);
236
237 // Find out how many nodes there are
238 unsigned n_node = nnode();
239
240 // Set up memory for the shape and test functions
241 Shape psif(n_node), testf(n_node);
242
243 // Loop over the integration points
244 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
245 {
246 // Get the integral weight
247 double w = integral_pt()->weight(ipt);
248
249 // Find the shape and test functions and return the Jacobian
250 // of the mapping
251 double J = shape_and_test_at_knot(ipt, psif, testf);
252
253 // Premultiply the weights and the Jacobian
254 double W = w * J;
255
256 // Need to find position to feed into Traction function
257 Vector<double> interpolated_x(Dim);
258
259 // Initialise to zero
260 for (unsigned i = 0; i < Dim; i++)
261 {
262 interpolated_x[i] = 0.0;
263 }
264
265 // Calculate velocities and derivatives
266 for (unsigned l = 0; l < n_node; l++)
267 {
268 // Loop over velocity components
269 for (unsigned i = 0; i < Dim; i++)
270 {
271 interpolated_x[i] += nodal_position(l, i) * psif[l];
272 }
273 }
274
275 // Get the local coordinate
276 s_local[0] = integral_pt()->knot(ipt, 0);
277
278 // Get bulk coordinate
279 s_bulk = this->local_coordinate_in_bulk(s_local);
280
281 // Upcast from GeneralisedElement to the present element
282 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(this->Bulk_element_pt);
283
284 // Get du_dphi from bulk element
285 dudphi = el_pt->interpolated_dudx_pnst(s_bulk, 0, 1);
286
287 // The contribution to the unweighted shear stress
288 shear_contribution += dudphi * W;
289 // dudphi*interpolated_x[0]*W;
290
291 } // End of loop over integration points
292
293 return shear_contribution;
294
295 } // End of get_shear_stress
296
297} // End of namespace oomph
298
299#endif
cstr elem_len * i
Definition: cfortran.h:603
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 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:5328
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
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
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:1432
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
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:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
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:3220
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
An OomphLibError object which should be thrown when an run-time error is encountered....
A class for elements that allow the imposition of an applied traction to the Navier–Stokes equations ...
void output(std::ostream &outfile)
Overload the output function.
double * Alpha_pt
Pointer to the angle alpha.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
double u(const unsigned &l, const unsigned &i)
local velocities
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Compute the element's residual Vector and the jacobian matrix Plus the mass matrix especially for eig...
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.
~PolarStressIntegralElement()
Destructor should not delete anything.
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to an imposed traction function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the jacobian.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
void set_boundary(int bound)
Function to set boundary.
double x(const unsigned &l, const unsigned &i)
local position
virtual int u_local_eqn(const unsigned &n, const unsigned &i)
Access function that returns the local equation numbers for velocity components. u_local_eqn(n,...
unsigned Dim
The highest dimension of the problem.
double get_shear_stress()
Function to calculate the shear stress along boundary.
PolarStressIntegralElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
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
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:60
//////////////////////////////////////////////////////////////////// ////////////////////////////////...