axisym_solid_traction_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 equations of elasticity
28
29#ifndef OOMPH_AXISYMM_SOLID_TRACTION_ELEMENTS_HEADER
30#define OOMPH_AXISYMM_SOLID_TRACTION_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37// OOMPH-LIB headers
38#include "../generic/Qelements.h"
39#include "../generic/hermite_elements.h"
40
41namespace oomph
42{
43 //======================================================================
44 /// A class for elements that allow the imposition of an applied traction
45 /// in the principle of virtual displacements.
46 /// The geometrical information can be read from the FaceGeometry<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 AxisymmetricSolidTractionElement : 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 const Vector<double>& n,
59 Vector<double>& result);
60
61 protected:
62 /// Return the surface traction force
63 void get_traction(const double& time,
64 const Vector<double>& x,
65 const Vector<double>& n,
66 Vector<double>& result) const
67 {
68 // If the function pointer is zero return zero
69 if (Traction_fct_pt == 0)
70 {
71 // Loop over dimensions and set body forces to zero
72 // It's axisymmetric, so we only have "two" dimensions
73 for (unsigned i = 0; i < 2; i++)
74 {
75 result[i] = 0.0;
76 }
77 }
78 // Otherwise call the function
79 else
80 {
81 (*Traction_fct_pt)(time, x, n, result);
82 }
83 }
84
85 public:
86 /// Constructor, which takes a "bulk" element and
87 /// the value of the index and its limit
89 const int& face_index)
90 : FaceGeometry<ELEMENT>(), FaceElement()
91 {
92 // Attach the geometrical information to the element. N.B. This function
93 // also assigns nbulk_value from the required_nvalue of the bulk element
94 element_pt->build_face_element(face_index, this);
95
96 // Set the body force function pointer to zero
98 }
99
100 /// Return the imposed traction pointer
101 void (*&traction_fct_pt())(const double&,
102 const Vector<double>&,
103 const Vector<double>&,
105 {
106 return Traction_fct_pt;
107 }
108
109 /// Return the residuals
111
112 /// Return the jacobian
114 DenseMatrix<double>& jacobian)
115 {
117 // Call the generic FD jacobian calculation
119 jacobian);
120 }
121
122 /// Overload the output function
123 void output(std::ostream& outfile)
124 {
125 FiniteElement::output(outfile);
126 }
127
128 /// Output function: x,y,[z],u,v,[w],p in tecplot format
129 void output(std::ostream& outfile, const unsigned& n_plot)
130 {
131 FiniteElement::output(outfile, n_plot);
132 }
133
134 /// Overload the output function
135 void output(FILE* file_pt)
136 {
137 FiniteElement::output(file_pt);
138 }
139
140 /// Output function: x,y,[z],u,v,[w],p in tecplot format
141 void output(FILE* file_pt, const unsigned& n_plot)
142 {
143 FiniteElement::output(file_pt, n_plot);
144 }
145 };
146
147
148 /// //////////////////////////////////////////////////////////////////////
149 /// //////////////////////////////////////////////////////////////////////
150 /// //////////////////////////////////////////////////////////////////////
151
152
153 //=======================================================================
154 /// Return the residuals for the AxisymmetricSolidTractionElements
155 //=======================================================================
156 template<class ELEMENT>
157 void AxisymmetricSolidTractionElement<
158 ELEMENT>::fill_in_contribution_to_residuals(Vector<double>& residuals)
159 {
160 // Find out how many nodes there are
161 unsigned n_node = nnode();
162 // Find out how many positional dofs there are
163 unsigned n_position_type = nnodal_position_type();
164
165 // Integer to hold the local equation number
166 int local_eqn = 0;
167
168 // Set up memory for the shape functions
169 // The surface is 1D, so we only have one local derivative
170 Shape psi(n_node, n_position_type);
171 DShape dpsids(n_node, n_position_type, 1);
172
173 // Set the value of n_intpt
174 unsigned n_intpt = integral_pt()->nweight();
175
176 // Loop over the integration points
177 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
178 {
179 // Get the integral weight
180 double w = integral_pt()->weight(ipt);
181
182 // Only need to call the local derivatives
183 dshape_local_at_knot(ipt, psi, dpsids);
184
185 // Calculate the global position and lagrangian coordinate
186 Vector<double> interpolated_x(2, 0.0), interpolated_xi(2, 0.0);
187 // Calculate the global and lagrangian derivtives wrt the local
188 // coordinates
189 Vector<double> interpolated_dxds(2, 0.0), interpolated_dxids(2, 0.0);
190
191 // Calculate displacements and derivatives
192 for (unsigned l = 0; l < n_node; l++)
193 {
194 // Loop over positional dofs
195 for (unsigned k = 0; k < n_position_type; k++)
196 {
197 // Loop over the number of lagrangian coordinates (2)
198 for (unsigned i = 0; i < 2; i++)
199 {
200 // Calculate the global position
201 interpolated_x[i] +=
202 nodal_position_gen(l, bulk_position_type(k), i) * psi(l, k);
203 interpolated_xi[i] +=
204 this->lagrangian_position_gen(l, bulk_position_type(k), i) *
205 psi(l, k);
206 // Calculate the derivatives of the global and lagrangian
207 // coordinates
208 interpolated_dxds[i] +=
209 nodal_position_gen(l, bulk_position_type(k), i) * dpsids(l, k, 0);
210 interpolated_dxids[i] +=
211 this->lagrangian_position_gen(l, bulk_position_type(k), i) *
212 dpsids(l, k, 0);
213 }
214 }
215 }
216
217 // Now calculate the entries of the deformed surface metric tensor
218 // Now find the local deformed metric tensor from the tangent Vectors
220 // The off-diagonal terms are Zero
221 A(0, 1) = A(1, 0) = 0.0;
222 // The diagonal terms are a little complicated
223 A(0, 0) =
224 (interpolated_dxds[0] - interpolated_x[1] * interpolated_dxids[1]) *
225 (interpolated_dxds[0] - interpolated_x[1] * interpolated_dxids[1]) +
226 (interpolated_dxds[1] + interpolated_x[0] * interpolated_dxids[1]) *
227 (interpolated_dxds[1] + interpolated_x[0] * interpolated_dxids[1]);
228
229
230 A(1, 1) = (interpolated_x[0] * sin(interpolated_xi[1]) +
231 interpolated_x[1] * cos(interpolated_xi[1])) *
232 (interpolated_x[0] * sin(interpolated_xi[1]) +
233 interpolated_x[1] * cos(interpolated_xi[1]));
234
235 // Premultiply the weights and the square-root of the determinant of
236 // the metric tensor
237 double W = w * sqrt(A(0, 0) * A(1, 1));
238
239 // Also find the normal -- just the cross product of the metric tensors
240 // but I want to express it in terms of e_r and e_theta components
241 // N.B. There is an issue at theta = 0,pi, where the normal is e_{r},
242 // but given that I never assemble it, should be OK!
243 // The minus sign is chosen to ensure that the normal is really outward
244 Vector<double> interpolated_normal(2);
245 // Component in the e_{r} direction
246 interpolated_normal[0] =
247 -1.0 *
248 (interpolated_x[0] * sin(interpolated_xi[1]) +
249 interpolated_x[1] * cos(interpolated_xi[1])) *
250 (interpolated_dxds[1] + interpolated_x[0] * interpolated_dxids[1]);
251 // Component in the e_{theta} direction
252 interpolated_normal[1] =
253 -1.0 *
254 (interpolated_x[0] * sin(interpolated_xi[1]) +
255 interpolated_x[1] * cos(interpolated_xi[1])) *
256 (interpolated_x[1] * interpolated_dxids[1] - interpolated_dxds[0]);
257
258 // If we're on the north or south face need to flip normal
259 if (s_fixed_value() == -1)
260 {
261 interpolated_normal[0] *= -1.0;
262 interpolated_normal[1] *= -1.0;
263 }
264
265 // Now adjust and scale the normal
266 double length = 0.0;
267 for (unsigned i = 0; i < 2; i++)
268 {
269 interpolated_normal[i] *= normal_sign();
270 length += interpolated_normal[i] * interpolated_normal[i];
271 }
272 for (unsigned i = 0; i < 2; i++)
273 {
274 interpolated_normal[i] /= sqrt(length);
275 }
276
277 // Now calculate the load
278 Vector<double> traction(2);
279
280 // Normal is outwards
281 get_traction(time(), interpolated_x, interpolated_normal, traction);
282
283 //=====LOAD TERMS FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
284
285 // Loop over the test functions, nodes of the element
286 for (unsigned l = 0; l < n_node; l++)
287 {
288 // Loop of types of dofs
289 for (unsigned k = 0; k < n_position_type; k++)
290 {
291 // Loop over the displacement components
292 for (unsigned i = 0; i < 2; i++)
293 {
294 local_eqn = this->position_local_eqn(l, bulk_position_type(k), i);
295 /*IF it's not a boundary condition*/
296 if (local_eqn >= 0)
297 {
298 // Add the loading terms to the residuals
299 residuals[local_eqn] -= traction[i] * psi(l, k) * W;
300 }
301 }
302 } // End of if not boundary condition
303 } // End of loop over shape functions
304 } // End of loop over integration points
305 }
306
307
308} // namespace oomph
309
310#endif
cstr elem_len * i
Definition: cfortran.h:603
A class for elements that allow the imposition of an applied traction in the principle of virtual dis...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
void get_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result) const
Return the surface traction force.
void output(std::ostream &outfile)
Overload the output function.
void output(FILE *file_pt)
Overload the output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Return the jacobian.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void(*&)(const double &, const Vector< double > &, const Vector< double > &, Vector< double > &) traction_fct_pt()
Return the imposed traction pointer.
void output(FILE *file_pt, const unsigned &n_plot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
AxisymmetricSolidTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function.
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...