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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for axisymmetric FoepplvonKarman elements
31 #ifndef OOMPH_AXISYM_DISPL_BASED_FOEPPLVONKARMAN_ELEMENTS_HEADER
32 #define OOMPH_AXISYM_DISPL_BASED_FOEPPLVONKARMAN_ELEMENTS_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 // OOMPH-LIB headers
41 #include "../generic/nodes.h"
42 #include "../generic/Qelements.h"
43 #include "../generic/oomph_utilities.h"
44 
45 namespace oomph
46 {
47 
48 //=============================================================
49 /// A class for all isoparametric elements that solve the
50 /// axisYm Foeppl von Karman equations in a displacement based formulation.
51 ///
52 /// This contains the generic maths. Shape functions, geometric
53 /// mapping etc. must get implemented in derived class.
54 //=============================================================
56  {
57 
58  public:
59 
60  /// \short Function pointer to pressure function fct(r,f(r)) --
61  /// r is a Vector!
62  typedef void (*AxisymFoepplvonKarmanPressureFctPt)(const double& r,
63  double& f);
64 
65  /// \short Constructor (must initialise the Pressure_fct_pt). Also
66  /// set physical parameters to their default values.
68  {}
69 
70  /// Broken copy constructor
72  {
73  BrokenCopy::broken_copy("AxisymFoepplvonKarmanEquations");
74  }
75 
76  /// Broken assignment operator
78  {
79  BrokenCopy::broken_assign("AxisymFoepplvonKarmanEquations");
80  }
81 
82  /// Poisson's ratio
83  const double &nu() const
84  {
85 #ifdef PARANOID
86  if (Nu_pt==0)
87  {
88  std::stringstream error_stream;
89  error_stream
90  << "Nu has not yet been set!" << std::endl;
91  throw OomphLibError(
92  error_stream.str(),
93  OOMPH_CURRENT_FUNCTION,
94  OOMPH_EXCEPTION_LOCATION);
95  }
96 #endif
97  return *Nu_pt;
98  }
99 
100  /// Pointer to Poisson's ratio
101  double* &nu_pt() {return Nu_pt;}
102 
103  /// FvK parameter
104  const double &eta() const {return *Eta_pt;}
105 
106  /// Pointer to FvK parameter
107  double* &eta_pt() {return Eta_pt;}
108 
109  /// \short Return the index at which the i-th unknown value
110  /// is stored. The default value, i, is appropriate for single-physics
111  /// problems. By default, these are:
112  /// 0: transverse displacement w
113  /// 1: laplacian w
114  /// 2: radial displacement u
115  /// In derived multi-physics elements, this function should be overloaded
116  /// to reflect the chosen storage scheme. Note that these equations require
117  /// that the unknown is always stored at the same index at each node.
118  virtual inline unsigned nodal_index_fvk(const unsigned& i=0)
119  const {return i;}
120 
121  /// Output with default number of plot points
122  void output(std::ostream &outfile)
123  {
124  const unsigned n_plot=5;
125  output(outfile,n_plot);
126  }
127 
128  /// \short Output FE representation of soln: r,w,sigma_r_r,sigma_phi_phi
129  /// at n_plot plot points
130  void output(std::ostream &outfile, const unsigned &n_plot);
131 
132  /// C_style output with default number of plot points
133  void output(FILE* file_pt)
134  {
135  const unsigned n_plot=5;
136  output(file_pt,n_plot);
137  }
138 
139  /// \short C-style output FE representation of soln: r,w at
140  /// n_plot plot points
141  void output(FILE* file_pt, const unsigned &n_plot);
142 
143  /// Output exact soln: r,w_exact at n_plot plot points
144  void output_fct(std::ostream &outfile, const unsigned &n_plot,
146 
147  /// \short Output exact soln: r,w_exact at
148  /// n_plot plot points (dummy time-dependent version to
149  /// keep intel compiler happy)
150  virtual void output_fct(std::ostream &outfile, const unsigned &n_plot,
151  const double& time,
153  exact_soln_pt)
154  {
155  throw OomphLibError(
156  "There is no time-dependent output_fct() for Foeppl von Karman"
157  "elements ",
158  OOMPH_CURRENT_FUNCTION,
159  OOMPH_EXCEPTION_LOCATION);
160  }
161 
162  /// Get error against and norm of exact solution
163  void compute_error(std::ostream &outfile,
165  double& error, double& norm);
166 
167 
168  /// Dummy, time dependent error checker
169  void compute_error(std::ostream &outfile,
171  const double& time, double& error, 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  {return Pressure_fct_pt;}
183 
184  /// Access function: Pointer to pressure function. Const version
186  {return Pressure_fct_pt;}
187 
188  /// \short Get pressure term at (Eulerian) position r. This function is
189  /// virtual to allow overloading in multi-physics problems where
190  /// the strength of the pressure function might be determined by
191  /// another system of equations.
192  inline virtual void get_pressure_fvk(const unsigned& ipt,
193  const double& r,
194  double& pressure) const
195  {
196  //If no pressure function has been set, return zero
197  if(Pressure_fct_pt==0)
198  {
199  pressure = 0.0;
200  }
201  else
202  {
203  // Get pressure strength
204  (*Pressure_fct_pt)(r,pressure);
205  }
206  }
207 
208  /// Get gradient of deflection: gradient[i] = dw/dr_i */
210  Vector<double>& gradient) const
211  {
212  //Find out how many nodes there are in the element
213  const unsigned n_node = nnode();
214 
215  //Get the index at which the unknown is stored
216  const unsigned w_nodal_index = nodal_index_fvk(0);
217 
218  //Set up memory for the shape and test functions
219  Shape psi(n_node);
220  DShape dpsidr(n_node,1);
221 
222  //Call the derivatives of the shape and test functions
223  dshape_eulerian(s,psi,dpsidr);
224 
225  //Initialise to zero
226  gradient[0] = 0.0;
227 
228  // Loop over nodes
229  for(unsigned l=0;l<n_node;l++)
230  {
231  gradient[0] += this->nodal_value(l,w_nodal_index)*dpsidr(l,0);
232  }
233  }
234 
235  /// Fill in the residuals with this element's contribution
237 
238 
239  // hierher Jacobian not yet implemented
240  //void fill_in_contribution_to_jacobian(Vector<double> &residuals,
241  // DenseMatrix<double> &jacobian);
242 
243 
244  /// \short Return FE representation of transverse displacement
245  inline double interpolated_w_fvk(const Vector<double> &s) const
246  {
247  //Find number of nodes
248  const unsigned n_node = nnode();
249 
250  //Get the index at which the transverse displacement unknown is stored
251  const unsigned w_nodal_index = nodal_index_fvk(0);
252 
253  //Local shape function
254  Shape psi(n_node);
255 
256  //Find values of shape function
257  shape(s,psi);
258 
259  //Initialise value of u
260  double interpolated_w = 0.0;
261 
262  //Loop over the local nodes and sum
263  for(unsigned l=0;l<n_node;l++)
264  {
265  interpolated_w += this->nodal_value(l,w_nodal_index)*psi[l];
266  }
267 
268  return(interpolated_w);
269  }
270 
271 
272  /// \short Return FE representation of radial displacement
273  inline double interpolated_u_fvk(const Vector<double> &s) const
274  {
275  //Find number of nodes
276  const unsigned n_node = nnode();
277 
278  //Get the index at which the radial displacement unknown is stored
279  const unsigned u_nodal_index = nodal_index_fvk(2);
280 
281  //Local shape function
282  Shape psi(n_node);
283 
284  //Find values of shape function
285  shape(s,psi);
286 
287  //Initialise value of u
288  double interpolated_u = 0.0;
289 
290  //Loop over the local nodes and sum
291  for(unsigned l=0;l<n_node;l++)
292  {
293  interpolated_u += this->nodal_value(l,u_nodal_index)*psi[l];
294  }
295 
296  return(interpolated_u);
297  }
298 
299 
300  /// \short Compute in-plane stresses. Return boolean to indicate success
301  /// (false if attempt to evaluate stresses at zero radius)
302  bool interpolated_stress(const Vector<double> &s, double& sigma_r_r,
303  double& sigma_phi_phi) const;
304 
305 
306  /// \short Self-test: Return 0 for OK
307  unsigned self_test();
308 
309 
310 // switch back on and test!
311 
312 
313  /// \short Sets a flag to signify that we are solving the linear,
314  /// pure bending equations, and pin all the nodal values that will
315  /// not be used in this case
317  {
318  // Set the boolean flag
319  Linear_bending_model = true;
320 
321  // Get the index of the first FvK nodal value
322  unsigned first_fvk_nodal_index = nodal_index_fvk();
323 
324  // Get the total number of FvK nodal values (assuming they are stored
325  // contiguously) at node 0 (it's the same at all nodes anyway)
326  unsigned total_fvk_nodal_indices = 3;
327 
328  // Get the number of nodes in this element
329  unsigned n_node = nnode();
330 
331  // Loop over the appropriate nodal indices
332  for(unsigned index=first_fvk_nodal_index+2;
333  index<first_fvk_nodal_index+total_fvk_nodal_indices;
334  index++)
335  {
336  // Loop over the nodes in the element
337  for(unsigned inod=0;inod<n_node;inod++)
338  {
339  // Pin the nodal value at the current index
340  node_pt(inod)->pin(index);
341  }
342  }
343  }
344 
345 
346  protected:
347 
348  /// \short Shape/test functions and derivs w.r.t. to global coords at
349  /// local coord. s; return Jacobian of mapping
351  (const Vector<double> &s,
352  Shape &psi,
353  DShape &dpsidr,
354  Shape &test,
355  DShape &dtestdr)
356  const=0;
357 
358 
359  /// \short Shape/test functions and derivs w.r.t. to global coords at
360  /// integration point ipt; return Jacobian of mapping
362  (const unsigned &ipt,
363  Shape &psi,
364  DShape &dpsidr,
365  Shape &test,
366  DShape &dtestdr)
367  const=0;
368 
369  /// Pointer to FvK parameter
370  double *Eta_pt;
371 
372  /// Pointer to pressure function:
374 
375  /// Pointer to Poisson's ratio
376  double* Nu_pt;
377 
378  private:
379 
380  /// Default value for physical constants
382 
383  /// \short Flag which stores whether we are using a linear,
384  /// pure bending model instead of the full non-linear Foeppl-von Karman
386 
387  };
388 
389 
390 
391 ///////////////////////////////////////////////////////////////////////////
392 ///////////////////////////////////////////////////////////////////////////
393 ///////////////////////////////////////////////////////////////////////////
394 
395 
396 //======================================================================
397 /// Axisym FoepplvonKarmanElement elements are 1D
398 /// Foeppl von Karman elements with isoparametric interpolation for the
399 /// function.
400 //======================================================================
401  template <unsigned NNODE_1D>
402  class AxisymFoepplvonKarmanElement : public virtual QElement<1,NNODE_1D>,
403  public virtual AxisymFoepplvonKarmanEquations
404  {
405 
406  private:
407 
408  /// \short Static int that holds the number of variables at
409  /// nodes: always the same
410  static const unsigned Initial_Nvalue;
411 
412  public:
413 
414  ///\short Constructor: Call constructors for QElement and
415  /// AxisymFoepplvonKarmanEquations
418  {}
419 
420  /// Broken copy constructor
423  {
424  BrokenCopy::broken_copy("AxisymFoepplvonKarmanElement");
425  }
426 
427  /// Broken assignment operator
429  {
430  BrokenCopy::broken_assign("AxisymFoepplvonKarmanElement");
431  }
432 
433  /// \short Required # of `values' (pinned or dofs)
434  /// at node n
435  inline unsigned required_nvalue(const unsigned &n) const
436  {return Initial_Nvalue;}
437 
438 
439  /// \short Output function:
440  /// r,w,sigma_r_r,sigma_phi_phi
441  void output(std::ostream &outfile)
443 
444  /// \short Output function:
445  /// r,w,sigma_r_r,sigma_phi_phi at n_plot plot points
446  void output(std::ostream &outfile, const unsigned &n_plot)
447  {AxisymFoepplvonKarmanEquations::output(outfile,n_plot);}
448 
449  /// \short C-style output function:
450  /// r,w
451  void output(FILE* file_pt)
453 
454  /// \short C-style output function:
455  /// r,w at n_plot plot points
456  void output(FILE* file_pt, const unsigned &n_plot)
457  {AxisymFoepplvonKarmanEquations::output(file_pt,n_plot);}
458 
459  /// \short Output function for an exact solution:
460  /// r,w_exact at n_plot plot points
461  void output_fct(std::ostream &outfile, const unsigned &n_plot,
464  output_fct(outfile,n_plot,exact_soln_pt);}
465 
466  /// \short Output function for a time-dependent exact solution.
467  /// r,w_exact at n_plot plot points
468  /// (Calls the steady version)
469  void output_fct(std::ostream &outfile, const unsigned &n_plot,
470  const double& time,
473  (outfile,n_plot,time,exact_soln_pt);}
474 
475 
476  protected:
477 
478  /// \short Shape, test functions & derivs. w.r.t. to global coords.
479  /// Return Jacobian.
481  const Vector<double> &s, Shape &psi, DShape &dpsidr,
482  Shape &test, DShape &dtestdr) const;
483 
484  /// \short Shape, test functions & derivs. w.r.t. to global coords. at
485  /// integration point ipt. Return Jacobian.
487  (const unsigned& ipt,
488  Shape &psi,
489  DShape &dpsidr,
490  Shape &test,
491  DShape &dtestdr)
492  const;
493 
494  };
495 
496 
497 
498 //Inline functions:
499 
500 //======================================================================
501 /// Define the shape functions and test functions and derivatives
502 /// w.r.t. global coordinates and return Jacobian of mapping.
503 ///
504 /// Galerkin: Test functions = shape functions
505 //======================================================================
506  template<unsigned NNODE_1D>
509  const Vector<double> &s,
510  Shape &psi,
511  DShape &dpsidr,
512  Shape &test,
513  DShape &dtestdr) const
514 
515  {
516  //Call the geometrical shape functions and derivatives
517  const double J = this->dshape_eulerian(s,psi,dpsidr);
518 
519  //Set the test functions equal to the shape functions
520  test = psi;
521  dtestdr= dpsidr;
522 
523  //Return the jacobian
524  return J;
525  }
526 
527 
528 //======================================================================
529 /// Define the shape functions and test functions and derivatives
530 /// w.r.t. global coordinates and return Jacobian of mapping.
531 ///
532 /// Galerkin: Test functions = shape functions
533 //======================================================================
534  template<unsigned NNODE_1D>
537  const unsigned &ipt,
538  Shape &psi,
539  DShape &dpsidr,
540  Shape &test,
541  DShape &dtestdr) const
542  {
543  //Call the geometrical shape functions and derivatives
544  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidr);
545 
546  //Set the pointers of the test functions
547  test = psi;
548  dtestdr = dpsidr;
549 
550  //Return the jacobian
551  return J;
552  }
553 
554 }
555 
556 #endif
557 
void output(FILE *file_pt)
C_style output with default number of plot points.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
AxisymFoepplvonKarmanEquations()
Constructor (must initialise the Pressure_fct_pt). Also set physical parameters to their default valu...
double *& nu_pt()
Pointer to Poisson's ratio.
AxisymFoepplvonKarmanElement()
Constructor: Call constructors for QElement and AxisymFoepplvonKarmanEquations.
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:3225
cstr elem_len * i
Definition: cfortran.h:607
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:2458
A general Finite Element class.
Definition: elements.h:1271
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.
AxisymFoepplvonKarmanPressureFctPt Pressure_fct_pt
Pointer to pressure function:
const double & eta() const
FvK parameter.
void operator=(const AxisymFoepplvonKarmanElement< NNODE_1D > &)
Broken assignment operator.
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 operator=(const AxisymFoepplvonKarmanEquations &)
Broken assignment operator.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
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. Return Jacobian.
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
const double & nu() const
Poisson's ratio.
void output(std::ostream &outfile)
Output with default number of plot points.
void output(FILE *file_pt)
C-style output function: r,w.
unsigned self_test()
Self-test: Return 0 for OK.
void(* AxisymFoepplvonKarmanPressureFctPt)(const double &r, double &f)
Function pointer to pressure function fct(r,f(r)) – r is a Vector!
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...
static char t char * s
Definition: cfortran.h:572
unsigned required_nvalue(const unsigned &n) const
Required # of `values' (pinned or dofs) at node n.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function: r,w at n_plot plot points.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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 compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Dummy, time dependent error checker.
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 ...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double interpolated_u_fvk(const Vector< double > &s) const
Return FE representation of radial displacement.
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.
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 ...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: r,w,sigma_r_r,sigma_phi_phi at n_plot plot points.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
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...
static double Default_Physical_Constant_Value
Default value for physical constants.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
double *& eta_pt()
Pointer to FvK parameter.
bool Linear_bending_model
Flag which stores whether we are using a linear, pure bending model instead of the full non-linear Fo...
double interpolated_w_fvk(const Vector< double > &s) const
Return FE representation of transverse displacement.
void get_gradient_of_deflection(const Vector< double > &s, Vector< double > &gradient) const
Get gradient of deflection: gradient[i] = dw/dr_i */.
AxisymFoepplvonKarmanEquations(const AxisymFoepplvonKarmanEquations &dummy)
Broken copy constructor.
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 ...
AxisymFoepplvonKarmanPressureFctPt pressure_fct_pt() const
Access function: Pointer to pressure function. Const version.
void use_linear_bending_model()
Sets a flag to signify that we are solving the linear, pure bending equations, and pin all the nodal ...
void output(std::ostream &outfile)
Output function: r,w,sigma_r_r,sigma_phi_phi.
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...
double * Nu_pt
Pointer to Poisson's ratio.
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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals with this element's contribution.
AxisymFoepplvonKarmanPressureFctPt & pressure_fct_pt()
Access function: Pointer to pressure function.