axisym_displ_based_fvk_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 axisymmetric FoepplvonKarman elements
27#ifndef OOMPH_AXISYM_DISPL_BASED_FOEPPLVONKARMAN_ELEMENTS_HEADER
28#define OOMPH_AXISYM_DISPL_BASED_FOEPPLVONKARMAN_ELEMENTS_HEADER
29
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36// OOMPH-LIB headers
37#include "../generic/nodes.h"
38#include "../generic/Qelements.h"
39#include "../generic/oomph_utilities.h"
40
41namespace oomph
42{
43 //=============================================================
44 /// A class for all isoparametric elements that solve the
45 /// axisYm Foeppl von Karman equations in a displacement based formulation.
46 ///
47 /// This contains the generic maths. Shape functions, geometric
48 /// mapping etc. must get implemented in derived class.
49 //=============================================================
51 {
52 public:
53 /// Function pointer to pressure function fct(r,f(r)) --
54 /// r is a Vector!
55 typedef void (*AxisymFoepplvonKarmanPressureFctPt)(const double& r,
56 double& f);
57
58 /// Constructor (must initialise the Pressure_fct_pt). Also
59 /// set physical parameters to their default values.
61
62 /// Broken copy constructor
64 const AxisymFoepplvonKarmanEquations& dummy) = delete;
65
66 /// Broken assignment operator
68
69 /// Poisson's ratio
70 const double& nu() const
71 {
72#ifdef PARANOID
73 if (Nu_pt == 0)
74 {
75 std::stringstream error_stream;
76 error_stream << "Nu has not yet been set!" << std::endl;
77 throw OomphLibError(
78 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
79 }
80#endif
81 return *Nu_pt;
82 }
83
84 /// Pointer to Poisson's ratio
85 double*& nu_pt()
86 {
87 return Nu_pt;
88 }
89
90 /// FvK parameter
91 const double& eta() const
92 {
93 return *Eta_pt;
94 }
95
96 /// Pointer to FvK parameter
97 double*& eta_pt()
98 {
99 return Eta_pt;
100 }
101
102 /// Return the index at which the i-th unknown value
103 /// is stored. The default value, i, is appropriate for single-physics
104 /// problems. By default, these are:
105 /// 0: transverse displacement w
106 /// 1: laplacian w
107 /// 2: radial displacement u
108 /// In derived multi-physics elements, this function should be overloaded
109 /// to reflect the chosen storage scheme. Note that these equations require
110 /// that the unknown is always stored at the same index at each node.
111 virtual inline unsigned nodal_index_fvk(const unsigned& i = 0) const
112 {
113 return i;
114 }
115
116 /// Output with default number of plot points
117 void output(std::ostream& outfile)
118 {
119 const unsigned n_plot = 5;
120 output(outfile, n_plot);
121 }
122
123 /// Output FE representation of soln: r,w,sigma_r_r,sigma_phi_phi
124 /// at n_plot plot points
125 void output(std::ostream& outfile, const unsigned& n_plot);
126
127 /// C_style output with default number of plot points
128 void output(FILE* file_pt)
129 {
130 const unsigned n_plot = 5;
131 output(file_pt, n_plot);
132 }
133
134 /// C-style output FE representation of soln: r,w at
135 /// n_plot plot points
136 void output(FILE* file_pt, const unsigned& n_plot);
137
138 /// Output exact soln: r,w_exact at n_plot plot points
139 void output_fct(std::ostream& outfile,
140 const unsigned& n_plot,
142
143 /// Output exact soln: r,w_exact at
144 /// n_plot plot points (dummy time-dependent version to
145 /// keep intel compiler happy)
146 virtual void output_fct(
147 std::ostream& outfile,
148 const unsigned& n_plot,
149 const double& time,
151 {
152 throw OomphLibError(
153 "There is no time-dependent output_fct() for Foeppl von Karman"
154 "elements ",
155 OOMPH_CURRENT_FUNCTION,
156 OOMPH_EXCEPTION_LOCATION);
157 }
158
159 /// Get error against and norm of exact solution
160 void compute_error(std::ostream& outfile,
162 double& error,
163 double& norm);
164
165
166 /// Dummy, time dependent error checker
167 void compute_error(std::ostream& outfile,
169 const double& time,
170 double& error,
171 double& norm)
172 {
173 throw OomphLibError(
174 "There is no time-dependent compute_error() for Foeppl von Karman"
175 "elements",
176 OOMPH_CURRENT_FUNCTION,
177 OOMPH_EXCEPTION_LOCATION);
178 }
179
180 /// Access function: Pointer to pressure function
182 {
183 return Pressure_fct_pt;
184 }
185
186 /// Access function: Pointer to pressure function. Const version
188 {
189 return Pressure_fct_pt;
190 }
191
192 /// Get pressure term at (Eulerian) position r. This function is
193 /// virtual to allow overloading in multi-physics problems where
194 /// the strength of the pressure function might be determined by
195 /// another system of equations.
196 inline virtual void get_pressure_fvk(const unsigned& ipt,
197 const double& r,
198 double& pressure) const
199 {
200 // If no pressure function has been set, return zero
201 if (Pressure_fct_pt == 0)
202 {
203 pressure = 0.0;
204 }
205 else
206 {
207 // Get pressure strength
208 (*Pressure_fct_pt)(r, pressure);
209 }
210 }
211
212 /// Get gradient of deflection: gradient[i] = dw/dr_i */
214 Vector<double>& gradient) const
215 {
216 // Find out how many nodes there are in the element
217 const unsigned n_node = nnode();
218
219 // Get the index at which the unknown is stored
220 const unsigned w_nodal_index = nodal_index_fvk(0);
221
222 // Set up memory for the shape and test functions
223 Shape psi(n_node);
224 DShape dpsidr(n_node, 1);
225
226 // Call the derivatives of the shape and test functions
227 dshape_eulerian(s, psi, dpsidr);
228
229 // Initialise to zero
230 gradient[0] = 0.0;
231
232 // Loop over nodes
233 for (unsigned l = 0; l < n_node; l++)
234 {
235 gradient[0] += this->nodal_value(l, w_nodal_index) * dpsidr(l, 0);
236 }
237 }
238
239 /// Fill in the residuals with this element's contribution
241
242
243 // hierher Jacobian not yet implemented
244 // void fill_in_contribution_to_jacobian(Vector<double> &residuals,
245 // DenseMatrix<double> &jacobian);
246
247
248 /// Return FE representation of transverse displacement
249 inline double interpolated_w_fvk(const Vector<double>& s) const
250 {
251 // Find number of nodes
252 const unsigned n_node = nnode();
253
254 // Get the index at which the transverse displacement unknown is stored
255 const unsigned w_nodal_index = nodal_index_fvk(0);
256
257 // Local shape function
258 Shape psi(n_node);
259
260 // Find values of shape function
261 shape(s, psi);
262
263 // Initialise value of u
264 double interpolated_w = 0.0;
265
266 // Loop over the local nodes and sum
267 for (unsigned l = 0; l < n_node; l++)
268 {
269 interpolated_w += this->nodal_value(l, w_nodal_index) * psi[l];
270 }
271
272 return (interpolated_w);
273 }
274
275
276 /// Return FE representation of radial displacement
277 inline double interpolated_u_fvk(const Vector<double>& s) const
278 {
279 // Find number of nodes
280 const unsigned n_node = nnode();
281
282 // Get the index at which the radial displacement unknown is stored
283 const unsigned u_nodal_index = nodal_index_fvk(2);
284
285 // Local shape function
286 Shape psi(n_node);
287
288 // Find values of shape function
289 shape(s, psi);
290
291 // Initialise value of u
292 double interpolated_u = 0.0;
293
294 // Loop over the local nodes and sum
295 for (unsigned l = 0; l < n_node; l++)
296 {
297 interpolated_u += this->nodal_value(l, u_nodal_index) * psi[l];
298 }
299
300 return (interpolated_u);
301 }
302
303
304 /// Compute in-plane stresses. Return boolean to indicate success
305 /// (false if attempt to evaluate stresses at zero radius)
307 double& sigma_r_r,
308 double& sigma_phi_phi) const;
309
310
311 /// Self-test: Return 0 for OK
312 unsigned self_test();
313
314
315 // switch back on and test!
316
317
318 /// Sets a flag to signify that we are solving the linear,
319 /// pure bending equations, and pin all the nodal values that will
320 /// not be used in this case
322 {
323 // Set the boolean flag
325
326 // Get the index of the first FvK nodal value
327 unsigned first_fvk_nodal_index = nodal_index_fvk();
328
329 // Get the total number of FvK nodal values (assuming they are stored
330 // contiguously) at node 0 (it's the same at all nodes anyway)
331 unsigned total_fvk_nodal_indices = 3;
332
333 // Get the number of nodes in this element
334 unsigned n_node = nnode();
335
336 // Loop over the appropriate nodal indices
337 for (unsigned index = first_fvk_nodal_index + 2;
338 index < first_fvk_nodal_index + total_fvk_nodal_indices;
339 index++)
340 {
341 // Loop over the nodes in the element
342 for (unsigned inod = 0; inod < n_node; inod++)
343 {
344 // Pin the nodal value at the current index
345 node_pt(inod)->pin(index);
346 }
347 }
348 }
349
350
351 protected:
352 /// Shape/test functions and derivs w.r.t. to global coords at
353 /// local coord. s; return Jacobian of mapping
355 const Vector<double>& s,
356 Shape& psi,
357 DShape& dpsidr,
358 Shape& test,
359 DShape& dtestdr) const = 0;
360
361
362 /// Shape/test functions and derivs w.r.t. to global coords at
363 /// integration point ipt; return Jacobian of mapping
365 const unsigned& ipt,
366 Shape& psi,
367 DShape& dpsidr,
368 Shape& test,
369 DShape& dtestdr) const = 0;
370
371 /// Pointer to FvK parameter
372 double* Eta_pt;
373
374 /// Pointer to pressure function:
376
377 /// Pointer to Poisson's ratio
378 double* Nu_pt;
379
380 /// Flag which stores whether we are using a linear,
381 /// pure bending model instead of the full non-linear Foeppl-von Karman
383
384 private:
385 /// Default value for physical constants
387 };
388
389
390 /// ////////////////////////////////////////////////////////////////////////
391 /// ////////////////////////////////////////////////////////////////////////
392 /// ////////////////////////////////////////////////////////////////////////
393
394
395 //======================================================================
396 /// Axisym FoepplvonKarmanElement elements are 1D
397 /// Foeppl von Karman elements with isoparametric interpolation for the
398 /// function.
399 //======================================================================
400 template<unsigned NNODE_1D>
402 : public virtual QElement<1, NNODE_1D>,
403 public virtual AxisymFoepplvonKarmanEquations
404 {
405 private:
406 /// Static int that holds the number of variables at
407 /// nodes: always the same
408 static const unsigned Initial_Nvalue;
409
410 public:
411 /// Constructor: Call constructors for QElement and
412 /// AxisymFoepplvonKarmanEquations
414 : QElement<1, NNODE_1D>(), AxisymFoepplvonKarmanEquations()
415 {
416 }
417
418 /// Broken copy constructor
420 const AxisymFoepplvonKarmanElement<NNODE_1D>& dummy) = delete;
421
422 /// Broken assignment operator
424
425 /// Required # of `values' (pinned or dofs)
426 /// at node n
427 inline unsigned required_nvalue(const unsigned& n) const
428 {
429 return Initial_Nvalue;
430 }
431
432
433 /// Output function:
434 /// r,w,sigma_r_r,sigma_phi_phi
435 void output(std::ostream& outfile)
436 {
438 }
439
440 /// Output function:
441 /// r,w,sigma_r_r,sigma_phi_phi at n_plot plot points
442 void output(std::ostream& outfile, const unsigned& n_plot)
443 {
445 }
446
447 /// C-style output function:
448 /// r,w
449 void output(FILE* file_pt)
450 {
452 }
453
454 /// C-style output function:
455 /// r,w at n_plot plot points
456 void output(FILE* file_pt, const unsigned& n_plot)
457 {
459 }
460
461 /// Output function for an exact solution:
462 /// r,w_exact at n_plot plot points
463 void output_fct(std::ostream& outfile,
464 const unsigned& n_plot,
466 {
468 outfile, n_plot, exact_soln_pt);
469 }
470
471 /// Output function for a time-dependent exact solution.
472 /// r,w_exact at n_plot plot points
473 /// (Calls the steady version)
474 void output_fct(std::ostream& outfile,
475 const unsigned& n_plot,
476 const double& time,
478 {
480 outfile, n_plot, time, exact_soln_pt);
481 }
482
483
484 protected:
485 /// Shape, test functions & derivs. w.r.t. to global coords.
486 /// Return Jacobian.
488 Shape& psi,
489 DShape& dpsidr,
490 Shape& test,
491 DShape& dtestdr) const;
492
493 /// Shape, test functions & derivs. w.r.t. to global coords. at
494 /// integration point ipt. Return Jacobian.
496 const unsigned& ipt,
497 Shape& psi,
498 DShape& dpsidr,
499 Shape& test,
500 DShape& dtestdr) const;
501 };
502
503
504 // Inline functions:
505
506 //======================================================================
507 /// Define the shape functions and test functions and derivatives
508 /// w.r.t. global coordinates and return Jacobian of mapping.
509 ///
510 /// Galerkin: Test functions = shape functions
511 //======================================================================
512 template<unsigned NNODE_1D>
513 double AxisymFoepplvonKarmanElement<
514 NNODE_1D>::dshape_and_dtest_eulerian_axisym_fvk(const Vector<double>& s,
515 Shape& psi,
516 DShape& dpsidr,
517 Shape& test,
518 DShape& dtestdr) const
519
520 {
521 // Call the geometrical shape functions and derivatives
522 const double J = this->dshape_eulerian(s, psi, dpsidr);
523
524 // Set the test functions equal to the shape functions
525 test = psi;
526 dtestdr = dpsidr;
527
528 // Return the jacobian
529 return J;
530 }
531
532
533 //======================================================================
534 /// Define the shape functions and test functions and derivatives
535 /// w.r.t. global coordinates and return Jacobian of mapping.
536 ///
537 /// Galerkin: Test functions = shape functions
538 //======================================================================
539 template<unsigned NNODE_1D>
542 Shape& psi,
543 DShape& dpsidr,
544 Shape& test,
545 DShape& dtestdr) const
546 {
547 // Call the geometrical shape functions and derivatives
548 const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidr);
549
550 // Set the pointers of the test functions
551 test = psi;
552 dtestdr = dpsidr;
553
554 // Return the jacobian
555 return J;
556 }
557
558} // namespace oomph
559
560#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,w at n_plot plot points.
double dshape_and_dtest_eulerian_axisym_fvk(const Vector< double > &s, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const
Shape, test functions & derivs. w.r.t. to global coords. Return Jacobian.
void output(FILE *file_pt)
C-style output function: r,w.
double dshape_and_dtest_eulerian_at_knot_axisym_fvk(const unsigned &ipt, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const
Shape, test functions & derivs. w.r.t. to global coords. at integration point ipt....
void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output function for a time-dependent exact solution. r,w_exact at n_plot plot points (Calls the stead...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,w,sigma_r_r,sigma_phi_phi at n_plot plot points.
AxisymFoepplvonKarmanElement(const AxisymFoepplvonKarmanElement< NNODE_1D > &dummy)=delete
Broken copy constructor.
void output(std::ostream &outfile)
Output function: r,w,sigma_r_r,sigma_phi_phi.
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
unsigned required_nvalue(const unsigned &n) const
Required # of ‘values’ (pinned or dofs) at node n.
void operator=(const AxisymFoepplvonKarmanElement< NNODE_1D > &)=delete
Broken assignment operator.
AxisymFoepplvonKarmanElement()
Constructor: Call constructors for QElement and AxisymFoepplvonKarmanEquations.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output function for an exact solution: r,w_exact at n_plot plot points.
A class for all isoparametric elements that solve the axisYm Foeppl von Karman equations in a displac...
const double & eta() const
FvK parameter.
AxisymFoepplvonKarmanPressureFctPt pressure_fct_pt() const
Access function: Pointer to pressure function. Const version.
AxisymFoepplvonKarmanPressureFctPt Pressure_fct_pt
Pointer to pressure function:
static double Default_Physical_Constant_Value
Default value for physical constants.
void operator=(const AxisymFoepplvonKarmanEquations &)=delete
Broken assignment operator.
double *& nu_pt()
Pointer to Poisson's ratio.
virtual void output_fct(std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,w_exact at n_plot plot points (dummy time-dependent version to keep intel compil...
void output(std::ostream &outfile)
Output with default number of plot points.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals with this element's contribution.
bool interpolated_stress(const Vector< double > &s, double &sigma_r_r, double &sigma_phi_phi) const
Compute in-plane stresses. Return boolean to indicate success (false if attempt to evaluate stresses ...
void output(FILE *file_pt)
C_style output with default number of plot points.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
unsigned self_test()
Self-test: Return 0 for OK.
double *& eta_pt()
Pointer to FvK parameter.
virtual double dshape_and_dtest_eulerian_at_knot_axisym_fvk(const unsigned &ipt, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
virtual unsigned nodal_index_fvk(const unsigned &i=0) const
Return the index at which the i-th unknown value is stored. The default value, i, is appropriate for ...
void(* AxisymFoepplvonKarmanPressureFctPt)(const double &r, double &f)
Function pointer to pressure function fct(r,f(r)) – r is a Vector!
void get_gradient_of_deflection(const Vector< double > &s, Vector< double > &gradient) const
Get gradient of deflection: gradient[i] = dw/dr_i *‍/.
AxisymFoepplvonKarmanPressureFctPt & pressure_fct_pt()
Access function: Pointer to pressure function.
AxisymFoepplvonKarmanEquations()
Constructor (must initialise the Pressure_fct_pt). Also set physical parameters to their default valu...
virtual double dshape_and_dtest_eulerian_axisym_fvk(const Vector< double > &s, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const =0
Shape/test functions and derivs w.r.t. to global coords at local coord. s; return Jacobian of mapping...
bool Linear_bending_model
Flag which stores whether we are using a linear, pure bending model instead of the full non-linear Fo...
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,w_exact at n_plot plot points.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void use_linear_bending_model()
Sets a flag to signify that we are solving the linear, pure bending equations, and pin all the nodal ...
double interpolated_u_fvk(const Vector< double > &s) const
Return FE representation of radial displacement.
double interpolated_w_fvk(const Vector< double > &s) const
Return FE representation of transverse displacement.
virtual void get_pressure_fvk(const unsigned &ipt, const double &r, double &pressure) const
Get pressure term at (Eulerian) position r. This function is virtual to allow overloading in multi-ph...
AxisymFoepplvonKarmanEquations(const AxisymFoepplvonKarmanEquations &dummy)=delete
Broken copy constructor.
const double & nu() const
Poisson's ratio.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
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 shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
Definition: elements.h:1765
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...