shell_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 KL shell elements
31 #ifndef OOMPH_KIRCHHOFF_LOVE_SHELL_ELEMENTS_HEADER
32 #define OOMPH_KIRCHHOFF_LOVE_SHELL_ELEMENTS_HEADER
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36  #include <oomph-lib-config.h>
37 #endif
38 
39 //OOMPH-LIB header
40 #include "../generic/hermite_elements.h"
41 #include "../generic/geom_objects.h"
42 #include "../generic/fsi.h"
43 #include "../generic/stored_shape_function_elements.h"
44 #include "../generic/matrices.h"
45 
46 namespace oomph
47 {
48 
49 //========================================================================
50 /// A class for elements that solves the equations of Kirchhoff Love shell
51 /// thin-shell theory
52 //========================================================================
54 {
55  private:
56 
57  /// Static default value for the Poisson's ratio
58  static double Default_nu_value;
59 
60  /// Static default value for the timescale ratio (1.0 for natural scaling)
61  static double Default_lambda_sq_value;
62 
63  /// Static default value for the thickness ratio
64  static double Default_h_value;
65 
66  /// Boolean flag to ignore membrane terms
68 
69  /// Pointer to Poisson's ratio
70  double *Nu_pt;
71 
72  /// Pointer to dimensionless wall thickness
73  double *H_pt;
74 
75  /// Pointer to timescale ratio (non-dimensional density)
76  double *Lambda_sq_pt;
77 
78  /// \short Pointer to membrane pre-stress terms -- should
79  /// probably generalise this to function pointers at some point
81 
82  /// Static default for prestress (set to zero)
83  static double Zero_prestress;
84 
85 protected:
86 
87  /// Invert a DIM by DIM matrix
88  inline double calculate_contravariant
89  (double A[2][2], double Aup[2][2]);
90 
91  /// Default load function (zero traction)
92  static void Zero_traction_fct(const Vector<double> &xi,
93  const Vector<double> &x,
94  const Vector<double> &N,
95  Vector<double> &load);
96 
97  /// \short Pointer to load vector function: Its arguments are:
98  /// Lagrangian coordinate, Eulerian coordinate, normal vector and
99  /// load vector itself (not all of the input arguments will be
100  /// required for all specific load functions but the list should
101  /// cover all cases)
103  const Vector<double> &x,
104  const Vector<double> &N,
105  Vector<double> &load);
106 
107 
108  /// \short Pointer to the GeomObject that specifies the undeformed midplane
109  /// of the shell
111 
112 
113  /// short Helper function to Return the residuals for
114  /// the equations of KL shell theory. This is used to prevent
115  /// a call of fill_in_contribution_to_residuals in
116  /// the function fill_in_contribution_to_jacobian that can
117  /// lead to virtual inheritance woes if this element is ever
118  /// used as part of a multi-physics element.
120 
121  /// \short Get the load vector for the computation of the rate of work
122  /// done by the load. Here we simply forward this to
123  /// load_vector(...) but allow it to be overloaded in derived classes
124  /// to allow the computation of the rate of work due to constituent
125  /// parts of the load vector (e.g. the fluid traction in an FSI
126  /// problem). Pass number of integration point (dummy),
127  /// Lagr. coordinate and normal vector and return the load vector
128  /// (not all of the input arguments will be
129  /// required for all specific load functions but the list should
130  /// cover all cases).
132  const unsigned &intpt,
133  const Vector<double>& xi,
134  const Vector<double>& x,
135  const Vector<double>& N,
136  Vector<double>& load)
137  {
138  load_vector(intpt, xi, x, N ,load);
139  }
140 
141 
142 public:
143 
144  /// Constructor: Initialise physical parameter values to defaults
146  {
147  //Set the default values of physical parameters
151 
152  // Don't ignore membrane terms
153  Ignore_membrane_terms=false;
154 
155  //Default load is zero traction
157 
158  // Default prestress is zero
159  Prestress_pt.resize(2,2);
164  }
165 
166  /// Access to the load vector function pointer
167  void (* &load_vector_fct_pt())(const Vector<double>& xi,
168  const Vector<double>& x,
169  const Vector<double>& N,
170  Vector<double>& load)
171  {return Load_vector_fct_pt;}
172 
173  /// \short Get the load vector: Pass number of integration point (dummy),
174  /// Lagr. coordinate and normal vector and return the load vector
175  /// (not all of the input arguments will be
176  /// required for all specific load functions but the list should
177  /// cover all cases). This function is virtual so it can be
178  /// overloaded for FSI.
179  virtual void load_vector(const unsigned &intpt,
180  const Vector<double>& xi,
181  const Vector<double>& x,
182  const Vector<double>& N,
183  Vector<double>& load)
184  {
185  Load_vector_fct_pt(xi,x,N,load);
186  }
187 
188 
189  /// \short Set pointer to (i,j)-th component of second Piola Kirchhoff
190  /// membrane prestress to specified value (automatically imposes
191  /// symmetry for off-diagonal entries)
192  void set_prestress_pt(const unsigned& i,
193  const unsigned& j,
194  double* value_pt)
195  {
196  Prestress_pt(i,j)=value_pt;
197  Prestress_pt(j,i)=value_pt;
198  }
199 
200  /// \short Return (i,j)-th component of second Piola Kirchhoff membrane
201  /// prestress
202  double prestress(const unsigned& i,
203  const unsigned& j)
204  {
205  return *Prestress_pt(i,j);
206  }
207 
208  /// Set to disable the calculation of membrane terms
210 
211  /// Set to renable the calculation of membrane terms (default)
213 
214  /// Return the Poisson's ratio
215  const double &nu() const {return *Nu_pt;}
216 
217  /// Return the wall thickness to undeformed radius ratio
218  const double &h() const {return *H_pt;}
219 
220  /// Return the timescale ratio (non-dimensional density)
221  const double& lambda_sq() const {return *Lambda_sq_pt;}
222 
223  /// Return a pointer to the Poisson ratio
224  double* &nu_pt() {return Nu_pt;}
225 
226  /// Return a pointer to the non-dim wall thickness
227  double* &h_pt() {return H_pt;}
228 
229  /// Return a pointer to timescale ratio (nondim density)
230  double* &lambda_sq_pt() {return Lambda_sq_pt;}
231 
232  /// \short Return a reference to geometric object that specifies the shell's
233  /// undeformed geometry
235 
236  /// \short Get normal vector
237  void get_normal(const Vector<double>& s, Vector<double>& N);
238 
239  ///Overload the standard fill in residuals contribution
241  {
242  //Simply call the shell residuals
244  }
245 
246  /// Return the jacobian is calculated by finite differences by default,
248  DenseMatrix<double> &jacobian);
249 
250  /// \short Get potential (strain) and kinetic energy of the element
251  void get_energy(double& pot_en, double& kin_en);
252 
253 
254  /// \short Get strain and bending tensors; returns pair comprising the
255  /// determinant of the undeformed (*.first) and deformed (*.second)
256  /// midplane metric tensor.
257  std::pair<double,double> get_strain_and_bend(const Vector<double>& s,
258  DenseDoubleMatrix& strain,
259  DenseDoubleMatrix& bend);
260 
261 
262  /// \short Get integral of instantaneous rate of work done on
263  /// the wall due to the load returned by the virtual
264  /// function load_vector_for_rate_of_work_computation(...).
265  /// In the current class
266  /// the latter function simply de-references the external
267  /// load but this can be overloaded in derived classes
268  /// (e.g. in FSI elements) to determine the rate of work done
269  /// by individual constituents of this load (e.g. the fluid load
270  /// in an FSI problem).
271  double load_rate_of_work();
272 
273  /// Generic FiniteElement output function
274  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
275 
276  /// Generic FiniteElement output function
277  void output(std::ostream &outfile, const unsigned &n_plot)
278  {FiniteElement::output(outfile,n_plot);}
279 
280  /// Generic FiniteElement output function
281  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
282 
283  /// Generic FiniteElement output function
284  void output(FILE* file_pt, const unsigned &n_plot)
285  {FiniteElement::output(file_pt,n_plot);}
286 
287 };
288 
289 
290 
291 //==================================================================
292 /// Matrix inversion for 2 dimensions
293 //==================================================================
295 (double A[2][2], double Aup[2][2])
296 {
297  //Calculate determinant
298  double det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
299 
300  //Calculate entries of the inverse
301  Aup[0][0] = A[1][1]/det;
302  Aup[0][1] = -A[0][1]/det;
303  Aup[1][0] = -A[1][0]/det;
304  Aup[1][1] = A[0][0]/det;
305 
306  //Return determinant
307  return(det);
308 }
309 
310 
311 
312 
313 
314 //=======================================================================
315 /// An element that solves the Kirchhoff-Love shell theory equations
316 /// using Hermite interpolation (displacements
317 /// and slopes are interpolated separately. The local and global
318 /// (Lagrangian) coordinates are not assumed to be aligned.
319 /// N.B. It will be DOG SLOW.
320 //=======================================================================
321 class HermiteShellElement : public virtual SolidQHermiteElement<2>,
323 {
324 public:
325 
326  /// Constructor, there are no internal data points
329  {
330  //Set the number of dimensions at each node (3D nodes on 2D surface)
332  }
333 
334  /// Output position veloc and accel.
335  void output_with_time_dep_quantities(std::ostream &outfile,
336  const unsigned &n_plot);
337 
338  /// Overload the output function
339  void output(std::ostream &outfile) {SolidQHermiteElement<2>::output(outfile);}
340 
341  /// Output function
342  void output(std::ostream &outfile, const unsigned &n_plot);
343 
344  /// Overload the output function
345  void output(FILE* file_pt) {SolidQHermiteElement<2>::output(file_pt);}
346 
347  /// Output function
348  void output(FILE* file_pt, const unsigned &n_plot);
349 
350 };
351 
352 
353 
354 //=========================================================================
355 /// An element that solves the Kirchhoff-Love shell theory equations
356 /// using Hermite interpolation (displacements
357 /// and slopes are interpolated separately. The local and global
358 /// (Lagrangian) coordinates are assumed to be aligned so that the
359 /// Jacobian of the mapping between these coordinates is diagonal.
360 /// This significantly simplifies (and speeds up) the computation of the
361 /// derivatives of the shape functions.
362 //=========================================================================
364  public SolidDiagQHermiteElement<2>
365 {
366 
367 public:
368 
369  /// Constructor, there are no internal data points
372 
373 };
374 
375 
376 
377 ////////////////////////////////////////////////////////////////////////
378 ////////////////////////////////////////////////////////////////////////
379 ////////////////////////////////////////////////////////////////////////
380 
381 
382 
383 //=======================================================================
384 /// Face geometry for the HermiteShell elements: 1D SolidQHermiteElement
385 //=======================================================================
386 template<>
388  public virtual SolidQHermiteElement<1>
389 {
390 
391 public:
392 
393  /// \short Constructor [this was only required explicitly
394  /// from gcc 4.5.2 onwards...]
396 
397 };
398 
399 
400 
401 
402 ////////////////////////////////////////////////////////////////////////
403 ////////////////////////////////////////////////////////////////////////
404 ////////////////////////////////////////////////////////////////////////
405 
406 
407 //=========================================================================
408 /// Diag Hermite Kirchhoff Love shell "upgraded" to a FSIWallElement
409 /// (and thus, by inheritance, a GeomObject), so it can be used in FSI.
410 //=========================================================================
412  public virtual FSIWallElement
413 {
414  private:
415 
416  /// \short Get the load vector for the computation of the rate of work
417  /// done by the load. Can switch between full load and fluid load only.
418  /// Overloads the version in the shell element base class.
419  /// Pass number of integration point (dummy)
420  /// Lagr. coordinate and normal vector and return the load vector
421  /// (not all of the input arguments will be
422  /// required for all specific load functions but the list should
423  /// cover all cases).
425  const unsigned &intpt,
426  const Vector<double>& xi,
427  const Vector<double>& x,
428  const Vector<double>& N,
429  Vector<double>& load)
430  {
431  /// Get fluid-only load vector
433  {
434  fluid_load_vector(intpt,N,load);
435  }
436  // Get full load vector as per default
437  else
438  {
439  load_vector(intpt,xi,x,N,load);
440  }
441  }
442 
443  ///Boolean flag to indicate whether the normal is directed into the fluid
445 
446  /// \short Boolean flag to indicate if rate-of-work by load is to be
447  ///based on the fluid traction only
449 
450  public:
451 
452  /// \short Constructor: Create shell element as FSIWallElement (and thus,
453  /// by inheritance, a GeomObject) with two Lagrangian coordinates
454  /// and 3 Eulerian coordinates. By default, we assume that the
455  /// normal vector computed by KirchhoffLoveShellEquations::get_normal(...)
456  /// points into the fluid.
457  /// If this is not the case, call the access function
458  /// FSIDiagHermiteShellElement::set_normal_pointing_out_of_fluid()
462  {
463  unsigned n_lagr=2;
464  unsigned n_dim=3;
465  setup_fsi_wall_element(n_lagr,n_dim);
466  }
467 
468  /// \short Destructor: empty
470 
471  /// \short Set the normal computed by
472  /// KirchhoffLoveShellEquations::get_normal(...) to point into the fluid
474 
475  /// \short Set the normal computed by
476  /// KirchhoffLoveShellEquations::get_normal(...) to point out of the fluid
478 
479 
480 /// \short Derivative of position vector w.r.t. the SolidFiniteElement's
481  /// Lagrangian coordinates; evaluated at current time.
483  const Vector<double>& s, DenseMatrix<double> &drdxi) const;
484 
485 
486  /// \short Get integral of instantaneous rate of work done on
487  /// the wall due to the fluid load returned by the function
488  /// fluid_load_vector(...).
490  {
492  double tmp=load_rate_of_work();
494  return tmp;
495  }
496 
497 
498  /// \short Get the load vector: Pass number of the integration point,
499  /// Lagr. coordinate, Eulerian coordinate and normal vector
500  /// and return the load vector. (Not all of the input arguments will be
501  /// required for all specific load functions but the list should
502  /// cover all cases). We first evaluate the load function defined via
503  /// KirchhoffLoveShellEquations::load_vector_fct_pt() -- this
504  /// represents the non-FSI load on the shell, e.g. an external
505  /// pressure load. Then we add to this the FSI load due to
506  /// the traction exerted by the adjacent FSIFluidElements, taking
507  /// the sign of the normal into account.
508  void load_vector(const unsigned& intpt,
509  const Vector<double>& xi,
510  const Vector<double>& x,
511  const Vector<double>& N,
512  Vector<double>& load)
513  {
514  //Initially call the standard Load_vector_fct_pt
515  Load_vector_fct_pt(xi,x,N,load);
516 
517  //Memory for the FSI load
518  Vector<double> fsi_load(3);
519 
520  //Get the fluid load on the wall stress scale
521  fluid_load_vector(intpt,N,fsi_load);
522 
523  //If the normal is outer to the fluid switch the direction
524  double sign = 1.0;
525  if (!Normal_points_into_fluid) {sign = -1.0;}
526 
527  //Add the FSI load to the load vector
528  for(unsigned i=0;i<3;i++)
529  {
530  load[i] += sign*fsi_load[i];
531  }
532  }
533 
534  /// \short Get the Jacobian and residuals. Wrapper to generic FSI version;
535  /// that catches the case when we replace the Jacobian by the
536  /// mass matrix (for the consistent assignment of initial conditions).
538  DenseMatrix<double> &jacobian)
539  {
540  //Call the basic shell jacobian
542  fill_in_contribution_to_jacobian(residuals,jacobian);
543  //Fill in the external interaction by finite differences
545  }
546 
547 
548  /// \short The number of "DOF types" that degrees of freedom in this element
549  /// are sub-divided into: Just the solid degrees of freedom themselves.
550  unsigned ndof_types() const
551  {
552  return 1;
553  }
554 
555  /// \short Create a list of pairs for all unknowns in this element,
556  /// so that the first entry in each pair contains the global equation
557  /// number of the unknown, while the second one contains the number
558  /// of the "DOF types" that this unknown is associated with.
559  /// (Function can obviously only be called if the equation numbering
560  /// scheme has been set up.)
562  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const;
563 
564 };
565 
566 
567 
568 ////////////////////////////////////////////////////////////////////////
569 ////////////////////////////////////////////////////////////////////////
570 ////////////////////////////////////////////////////////////////////////
571 
572 
573 //======================================================================
574 /// Element that allows the imposition of boundary
575 /// conditions for a shell that is clamped to a 2D plane
576 /// that is specified by its normal. Constraint is applied by
577 /// a Lagrange multiplier.
578 /// \n
579 /// \b Note \b 1: Note that the introduction of the Lagrange
580 /// multiplier adds two additional values (relative to the number
581 /// of values before the addition of the FaceElement) to the nodes.
582 /// This ensures that nodes that are shared by adjacent FaceElements
583 /// are not resized repeatedly but also means that this won't work
584 /// if two "edges" of the shell (that share a node) are subject
585 /// to different constraints, each applied with its own
586 /// independent Lagrange multiplier. In such cases a modified
587 /// version of this class must be written.
588 /// \n
589 /// \b Note \b 2: The FaceGeometry for a HermiteShellElement is
590 /// the 1D two-node element SolidQHermiteElement<1> which has
591 /// four shape functions (two nodes, two types -- representing the
592 /// shape functions that interpolate the value and the derivative).
593 /// These are the "correct" shape functions for the interpolation
594 /// of the Lagrange multiplier and the isoparametric representation
595 /// of the geometry. However, when applying the contribution from the
596 /// constraint equation to the bulk equations, we have to take
597 /// all four types of dof into account so the element has to
598 /// reset the number of positional dofs to four. To avoid any clashes
599 /// we overload (the relevant subset of) the access functions to the
600 /// shape functions and their derivatives and set the shape functions
601 /// associated with the spurious positional dofs to zero. This is a bit
602 /// hacky but the only way (?) this can be done...
603 //======================================================================
605  public virtual FaceGeometry<HermiteShellElement>,
606  public virtual SolidFaceElement
607 {
608 
609 public:
610 
611  /// \short Constructor, takes the pointer to the "bulk" element and the
612  /// face index.
614  FiniteElement* const &bulk_el_pt, const int &face_index);
615 
616 
617  ///\short Broken empty constructor
619  {
620  throw OomphLibError(
621  "Don't call empty constructor for ClampedHermiteShellBoundaryConditionElement",
622  OOMPH_CURRENT_FUNCTION,
623  OOMPH_EXCEPTION_LOCATION);
624  }
625 
626  /// Broken copy constructor
629  {
631  "ClampedHermiteShellBoundaryConditionElement");
632  }
633 
634  /// Broken assignment operator
636  {
638  "ClampedHermiteShellBoundaryConditionElement");
639  }
640 
641  /// \short Set normal vector to clamping plane
642  void set_symmetry_line(const Vector<double>& normal_to_clamping_plane)
643  {
644  Normal_to_clamping_plane[0]=normal_to_clamping_plane[0];
645  Normal_to_clamping_plane[1]=normal_to_clamping_plane[1];
646  Normal_to_clamping_plane[2]=normal_to_clamping_plane[2];
647  }
648 
649  /// Fill in the element's contribution to its residual vector
651 
652 
653  //////////////////////////////////////////////////////////////////
654  // Note: We should also overload all other versions of shape(...)
655  // and dshape(...) but these are the only ones that are used.
656  //////////////////////////////////////////////////////////////////
657 
658 
659  /// \short Calculate the geometric shape functions
660  /// at local coordinate s. Set any "superfluous" shape functions to zero.
661  void shape(const Vector<double> &s, Shape &psi) const
662  {
663  // Initialise all of them to zero
664  unsigned n=psi.nindex1();
665  unsigned m=psi.nindex2();
666  for (unsigned i=0;i<n;i++)
667  {
668  for (unsigned j=0;j<m;j++)
669  {
670  psi(i,j)=0.0;
671  }
672  }
674  }
675 
676 
677  /// \short Calculate the geometric shape functions
678  /// at local coordinate s. Set any "superfluous" shape functions to zero.
679  void dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids) const
680  {
681  // Initialise all of them to zero
682  unsigned n=psi.nindex1();
683  unsigned m=psi.nindex2();
684  for (unsigned i=0;i<n;i++)
685  {
686  for (unsigned j=0;j<m;j++)
687  {
688  psi(i,j)=0.0;
689  }
690  }
691  unsigned n1=dpsids.nindex1();
692  unsigned n2=dpsids.nindex2();
693  unsigned n3=dpsids.nindex3();
694  for (unsigned i=0;i<n1;i++)
695  {
696  for (unsigned j=0;j<n2;j++)
697  {
698  for (unsigned k=0;k<n3;k++)
699  {
700  dpsids(i,j,k)=0.0;
701  }
702  }
703  }
705  }
706 
707  /// Output function -- forward to broken version in FiniteElement
708  /// until somebody decides what exactly they want to plot here...
709  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
710 
711  /// \short Output function
712  void output(std::ostream &outfile, const unsigned &n_plot);
713 
714  /// C-style output function -- forward to broken version in FiniteElement
715  /// until somebody decides what exactly they want to plot here...
716  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
717 
718  /// \short C-style output function -- forward to broken version in
719  /// FiniteElement until somebody decides what exactly they want to plot
720  /// here...
721  void output(FILE* file_pt, const unsigned &n_plot)
722  {FiniteElement::output(file_pt,n_plot);}
723 
724  /// \short The number of "DOF types" that degrees of freedom in this element
725  /// are sub-divided into: Just the solid degrees of freedom themselves.
726  unsigned ndof_types() const
727  {
728  return 1;
729  }
730 
731  /// \short Create a list of pairs for all unknowns in this element,
732  /// so that the first entry in each pair contains the global equation
733  /// number of the unknown, while the second one contains the number
734  /// of the "DOF types" that this unknown is associated with.
735  /// (Function can obviously only be called if the equation numbering
736  /// scheme has been set up.)
738  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const;
739 
740 
741 private:
742 
743  /// \short Normal vector to the clamping plane
745 
746 };
747 
748 
749 
750 }
751 
752 #endif
753 
754 
755 
756 
double load_rate_of_work()
Get integral of instantaneous rate of work done on the wall due to the load returned by the virtual f...
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1234
virtual void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Get the Jacobian and residuals. Wrapper to generic FSI version; that catches the case when we replace...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the element's contribution to its residual vector.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
static double Default_lambda_sq_value
Static default value for the timescale ratio (1.0 for natural scaling)
void enable_membrane_terms()
Set to renable the calculation of membrane terms (default)
This is a base class for all SolidFiniteElements that participate in FSI computations. These elements provide interfaces and generic funcionality for the two additional roles that SolidFiniteElements play in FSI problems:They parameterise the domain boundary for the fluid domain. To allow them to play this role, FSIWallElements are derived from the SolidFiniteElement and the GeomObject class, indicating that the every specific FSIWallElement must implement the pure virtual function GeomObject::position(...) which should compute the position vector to a point in the SolidFiniteElement, parametrised by its local coordinates.In FSI problems fluid exerts a traction onto the wall and this traction must be added to any other load terms (such as an external pressure acting on an elastic pipe) that are already applied to the SolidFiniteElements by other means.
Definition: fsi.h:210
void output_with_time_dep_quantities(std::ostream &outfile, const unsigned &n_plot)
Output position veloc and accel.
unsigned long nindex3() const
Return the range of index 3 of the derivatives of the shape functions.
Definition: shape.h:495
bool Normal_points_into_fluid
Boolean flag to indicate whether the normal is directed into the fluid.
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:2872
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Just the soli...
double * Nu_pt
Pointer to Poisson's ratio.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Just the soli...
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void fill_in_jacobian_from_external_interaction_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from all external interaction degrees of freedom (geometr...
static void Zero_traction_fct(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Return the jacobian is calculated by finite differences by default,.
void set_normal_pointing_out_of_fluid()
Set the normal computed by KirchhoffLoveShellEquations::get_normal(...) to point out of the fluid...
double *& lambda_sq_pt()
Return a pointer to timescale ratio (nondim density)
A general Finite Element class.
Definition: elements.h:1271
void(*&)(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load) load_vector_fct_pt()
Access to the load vector function pointer.
void set_normal_pointing_into_fluid()
Set the normal computed by KirchhoffLoveShellEquations::get_normal(...) to point into the fluid...
const double & h() const
Return the wall thickness to undeformed radius ratio.
void setup_fsi_wall_element(const unsigned &nlagr_solid, const unsigned &ndim_fluid)
Setup: Assign storage – pass the Eulerian dimension of the "adjacent" fluid elements and the number o...
Definition: fsi.cc:124
unsigned long nindex2() const
Return the range of index 2 of the derivatives of the shape functions.
Definition: shape.h:492
const double & lambda_sq() const
Return the timescale ratio (non-dimensional density)
void output(FILE *file_pt)
Overload the output function.
bool Ignore_membrane_terms
Boolean flag to ignore membrane terms.
void dposition_dlagrangian_at_local_coordinate(const Vector< double > &s, DenseMatrix< double > &drdxi) const
Derivative of position vector w.r.t. the SolidFiniteElement's.
static double Zero_prestress
Static default for prestress (set to zero)
double calculate_contravariant(double A[2][2], double Aup[2][2])
Invert a DIM by DIM matrix.
void output(std::ostream &outfile, const unsigned &n_plot)
Generic FiniteElement output function.
GeomObject * Undeformed_midplane_pt
Pointer to the GeomObject that specifies the undeformed midplane of the shell.
std::pair< double, double > get_strain_and_bend(const Vector< double > &s, DenseDoubleMatrix &strain, DenseDoubleMatrix &bend)
Get strain and bending tensors; returns pair comprising the determinant of the undeformed (*...
static double Default_nu_value
Static default value for the Poisson's ratio.
static double Default_h_value
Static default value for the thickness ratio.
void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Calculate the geometric shape functions at local coordinate s. Set any "superfluous" shape functions ...
void output(std::ostream &outfile)
Overload the output function.
bool Compute_rate_of_work_by_load_with_fluid_load_only
Boolean flag to indicate if rate-of-work by load is to be based on the fluid traction only...
DenseMatrix< double * > Prestress_pt
Pointer to membrane pre-stress terms – should probably generalise this to function pointers at some p...
static char t char * s
Definition: cfortran.h:572
void disable_membrane_terms()
Set to disable the calculation of membrane terms.
virtual void load_vector(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Get the load vector: Pass number of integration point (dummy), Lagr. coordinate and normal vector and...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Overload the standard fill in residuals contribution.
void shape(const Vector< double > &s, Shape &psi) const
Calculate the geometric shape functions at local coordinate s. Set any "superfluous" shape functions ...
void(* Load_vector_fct_pt)(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Pointer to load vector function: Its arguments are: Lagrangian coordinate, Eulerian coordinate...
void set_prestress_pt(const unsigned &i, const unsigned &j, double *value_pt)
Set pointer to (i,j)-th component of second Piola Kirchhoff membrane prestress to specified value (au...
void output(FILE *file_pt, const unsigned &n_plot)
Generic FiniteElement output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function – forward to broken version in FiniteElement until somebody decides what exac...
ClampedHermiteShellBoundaryConditionElement(const ClampedHermiteShellBoundaryConditionElement &dummy)
Broken copy constructor.
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:549
void fluid_load_vector(const unsigned &intpt, const Vector< double > &N, Vector< double > &load)
Get FE Jacobian by systematic finite differencing w.r.t. nodal positition Data, internal and external...
Definition: fsi.cc:166
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
unsigned nindex1() const
Return the range of index 1 of the shape function object.
Definition: shape.h:262
HermiteShellElement()
Constructor, there are no internal data points.
FSIDiagHermiteShellElement()
Constructor: Create shell element as FSIWallElement (and thus, by inheritance, a GeomObject) with two...
DiagHermiteShellElement()
Constructor, there are no internal data points.
double * H_pt
Pointer to dimensionless wall thickness.
void load_vector(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Get the load vector: Pass number of the integration point, Lagr. coordinate, Eulerian coordinate and ...
Vector< double > Normal_to_clamping_plane
Normal vector to the clamping plane.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
double prestress(const unsigned &i, const unsigned &j)
Return (i,j)-th component of second Piola Kirchhoff membrane prestress.
virtual void load_vector_for_rate_of_work_computation(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Get the load vector for the computation of the rate of work done by the load. Can switch between full...
void output(std::ostream &outfile)
Overload the output function.
void fill_in_contribution_to_residuals_shell(Vector< double > &residuals)
Return the residuals for the equations of KL shell theory.
const double & nu() const
Return the Poisson's ratio.
void set_nodal_dimension(const unsigned &nodal_dim)
Set the dimension of the nodes in the element. This will typically only be required when constructing...
Definition: elements.h:1348
double *& nu_pt()
Return a pointer to the Poisson ratio.
double *& h_pt()
Return a pointer to the non-dim wall thickness.
double fluid_load_rate_of_work()
Get integral of instantaneous rate of work done on the wall due to the fluid load returned by the fun...
ClampedHermiteShellBoundaryConditionElement()
Broken empty constructor.
GeomObject *& undeformed_midplane_pt()
Return a reference to geometric object that specifies the shell's undeformed geometry.
KirchhoffLoveShellEquations()
Constructor: Initialise physical parameter values to defaults.
unsigned long nindex1() const
Return the range of index 1 of the derivatives of the shape functions.
Definition: shape.h:489
double * Lambda_sq_pt
Pointer to timescale ratio (non-dimensional density)
void output(std::ostream &outfile)
Generic FiniteElement output function.
virtual void load_vector_for_rate_of_work_computation(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Get the load vector for the computation of the rate of work done by the load. Here we simply forward ...
void operator=(const ClampedHermiteShellBoundaryConditionElement &)
Broken assignment operator.
FaceGeometry()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
void get_energy(double &pot_en, double &kin_en)
Get potential (strain) and kinetic energy of the element.
void output(FILE *file_pt)
Generic FiniteElement output function.
SolidFiniteElement class.
Definition: elements.h:3320
unsigned nindex2() const
Return the range of index 2 of the shape function object.
Definition: shape.h:265
~FSIDiagHermiteShellElement()
Destructor: empty.
void get_normal(const Vector< double > &s, Vector< double > &N)
Get normal vector.
void set_symmetry_line(const Vector< double > &normal_to_clamping_plane)
Set normal vector to clamping plane.
void resize(const unsigned long &n)
Definition: matrices.h:511