linear_elasticity_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: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
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 general linear elasticity elements
31 
32 //Include guards to prevent multiple inclusion of the header
33 #ifndef OOMPH_LINEAR_ELASTICITY_ELEMENTS_HEADER
34 #define OOMPH_LINEAR_ELASTICITY_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 
42 //OOMPH-LIB headers
43 #include "../generic/Qelements.h"
44 #include "../generic/mesh.h"
45 #include "../generic/hermite_elements.h"
46 #include "./elasticity_tensor.h"
47 #include "../generic/projection.h"
48 
49 namespace oomph
50 {
51 //=======================================================================
52 /// A base class for elements that solve the equations of linear
53 /// elasticity in Cartesian coordinates.
54 /// Combines a few generic functions that are shared by
55 /// LinearElasticityEquations
56 /// and LinearElasticityEquationsWithPressure (hierher: The latter
57 /// don't exist yet but will be written as soon as somebody needs them...)
58 //=======================================================================
59  template <unsigned DIM>
61  {
62  public:
63 
64  /// \short Return the index at which the i-th unknown displacement
65  /// component is stored. The default value, i, is appropriate for
66  /// single-physics problems.
67  virtual inline unsigned u_index_linear_elasticity(const unsigned i) const
68  {return i;}
69 
70  /// d^2u/dt^2 at local node n
71  double d2u_dt2_linear_elasticity(const unsigned &n,
72  const unsigned &i) const
73  {
74  // Get the timestepper
76 
77  // Storage for the derivative - initialise to 0
78  double d2u_dt2=0.0;
79 
80  // If we are doing an unsteady solve then calculate the derivative
81  if(!time_stepper_pt->is_steady())
82  {
83  // Get the nodal index
84  const unsigned u_nodal_index=u_index_linear_elasticity(i);
85 
86  // Get the number of values required to represent history
87  const unsigned n_time=time_stepper_pt->ntstorage();
88 
89  // Loop over history values
90  for(unsigned t=0;t<n_time;t++)
91  {
92  //Add the contribution to the derivative
93  d2u_dt2+=time_stepper_pt->weight(2,t)*nodal_value(t,n,u_nodal_index);
94  }
95  }
96 
97  return d2u_dt2;
98  }
99 
100  /// Compute vector of FE interpolated displacement u at local coordinate s
102  Vector<double>& disp)
103  const
104  {
105  //Find number of nodes
106  unsigned n_node = nnode();
107 
108  //Local shape function
109  Shape psi(n_node);
110 
111  //Find values of shape function
112  shape(s,psi);
113 
114  for (unsigned i=0;i<DIM;i++)
115  {
116  //Index at which the nodal value is stored
117  unsigned u_nodal_index = u_index_linear_elasticity(i);
118 
119  //Initialise value of u
120  disp[i] = 0.0;
121 
122  //Loop over the local nodes and sum
123  for(unsigned l=0;l<n_node;l++)
124  {
125  disp[i] += nodal_value(l,u_nodal_index)*psi[l];
126  }
127  }
128  }
129 
130  /// Return FE interpolated displacement u[i] at local coordinate s
132  const unsigned &i) const
133  {
134  //Find number of nodes
135  unsigned n_node = nnode();
136 
137  //Local shape function
138  Shape psi(n_node);
139 
140  //Find values of shape function
141  shape(s,psi);
142 
143  //Get nodal index at which i-th velocity is stored
144  unsigned u_nodal_index = u_index_linear_elasticity(i);
145 
146  //Initialise value of u
147  double interpolated_u = 0.0;
148 
149  //Loop over the local nodes and sum
150  for(unsigned l=0;l<n_node;l++)
151  {
152  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
153  }
154 
155  return(interpolated_u);
156  }
157 
158 
159  /// \short Function pointer to function that specifies the body force
160  /// as a function of the Cartesian coordinates and time FCT(t,x,b) --
161  /// x and b are Vectors!
162  typedef void (*BodyForceFctPt)(const double& t,
163  const Vector<double>& x,
164  Vector<double>& b);
165 
166  /// \short Constructor: Set null pointers for constitutive law and for
167  /// isotropic growth function. Set physical parameter values to
168  /// default values, switch on inertia and set body force to zero.
171  Body_force_fct_pt(0) {}
172 
173  /// Return the pointer to the elasticity_tensor
175 
176  /// Access function to the entries in the elasticity tensor
177  inline double E(const unsigned &i,const unsigned &j,
178  const unsigned &k, const unsigned &l) const
179  {
180  return (*Elasticity_tensor_pt)(i,j,k,l);
181  }
182 
183  ///Access function for timescale ratio (nondim density)
184  const double& lambda_sq() const {return *Lambda_sq_pt;}
185 
186  /// Access function for pointer to timescale ratio (nondim density)
187  double* &lambda_sq_pt() {return Lambda_sq_pt;}
188 
189  /// Access function: Pointer to body force function
191 
192  /// Access function: Pointer to body force function (const version)
194 
195 
196  /// Switch on solid inertia
197  void enable_inertia() {Unsteady=true;}
198 
199  /// Switch off solid inertia
200  void disable_inertia() {Unsteady=false;}
201 
202  ///Access function to flag that switches inertia on/off (const version)
203  bool is_inertia_enabled() const {return Unsteady;}
204 
205  ///Pin the element's redundant solid pressures (needed for refinement)
207 
208  /// \short Loop over all elements in Vector (which typically contains
209  /// all the elements in a refineable solid mesh) and pin the nodal solid
210  /// pressure degrees of freedom that are not being used. Function uses
211  /// the member function
212  /// - \c LinearElasticityEquationsBase<DIM>::
213  /// pin_elemental_redundant_nodal_pressure_dofs()
214  /// .
215  /// which is empty by default and should be implemented for
216  /// elements with nodal solid pressure degrees of freedom
217  /// (e.g. linear elasticity elements with continuous pressure interpolation.)
219  const Vector<GeneralisedElement*>& element_pt)
220  {
221  // Loop over all elements and call the function that pins their
222  // unused nodal solid pressure data
223  unsigned n_element = element_pt.size();
224  for(unsigned e=0;e<n_element;e++)
225  {
226  dynamic_cast<LinearElasticityEquationsBase<DIM>*>(element_pt[e])->
228  }
229  }
230 
231  /// \short Return the Cauchy stress tensor, as calculated
232  /// from the elasticity tensor at specified local coordinate
233  /// Virtual so separaete versions can (and must!) be provided
234  /// for displacement and pressure-displacement formulations.
235  virtual void get_stress(const Vector<double> &s,
236  DenseMatrix<double> &sigma) const=0;
237 
238  /// \short Return the strain tensor
239  void get_strain(const Vector<double> &s, DenseMatrix<double> &strain) const;
240 
241  /// \short Evaluate body force at Eulerian coordinate x at present time
242  /// (returns zero vector if no body force function pointer has been set)
243  inline void body_force(const Vector<double>& x,
244  Vector<double>& b) const
245  {
246  //If no function has been set, return zero vector
247  if(Body_force_fct_pt==0)
248  {
249  // Get spatial dimension of element
250  unsigned n=dim();
251  for (unsigned i=0;i<n;i++)
252  {
253  b[i] = 0.0;
254  }
255  }
256  else
257  {
258  // Get time from timestepper of first node (note that this must
259  // work -- body force only makes sense for elements that can be
260  // deformed and given that the deformation of solid finite elements
261  // is controlled by their nodes, nodes must exist!)
262  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
263 
264  // Now evaluate the body force
265  (*Body_force_fct_pt)(time,x,b);
266  }
267  }
268 
269 
270 
271 
272  /// \short The number of "DOF types" that degrees of freedom in this element
273  /// are sub-divided into: for now lump them all into one DOF type.
274  /// Can be adjusted later
275  unsigned ndof_types() const
276  {
277  return 2;
278  //return 1;
279  }
280 
281  /// \short Create a list of pairs for all unknowns in this element,
282  /// so that the first entry in each pair contains the global equation
283  /// number of the unknown, while the second one contains the number
284  /// of the "DOF types" that this unknown is associated with.
285  /// (Function can obviously only be called if the equation numbering
286  /// scheme has been set up.)
288  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
289  {
290 
291  // temporary pair (used to store dof lookup prior to being added
292  // to list)
293  std::pair<unsigned long,unsigned> dof_lookup;
294 
295  // number of nodes
296  const unsigned n_node = this->nnode();
297 
298  //Integer storage for local unknown
299  int local_unknown=0;
300 
301  //Loop over the nodes
302  for(unsigned n=0;n<n_node;n++)
303  {
304  //Loop over dimension
305  for(unsigned i=0;i<DIM;i++)
306  {
307  //If the variable is free
308  local_unknown = nodal_local_eqn(n,i);
309 
310  // ignore pinned values
311  if (local_unknown >= 0)
312  {
313  // store dof lookup in temporary pair: First entry in pair
314  // is global equation number; second entry is dof type
315  dof_lookup.first = this->eqn_number(local_unknown);
316  dof_lookup.second = i;
317  //dof_lookup.second = DIM;
318 
319  // add to list
320  dof_lookup_list.push_front(dof_lookup);
321 
322  }
323  }
324  }
325  }
326 
327 
328  protected:
329 
330  /// Pointer to the elasticity tensor
332 
333  /// Timescale ratio (non-dim. density)
334  double* Lambda_sq_pt;
335 
336  /// Flag that switches inertia on/off
337  bool Unsteady;
338 
339  /// Pointer to body force function
341 
342  /// Static default value for timescale ratio (1.0 -- for natural scaling)
343  static double Default_lambda_sq_value;
344 
345  };
346 
347 
348 ///////////////////////////////////////////////////////////////////////
349 ///////////////////////////////////////////////////////////////////////
350 ///////////////////////////////////////////////////////////////////////
351 
352 
353 //=======================================================================
354 /// A class for elements that solve the equations of linear elasticity
355 /// in cartesian coordinates.
356 //=======================================================================
357  template <unsigned DIM>
359  {
360  public:
361 
362  /// \short Constructor
364 
365  /// Number of values required at node n.
366  unsigned required_nvalue(const unsigned &n) const {return DIM;}
367 
368  /// \short Return the residuals for the solid equations (the discretised
369  /// principle of virtual displacements)
371  {
374  }
375 
376  /// The jacobian is calculated by finite differences by default,
377  /// We need only to take finite differences w.r.t. positional variables
378  /// For this element
380  DenseMatrix<double> &jacobian)
381  {
382  //Add the contribution to the residuals
384  residuals,jacobian,1);
385  }
386 
387  /// \short Return the Cauchy stress tensor, as calculated
388  /// from the elasticity tensor at specified local coordinate
389  void get_stress(const Vector<double> &s,
390  DenseMatrix<double> &sigma) const;
391 
392 
393  ///Output exact solution x,y,[z],u,v,[w]
394  void output_fct(std::ostream &outfile,
395  const unsigned &nplot,
397 
398  ///Output exact solution x,y,[z],u,v,[w] (unsteady version)
399  void output_fct(std::ostream &outfile,
400  const unsigned &nplot,
401  const double &time,
403 
404  /// Output: x,y,[z],u,v,[w]
405  void output(std::ostream &outfile)
406  {
407  unsigned n_plot=5;
408  output(outfile,n_plot);
409  }
410 
411  /// Output: x,y,[z],u,v,[w]
412  void output(std::ostream &outfile, const unsigned &n_plot);
413 
414 
415  /// C-style output: x,y,[z],u,v,[w]
416  void output(FILE* file_pt)
417  {
418  unsigned n_plot=5;
419  output(file_pt,n_plot);
420  }
421 
422  /// Output: x,y,[z],u,v,[w]
423  void output(FILE* file_pt, const unsigned &n_plot);
424 
425  /// \short Validate against exact solution.
426  /// Solution is provided via function pointer.
427  /// Plot at a given number of plot points and compute L2 error
428  /// and L2 norm of displacement solution over element
429  void compute_error(std::ostream &outfile,
431  double& error, double& norm);
432 
433  /// \short Validate against exact solution.
434  /// Solution is provided via function pointer.
435  /// Plot at a given number of plot points and compute L2 error
436  /// and L2 norm of displacement solution over element
437  /// (unsteady version)
438  void compute_error(std::ostream &outfile,
440  const double& time, double& error, double& norm);
441 
442  private:
443 
444 
445  /// \short Private helper function to compute residuals and (if requested
446  /// via flag) also the Jacobian matrix.
448  Vector<double> &residuals,DenseMatrix<double> &jacobian,unsigned flag);
449 
450  };
451 
452 
453 ////////////////////////////////////////////////////////////////////////
454 ////////////////////////////////////////////////////////////////////////
455 ////////////////////////////////////////////////////////////////////////
456 
457 
458 //===========================================================================
459 /// An Element that solves the equations of linear elasticity
460 /// in Cartesian coordinates, using QElements for the geometry
461 //============================================================================
462  template<unsigned DIM, unsigned NNODE_1D>
463  class QLinearElasticityElement : public virtual QElement<DIM,NNODE_1D>,
464  public virtual LinearElasticityEquations<DIM>
465  {
466  public:
467 
468  /// Constructor
469  QLinearElasticityElement() : QElement<DIM,NNODE_1D>(),
470  LinearElasticityEquations<DIM>() { }
471 
472  ///Output exact solution x,y,[z],u,v,[w]
473  void output_fct(std::ostream &outfile,
474  const unsigned &nplot,
476  {
477  LinearElasticityEquations<DIM>::output_fct(outfile,nplot,exact_soln_pt);
478  }
479 
480  /// Output function
481  void output(std::ostream &outfile)
483 
484  /// Output function
485  void output(std::ostream &outfile, const unsigned &n_plot)
486  {LinearElasticityEquations<DIM>::output(outfile,n_plot);}
487 
488 
489  /// C-style output function
490  void output(FILE* file_pt)
492 
493  /// C-style output function
494  void output(FILE* file_pt, const unsigned &n_plot)
495  {LinearElasticityEquations<DIM>::output(file_pt,n_plot);}
496 
497  };
498 
499 
500 //============================================================================
501 /// FaceGeometry of a linear 2D QLinearElasticityElement element
502 //============================================================================
503  template<>
505  public virtual QElement<1,2>
506  {
507  public:
508  /// Constructor must call the constructor of the underlying solid element
509  FaceGeometry() : QElement<1,2>() {}
510  };
511 
512 
513 
514 //============================================================================
515 /// FaceGeometry of a quadratic 2D QLinearElasticityElement element
516 //============================================================================
517  template<>
519  public virtual QElement<1,3>
520  {
521  public:
522  /// Constructor must call the constructor of the underlying element
523  FaceGeometry() : QElement<1,3>() {}
524  };
525 
526 
527 
528 //============================================================================
529 /// FaceGeometry of a cubic 2D QLinearElasticityElement element
530 //============================================================================
531  template<>
533  public virtual QElement<1,4>
534  {
535  public:
536  /// Constructor must call the constructor of the underlying element
537  FaceGeometry() : QElement<1,4>() {}
538  };
539 
540 
541 //============================================================================
542 /// FaceGeometry of a linear 3D QLinearElasticityElement element
543 //============================================================================
544  template<>
546  public virtual QElement<2,2>
547  {
548  public:
549  /// Constructor must call the constructor of the underlying element
550  FaceGeometry() : QElement<2,2>() {}
551  };
552 
553 //============================================================================
554 /// FaceGeometry of a quadratic 3D QLinearElasticityElement element
555 //============================================================================
556  template<>
558  public virtual QElement<2,3>
559  {
560  public:
561  /// Constructor must call the constructor of the underlying element
562  FaceGeometry() : QElement<2,3>() {}
563  };
564 
565 
566 //============================================================================
567 /// FaceGeometry of a cubic 3D QLinearElasticityElement element
568 //============================================================================
569  template<>
571  public virtual QElement<2,4>
572  {
573  public:
574  /// Constructor must call the constructor of the underlying element
575  FaceGeometry() : QElement<2,4>() {}
576  };
577 
578 ////////////////////////////////////////////////////////////////////
579 ////////////////////////////////////////////////////////////////////
580 ////////////////////////////////////////////////////////////////////
581 
582 
583 //==========================================================
584 /// Linear elasticity upgraded to become projectable
585 //==========================================================
586  template<class LINEAR_ELAST_ELEMENT>
588  public virtual ProjectableElement<LINEAR_ELAST_ELEMENT>
589  {
590 
591  public:
592 
593  /// \short Constructor [this was only required explicitly
594  /// from gcc 4.5.2 onwards...]
596 
597 
598  /// \short Specify the values associated with field fld.
599  /// The information is returned in a vector of pairs which comprise
600  /// the Data object and the value within it, that correspond to field fld.
601  /// In the underlying linear elasticity elements the
602  /// the displacements are stored at the nodal values
604  {
605  // Create the vector
606  Vector<std::pair<Data*,unsigned> > data_values;
607 
608  // Loop over all nodes and extract the fld-th nodal value
609  unsigned nnod=this->nnode();
610  for (unsigned j=0;j<nnod;j++)
611  {
612  // Add the data value associated with the displacement components
613  data_values.push_back(std::make_pair(this->node_pt(j),fld));
614  }
615 
616  // Return the vector
617  return data_values;
618 
619  }
620 
621  /// \short Number of fields to be projected: dim, corresponding to
622  /// the displacement components
624  {
625  return this->dim();
626  }
627 
628  /// \short Number of history values to be stored for fld-th field.
629  /// (includes present value!)
630  unsigned nhistory_values_for_projection(const unsigned &fld)
631  {
632 #ifdef PARANOID
633  if (fld>1)
634  {
635  std::stringstream error_stream;
636  error_stream
637  << "Elements only store two fields so fld can't be"
638  << " " << fld << std::endl;
639  throw OomphLibError(
640  error_stream.str(),
641  OOMPH_CURRENT_FUNCTION,
642  OOMPH_EXCEPTION_LOCATION);
643  }
644 #endif
645  return this->node_pt(0)->ntstorage();
646  }
647 
648  ///\short Number of positional history values: Read out from
649  /// positional timestepper (Note: count includes current value!)
651  {
652  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
653  }
654 
655  /// \short Return Jacobian of mapping and shape functions of field fld
656  /// at local coordinate s
657  double jacobian_and_shape_of_field(const unsigned &fld,
658  const Vector<double> &s,
659  Shape &psi)
660  {
661  unsigned n_dim=this->dim();
662  unsigned n_node=this->nnode();
663  DShape dpsidx(n_node,n_dim);
664 
665  // Call the derivatives of the shape functions and return
666  // the Jacobian
667  return this->dshape_eulerian(s,psi,dpsidx);
668  }
669 
670 
671 
672  /// \short Return interpolated field fld at local coordinate s, at time level
673  /// t (t=0: present; t>0: history values)
674  double get_field(const unsigned &t,
675  const unsigned &fld,
676  const Vector<double>& s)
677  {
678  unsigned n_node=this->nnode();
679 
680 /* #ifdef PARANOID */
681 /* unsigned n_dim=this->node_pt(0)->ndim(); */
682 /* #endif */
683 
684  //Local shape function
685  Shape psi(n_node);
686 
687  //Find values of shape function
688  this->shape(s,psi);
689 
690  //Initialise value of u
691  double interpolated_u = 0.0;
692 
693  //Sum over the local nodes at that time
694  for(unsigned l=0;l<n_node;l++)
695  {
696 // over-zealous I think. This will quietly do the right thing
697 // even if there are additional degrees of freedom floating around.
698 /* #ifdef PARANOID */
699 /* unsigned nvalue=this->node_pt(l)->nvalue(); */
700 /* if (nvalue!=n_dim) */
701 /* { */
702 /* std::stringstream error_stream; */
703 /* error_stream */
704 /* << "Current implementation only works for non-resized nodes\n" */
705 /* << "but nvalue= " << nvalue << "!= dim = " << n_dim << std::endl; */
706 /* throw OomphLibError( */
707 /* error_stream.str(), */
708 /* OOMPH_CURRENT_FUNCTION, */
709 /* OOMPH_EXCEPTION_LOCATION); */
710 /* } */
711 /* #endif */
712  interpolated_u += this->nodal_value(t,l,fld)*psi[l];
713  }
714  return interpolated_u;
715  }
716 
717 
718  ///Return number of values in field fld
719  unsigned nvalue_of_field(const unsigned &fld)
720  {
721  return this->nnode();
722  }
723 
724 
725  ///Return local equation number of value j in field fld.
726  int local_equation(const unsigned &fld,
727  const unsigned &j)
728  {
729 // over-zealous I think. This will quietly do the right thing
730 // even if there are additional degrees of freedom floating around.
731 /* #ifdef PARANOID */
732 /* unsigned n_dim=this->node_pt(0)->ndim(); */
733 /* unsigned nvalue=this->node_pt(j)->nvalue(); */
734 /* if (nvalue!=n_dim) */
735 /* { */
736 /* std::stringstream error_stream; */
737 /* error_stream */
738 /* << "Current implementation only works for non-resized nodes\n" */
739 /* << "but nvalue= " << nvalue << "!= dim = " << n_dim << std::endl; */
740 /* throw OomphLibError( */
741 /* error_stream.str(), */
742 /* OOMPH_CURRENT_FUNCTION, */
743 /* OOMPH_EXCEPTION_LOCATION); */
744 /* } */
745 /* #endif */
746  return this->nodal_local_eqn(j,fld);
747  }
748 
749 
750  };
751 
752 
753 //=======================================================================
754 /// Face geometry for element is the same as that for the underlying
755 /// wrapped element
756 //=======================================================================
757  template<class ELEMENT>
759  : public virtual FaceGeometry<ELEMENT>
760  {
761  public:
762  FaceGeometry() : FaceGeometry<ELEMENT>() {}
763  };
764 
765 
766 //=======================================================================
767 /// Face geometry of the Face Geometry for element is the same as
768 /// that for the underlying wrapped element
769 //=======================================================================
770  template<class ELEMENT>
772  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
773  {
774  public:
776  };
777 
778 
779 }
780 
781 #endif
782 
783 
784 
785 
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: for now lump ...
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
virtual void fill_in_generic_contribution_to_residuals_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Private helper function to compute residuals and (if requested via flag) also the Jacobian matrix...
virtual unsigned u_index_linear_elasticity(const unsigned i) const
Return the index at which the i-th unknown displacement component is stored. The default value...
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
FaceGeometry()
Constructor must call the constructor of the underlying element.
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
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
double d2u_dt2_linear_elasticity(const unsigned &n, const unsigned &i) const
d^2u/dt^2 at local node n
ProjectableLinearElasticityElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
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 jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
virtual void pin_elemental_redundant_nodal_solid_pressures()
Pin the element's redundant solid pressures (needed for refinement)
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
void(* BodyForceFctPt)(const double &t, const Vector< double > &x, Vector< double > &b)
Function pointer to function that specifies the body force as a function of the Cartesian coordinates...
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
void enable_inertia()
Switch on solid inertia.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (includes present value!) ...
ElasticityTensor * Elasticity_tensor_pt
Pointer to the elasticity tensor.
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
unsigned nfields_for_projection()
Number of fields to be projected: dim, corresponding to the displacement components.
LinearElasticityEquationsBase()
Constructor: Set null pointers for constitutive law and for isotropic growth function. Set physical parameter values to default values, switch on inertia and set body force to zero.
e
Definition: cfortran.h:575
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Validate against exact solution. Solution is provided via function pointer. Plot at a given number of...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma) const
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...
FaceGeometry()
Constructor must call the constructor of the underlying element.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
static void pin_redundant_nodal_solid_pressures(const Vector< GeneralisedElement * > &element_pt)
Loop over all elements in Vector (which typically contains all the elements in a refineable solid mes...
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
bool Unsteady
Flag that switches inertia on/off.
void output(FILE *file_pt)
C-style output: x,y,[z],u,v,[w].
void disable_inertia()
Switch off solid inertia.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values: Read out from positional timestepper (Note: count includes curre...
void interpolated_u_linear_elasticity(const Vector< double > &s, Vector< double > &disp) const
Compute vector of FE interpolated displacement u at local coordinate s.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals for the solid equations (the discretised principle of virtual displacements) ...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution x,y,[z],u,v,[w].
virtual void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma) const =0
Return the Cauchy stress tensor, as calculated from the elasticity tensor at specified local coordina...
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution x,y,[z],u,v,[w].
static char t char * s
Definition: cfortran.h:572
void output(FILE *file_pt)
C-style output function.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double E(const unsigned &i, const unsigned &j, const unsigned &k, const unsigned &l) const
Access function to the entries in the elasticity tensor.
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
BodyForceFctPt body_force_fct_pt() const
Access function: Pointer to body force function (const version)
void body_force(const Vector< double > &x, Vector< double > &b) const
Evaluate body force at Eulerian coordinate x at present time (returns zero vector if no body force fu...
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
FaceGeometry()
Constructor must call the constructor of the underlying element.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
FaceGeometry()
Constructor must call the constructor of the underlying element.
void output(std::ostream &outfile)
Output function.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
bool is_inertia_enabled() const
Access function to flag that switches inertia on/off (const version)
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
double interpolated_u_linear_elasticity(const Vector< double > &s, const unsigned &i) const
Return FE interpolated displacement u[i] at local coordinate s.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
void output(std::ostream &outfile)
Output: x,y,[z],u,v,[w].
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
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...
FaceGeometry()
Constructor must call the constructor of the underlying element.
Linear elasticity upgraded to become projectable.
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
BodyForceFctPt Body_force_fct_pt
Pointer to body force function.
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 – for natural scaling)
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
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...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
unsigned required_nvalue(const unsigned &n) const
Number of values required at node n.
ElasticityTensor *& elasticity_tensor_pt()
Return the pointer to the elasticity_tensor.