solid_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 May 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 solid mechanics elements
31 
32 //Include guards to prevent multiple inclusion of the header
33 #ifndef OOMPH_ELASTICITY_ELEMENTS_HEADER
34 #define OOMPH_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/Telements.h"
45 #include "../generic/mesh.h"
46 #include "../generic/hermite_elements.h"
47 #include "../constitutive/constitutive_laws.h"
48 #include "../generic/error_estimator.h"
49 #include "../generic/projection.h"
50 
51 
52 namespace oomph
53 {
54 
55 
56 //=======================================================================
57 /// A base class for elements that solve the equations of solid mechanics,
58 /// based on the principle of virtual displacements in Cartesian coordinates.
59 /// Combines a few generic functions that are shared by PVDEquations
60 /// and PVDEquationsWithPressure.
61 //=======================================================================
62  template <unsigned DIM>
63  class PVDEquationsBase : public virtual SolidFiniteElement
64  {
65  private:
66 
67  /// \short Static "magic" number that indicates that the solid pressure
68  /// is not stored at a node
70 
71  public:
72 
73  /// \short Function pointer to function that specifies the isotropic
74  /// growth as a function of the Lagrangian coordinates FCT(xi,gamma(xi)) --
75  /// xi is a Vector!
76  typedef void (*IsotropicGrowthFctPt)
77  (const Vector<double>& xi, double& gamma);
78 
79  /// \short Function pointer to function that specifies the pre-stress
80  /// sigma_0(i,j) as a function of the Lagrangian coordinates
81  /// FCT(i,j,xi) -- xi is a Vector!
82  typedef double (*PrestressFctPt)
83  (const unsigned& i, const unsigned& j, const Vector<double>& xi);
84 
85  /// \short Function pointer to function that specifies the body force
86  /// as a function of the Lagrangian coordinates and time FCT(t,xi,b) --
87  /// xi and b are Vectors!
88  typedef void (*BodyForceFctPt)(const double& t,
89  const Vector<double>& xi,
90  Vector<double>& b);
91 
92  /// \short Constructor: Set null pointers for constitutive law and for
93  /// isotropic growth function. Set physical parameter values to
94  /// default values, enable inertia and set body force to zero.
95  /// Default evaluation of Jacobian: analytically rather than by FD.
100  {}
101 
102  /// Return the constitutive law pointer
104 
105 
106  ///Access function for timescale ratio (nondim density)
107  const double& lambda_sq() const {return *Lambda_sq_pt;}
108 
109 
110  /// Access function for pointer to timescale ratio (nondim density)
111  double* &lambda_sq_pt() {return Lambda_sq_pt;}
112 
113 
114  /// Access function: Pointer to isotropic growth function
116  {return Isotropic_growth_fct_pt;}
117 
118  /// Access function: Pointer to pre-stress function
120  {return Prestress_fct_pt;}
121 
122  /// Access function: Pointer to isotropic growth function (const version)
124  {return Isotropic_growth_fct_pt;}
125 
126  /// Access function: Pointer to body force function
128 
129  /// Access function: Pointer to body force function (const version)
131 
132  /// Switch on solid inertia
133  void enable_inertia() {Unsteady=true;}
134 
135  /// Switch off solid inertia
136  void disable_inertia() {Unsteady=false;}
137 
138  ///Access function to flag that switches inertia on/off (const version)
139  bool is_inertia_enabled() const {return Unsteady;}
140 
141  /// \short Return the number of solid pressure degrees of freedom
142  /// Default is that there are no solid pressures
143  virtual unsigned npres_solid() const {return 0;}
144 
145  /// \short Return the local degree of freedom associated with the
146  /// i-th solid pressure. Default is that there are none.
147  virtual int solid_p_local_eqn(const unsigned &i)const {return -1;}
148 
149  /// \short Return the index at which the solid pressure is stored if it
150  /// is stored at the nodes. If not stored at the nodes this will return
151  /// a negative number.
152  virtual int solid_p_nodal_index() const
154 
155 
156  /// \short Unpin all solid pressure dofs in the element
157  virtual void unpin_elemental_solid_pressure_dofs()=0;
158 
159  ///Pin the element's redundant solid pressures (needed for refinement)
161 
162  /// \short Loop over all elements in Vector (which typically contains
163  /// all the elements in a refineable solid mesh) and pin the nodal solid
164  /// pressure degrees of freedom that are not being used. Function uses
165  /// the member function
166  /// - \c PVDEquationsBase<DIM>::
167  /// pin_elemental_redundant_nodal_pressure_dofs()
168  /// .
169  /// which is empty by default and should be implemented for
170  /// elements with nodal solid pressure degrees of freedom
171  /// (e.g. solid elements with continuous pressure interpolation.)
173  const Vector<GeneralisedElement*>& element_pt)
174  {
175  // Loop over all elements
176  unsigned n_element = element_pt.size();
177  for(unsigned e=0;e<n_element;e++)
178  {
179  dynamic_cast<PVDEquationsBase<DIM>*>(element_pt[e])->
181  }
182  }
183 
184  /// \short Unpin all pressure dofs in elements listed in vector.
186  element_pt)
187  {
188  // Loop over all elements
189  unsigned n_element = element_pt.size();
190  for(unsigned e=0;e<n_element;e++)
191  {
192  dynamic_cast<PVDEquationsBase<DIM>*>(element_pt[e])->
194  }
195  }
196 
197  /// \short Return the 2nd Piola Kirchoff stress tensor, as calculated
198  /// from the constitutive law at specified local coordinate
199  /// (needed by \c get_principal_stress(...), so I'm afraid I will
200  /// have to insist that you implement it...
201  virtual void get_stress(const Vector<double> &s,
202  DenseMatrix<double> &sigma)=0;
203 
204  /// \short Return the strain tensor
205  void get_strain(const Vector<double> &s, DenseMatrix<double> &strain) const;
206 
207  /// Get potential (strain) and kinetic energy
208  void get_energy(double &pot_en, double &kin_en);
209 
210  /// \short Return the deformed covariant basis vectors
211  /// at specified local coordinate: \c def_covariant_basis(i,j)
212  /// is the j-th component of the i-th basis vector.
215  def_covariant_basis);
216 
217 
218  /// \short Compute principal stress vectors and (scalar) principal stresses
219  /// at specified local coordinate. \c principal_stress_vector(i,j)
220  /// is the j-th component of the i-th principal stress vector.
222  DenseMatrix<double>& principal_stress_vector,
223  Vector<double>& principal_stress);
224 
225 
226  /// \short Evaluate isotropic growth function at Lagrangian coordinate xi
227  /// and/or local coordinate s.
228  /// (returns 1, i.e. no growth, if no function pointer has been set)
229  /// This function is virtual to allow overloading in multi-physics
230  /// problems where the growth function might be determined by
231  /// another system of equations
232  virtual inline void get_isotropic_growth(const unsigned& ipt,
233  const Vector<double> &s,
234  const Vector<double>& xi,
235  double& gamma) const
236  {
237  //If no function has been set, return 1
238  if(Isotropic_growth_fct_pt==0) {gamma = 1.0;}
239  else
240  {
241  // Get isotropic growth
242  (*Isotropic_growth_fct_pt)(xi,gamma);
243  }
244  }
245 
246 
247  /// \short Evaluate body force at Lagrangian coordinate xi at present time
248  /// (returns zero vector if no body force function pointer has been set)
249  inline void body_force(const Vector<double>& xi,
250  Vector<double>& b) const
251  {
252  //If no function has been set, return zero vector
253  if(Body_force_fct_pt==0)
254  {
255  // Get spatial dimension of element
256  unsigned n=dim();
257  for (unsigned i=0;i<n;i++)
258  {
259  b[i] = 0.0;
260  }
261  }
262  else
263  {
264  // Get time from timestepper of first node (note that this must
265  // work -- body force only makes sense for elements that can be
266  // deformed and given that the deformation of solid finite elements
267  // is controlled by their nodes, nodes must exist!)
268  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
269 
270  // Now evaluate the body force
271  (*Body_force_fct_pt)(time,xi,b);
272  }
273 
274  }
275 
276 
277 
278 
279  /// \short returns the number of DOF types associated with this element.
280  unsigned ndof_types() const
281  {
282  return DIM;
283  }
284 
285  /// \short Create a list of pairs for all unknowns in this element,
286  /// so that the first entry in each pair contains the global equation
287  /// number of the unknown, while the second one contains the number
288  /// of the "DOF" that this unknown is associated with.
289  /// (Function can obviously only be called if the equation numbering
290  /// scheme has been set up.)
291  /// E.g. in a 3D problem there are 3 types of DOF:
292  /// 0 - x displacement
293  /// 1 - y displacement
294  /// 2 - z displacement
296  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
297  {
298  // temporary pair (used to store dof lookup prior to being added to list
299  std::pair<unsigned,unsigned> dof_lookup;
300 
301  // number of nodes
302  const unsigned n_node = this->nnode();
303 
304  //Get the number of position dofs and dimensions at the node
305  const unsigned n_position_type = nnodal_position_type();
306  const unsigned nodal_dim = nodal_dimension();
307 
308  //Integer storage for local unknown
309  int local_unknown=0;
310 
311  //Loop over the nodes
312  for(unsigned n=0;n<n_node;n++)
313  {
314 
315  //Loop over position dofs
316  for(unsigned k=0;k<n_position_type;k++)
317  {
318  //Loop over dimension
319  for(unsigned i=0;i<nodal_dim;i++)
320  {
321  //If the variable is free
322  local_unknown = position_local_eqn(n,k,i);
323 
324  // ignore pinned values
325  if (local_unknown >= 0)
326  {
327  // store dof lookup in temporary pair: First entry in pair
328  // is global equation number; second entry is dof type
329  dof_lookup.first = this->eqn_number(local_unknown);
330  dof_lookup.second = i;
331 
332  // add to list
333  dof_lookup_list.push_front(dof_lookup);
334 
335  }
336  }
337  }
338  }
339  }
340 
341  /// Set Jacobian to be evaluated by FD? Else: Analytically.
343 
344  /// Set Jacobian to be evaluated analytically Else: by FD
346 
347  /// Return the flag indicating whether the jacobian is evaluated by fd
349 
350  /// \short Return (i,j)-th component of second Piola Kirchhoff membrane
351  /// prestress at Lagrangian coordinate xi
352  double prestress(const unsigned& i,
353  const unsigned& j,
354  const Vector<double> xi)
355  {
356  if (Prestress_fct_pt==0)
357  {
358  return 0.0;
359  }
360  else
361  {
362  return (*Prestress_fct_pt)(i,j,xi);
363  }
364  }
365 
366  protected:
367 
368  /// Pointer to isotropic growth function
370 
371  /// Pointer to prestress function
373 
374  /// Pointer to the constitutive law
376 
377  /// Timescale ratio (non-dim. density)
378  double* Lambda_sq_pt;
379 
380  /// Flag that switches inertia on/off
381  bool Unsteady;
382 
383  /// Pointer to body force function
385 
386  /// Static default value for timescale ratio (1.0 -- for natural scaling)
387  static double Default_lambda_sq_value;
388 
389  /// Use FD to evaluate Jacobian
391 
392  };
393 
394 
395 //=======================================================================
396 /// A class for elements that solve the equations of solid mechanics, based
397 /// on the principle of virtual displacements in cartesian coordinates.
398 //=======================================================================
399 template <unsigned DIM>
400 class PVDEquations : public virtual PVDEquationsBase<DIM>
401 {
402 
403  public:
404 
405  /// \short Constructor
407 
408  /// \short Return the 2nd Piola Kirchoff stress tensor, as calculated
409  /// from the constitutive law at specified local coordinate
410  void get_stress(const Vector<double> &s, DenseMatrix<double> &sigma);
411 
412  /// \short Fill in the residuals for the solid equations (the discretised
413  /// principle of virtual displacements)
415  {
418  }
419 
420  /// \short Fill in contribution to Jacobian (either by FD or analytically,
421  /// control this via evaluate_jacobian_by_fd()
423  DenseMatrix<double> &jacobian)
424  {
425 
426  //Solve for the consistent acceleration in Newmark scheme?
428  {
429  // Add the contribution to the residuals -- these are the
430  // full residuals of whatever solid equations we're solving
431  this->fill_in_contribution_to_residuals(residuals);
432 
433  // Jacobian is not the Jacobian associated with these
434  // residuals (treating the postions as unknowns)
435  // but the derivatives w.r.t. to the discrete generalised
436  // accelerations in the Newmark scheme -- the Jacobian
437  // is therefore the associated mass matrix, multiplier
438  // by suitable scaling factors
439  this->fill_in_jacobian_for_newmark_accel(jacobian);
440  return;
441  }
442 
443  // Are we assigning a solid initial condition?
444  if (this->Solid_ic_pt!=0)
445  {
446  this->fill_in_jacobian_for_solid_ic(residuals,jacobian);
447  return;
448  }
449 
450 
451  // Use FD
452  if ((this->Evaluate_jacobian_by_fd))
453  {
454  //Add the contribution to the residuals from this element
457 
458  //Get the solid entries in the jacobian using finite differences
460  }
461  // Do it analytically
462  else
463  {
464  fill_in_generic_contribution_to_residuals_pvd(residuals,jacobian,1);
465  }
466  }
467 
468 
469  /// Output: x,y,[z],xi0,xi1,[xi2],gamma
470  void output(std::ostream &outfile)
471  {
472  unsigned n_plot=5;
473  output(outfile,n_plot);
474  }
475 
476  /// Output: x,y,[z],xi0,xi1,[xi2],gamma
477  void output(std::ostream &outfile, const unsigned &n_plot);
478 
479 
480  /// C-style output: x,y,[z],xi0,xi1,[xi2],gamma
481  void output(FILE* file_pt)
482  {
483  unsigned n_plot=5;
484  output(file_pt,n_plot);
485  }
486 
487  /// Output: x,y,[z],xi0,xi1,[xi2],gamma
488  void output(FILE* file_pt, const unsigned &n_plot);
489 
490 
491  /// \short Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress
492  /// components
493  void extended_output(std::ostream &outfile, const unsigned &n_plot);
494 
495 
496  protected:
497 
498  /// \short Compute element residual Vector only (if flag=and/or element
499  /// Jacobian matrix
501  Vector<double> &residuals, DenseMatrix<double> &jacobian,
502  const unsigned& flag);
503 
504  /// \short Return the 2nd Piola Kirchhoff stress tensor, as
505  /// calculated from the constitutive law: Pass metric tensors in the
506  /// stress free and current configurations.
507  inline void get_stress(const DenseMatrix<double> &g,
508  const DenseMatrix<double> &G,
509  DenseMatrix<double> &sigma)
510  {
511 #ifdef PARANOID
512  //If the pointer to the constitutive law hasn't been set, issue an error
513  if(this->Constitutive_law_pt==0)
514  {
515  //Write an error message
516  std::string error_message =
517  "Elements derived from PVDEquations must have a constitutive law:\n";
518  error_message +=
519  "set one using the constitutive_law_pt() member function";
520  //Throw the error
521  throw OomphLibError(error_message,
522  OOMPH_CURRENT_FUNCTION,
523  OOMPH_EXCEPTION_LOCATION);
524  }
525 #endif
526  this->Constitutive_law_pt
528  }
529 
530  /// \short Return the derivatives of the 2nd Piola Kirchhoff stress tensor,
531  /// as calculated from the constitutive law: Pass metric tensors in the
532  /// stress free and current configurations and the current value of the
533  /// the stress tensor.
535  const DenseMatrix<double> &G,
536  const DenseMatrix<double> &sigma,
537  RankFourTensor<double> &d_sigma_dG)
538  {
539 #ifdef PARANOID
540  //If the pointer to the constitutive law hasn't been set, issue an error
541  if(this->Constitutive_law_pt==0)
542  {
543  //Write an error message
544  std::string error_message =
545  "Elements derived from PVDEquations must have a constitutive law:\n";
546  error_message +=
547  "set one using the constitutive_law_pt() member function";
548  //Throw the error
549  throw OomphLibError(error_message,
550  OOMPH_CURRENT_FUNCTION,
551  OOMPH_EXCEPTION_LOCATION);
552  }
553 #endif
554  //Only bother with the symmetric part by passing false as last entry
555  this->Constitutive_law_pt
556  ->calculate_d_second_piola_kirchhoff_stress_dG(g,G,sigma,d_sigma_dG,
557  false);
558  }
559 
560 
561 
562  private:
563 
564  /// Unpin all solid pressure dofs -- empty as there are no pressures
566 
567  };
568 
569 
570 //===========================================================================
571 /// An Element that solves the solid mechanics equations, based on
572 /// the principle of virtual displacements in Cartesian coordinates,
573 /// using SolidQElements for the interpolation of the variable positions.
574 //============================================================================
575  template<unsigned DIM, unsigned NNODE_1D>
576  class QPVDElement : public virtual SolidQElement<DIM,NNODE_1D>,
577  public virtual PVDEquations<DIM>
578  {
579  public:
580 
581  /// Constructor, there are no internal data points
582  QPVDElement() : SolidQElement<DIM,NNODE_1D>(), PVDEquations<DIM>() { }
583 
584  /// Output function
585  void output(std::ostream &outfile) {PVDEquations<DIM>::output(outfile);}
586 
587  /// Output function
588  void output(std::ostream &outfile, const unsigned &n_plot)
589  {PVDEquations<DIM>::output(outfile,n_plot);}
590 
591 
592  /// C-style output function
593  void output(FILE* file_pt) {PVDEquations<DIM>::output(file_pt);}
594 
595  /// C-style output function
596  void output(FILE* file_pt, const unsigned &n_plot)
597  {PVDEquations<DIM>::output(file_pt,n_plot);}
598 
599  };
600 
601 
602 //============================================================================
603 /// FaceGeometry of a 2D QPVDElement element
604 //============================================================================
605  template<unsigned NNODE_1D>
606  class FaceGeometry<QPVDElement<2,NNODE_1D> > :
607  public virtual SolidQElement<1,NNODE_1D>
608  {
609  public:
610  /// Constructor must call the constructor of the underlying solid element
611  FaceGeometry() : SolidQElement<1,NNODE_1D>() {}
612  };
613 
614 
615  //==============================================================
616 /// FaceGeometry of the FaceGeometry of the 2D QPVDElement
617 //==============================================================
618 template<unsigned NNODE_1D>
619 class FaceGeometry<FaceGeometry<QPVDElement<2,NNODE_1D> > >:
620  public virtual PointElement
621 {
622  public:
623  //Make sure that we call the constructor of the SolidQElement
624  //Only the Intel compiler seems to need this!
626 };
627 
628 
629 //============================================================================
630 /// FaceGeometry of a 3D QPVDElement element
631 //============================================================================
632 template<unsigned NNODE_1D>
633  class FaceGeometry<QPVDElement<3,NNODE_1D> > :
634  public virtual SolidQElement<2,NNODE_1D>
635  {
636  public:
637  /// Constructor must call the constructor of the underlying solid element
638  FaceGeometry() : SolidQElement<2,NNODE_1D>() {}
639  };
640 
641 //============================================================================
642 /// FaceGeometry of FaceGeometry of a 3D QPVDElement element
643 //============================================================================
644  template<unsigned NNODE_1D>
645  class FaceGeometry<FaceGeometry<QPVDElement<3,NNODE_1D> > > :
646  public virtual SolidQElement<1,NNODE_1D>
647  {
648  public:
649  /// Constructor must call the constructor of the underlying solid element
650  FaceGeometry() : SolidQElement<1,NNODE_1D>() {}
651  };
652 
653 //===========================================================================
654 /// An Element that solves the principle of virtual diplacements
655 /// using Hermite interpolation for the variable positions.
656 //============================================================================
657  template<unsigned DIM>
658  class HermitePVDElement : public virtual SolidQHermiteElement<DIM>,
659  public virtual PVDEquations<DIM>
660  {
661 
662  public:
663 
664  /// Constructor, there are no internal data points
666  PVDEquations<DIM>() { }
667 
668  /// SolidQHermiteElement output function
669  void output(std::ostream &outfile)
671 
672  /// SolidQHermiteElement output function
673  void output(std::ostream &outfile, const unsigned &n_plot)
674  {SolidQHermiteElement<DIM>::output(outfile,n_plot);}
675 
676  /// C-style SolidQHermiteElement output function
677  void output(FILE* file_pt) {SolidQHermiteElement<DIM>::output(file_pt);}
678 
679  /// C-style SolidQHermiteElement output function
680  void output(FILE* file_pt, const unsigned &n_plot)
681  {SolidQHermiteElement<DIM>::output(file_pt,n_plot);}
682 
683  };
684 
685 
686 
687 
688 //==========================================================
689 /// PVDElementWithContinuousPressure upgraded to become projectable
690 //==========================================================
691  template<class PVD_ELEMENT>
693  public virtual ProjectableElement<PVD_ELEMENT>
694  {
695 
696  public:
697 
698  /// \short Constructor [this was only required explicitly
699  /// from gcc 4.5.2 onwards...]
701 
702 
703  /// \short Specify the values associated with field fld.
704  /// The information is returned in a vector of pairs which comprise
705  /// the Data object and the value within it, that correspond to field fld.
706  /// In the underlying PVD elements there are no field values
708  {
709  // Create the vector
710  Vector<std::pair<Data*,unsigned> > data_values;
711 
712  // Loop over all vertex nodes
713  //const unsigned n_solid_pres=this->npres_solid();
714  //for(unsigned j=0;j<n_solid_pres;j++)
715  // {
716  // // Add the data value associated with the pressure components
717  // unsigned vertex_index=this->Pconv[j];
718  // data_values.push_back(std::make_pair(this->node_pt(vertex_index),0));
719  // }
720 
721  // Return the vector
722  return data_values;
723  }
724 
725  /// \short Number of fields to be projected: 0
726  unsigned nfields_for_projection() {return 0;}
727 
728  /// \short Number of history values to be stored for fld-th field
729  /// (Includes the current value!). No nodal data.
730  unsigned nhistory_values_for_projection(const unsigned &fld)
731  {return 0;}
732 
733  ///\short Number of positional history values (Includes the current value!)
735  {
736  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
737  }
738 
739  /// \short Return Jacobian of mapping and shape functions of field fld
740  /// at local coordinate s
741  double jacobian_and_shape_of_field(const unsigned &fld,
742  const Vector<double> &s,
743  Shape &psi)
744  {
745  //Return the Jacobian of the eulerian mapping
746  return this->J_eulerian(s);
747  }
748 
749  /// \short Return interpolated field fld at local coordinate s, at time level
750  /// t (t=0: present; t>0: history values)
751  double get_field(const unsigned &t,
752  const unsigned &fld,
753  const Vector<double>& s)
754  {
755  //Dummy return
756  return 0.0;
757  }
758 
759 
760  ///Return number of values in field fld
761  unsigned nvalue_of_field(const unsigned &fld)
762  {
763  return 0;
764  }
765 
766 
767  ///Return local equation number of value j in field fld.
768  int local_equation(const unsigned &fld,
769  const unsigned &j)
770  {
771  return -1;
772  }
773 
774  };
775 
776 
777 //=======================================================================
778 /// Face geometry for element is the same as that for the underlying
779 /// wrapped element
780 //=======================================================================
781  template<class ELEMENT>
783  : public virtual FaceGeometry<ELEMENT>
784  {
785  public:
786  FaceGeometry() : FaceGeometry<ELEMENT>() {}
787  };
788 
789 
790 //=======================================================================
791 /// Face geometry of the Face Geometry for element is the same as
792 /// that for the underlying wrapped element
793 //=======================================================================
794  template<class ELEMENT>
796  ProjectablePVDElement<ELEMENT> > >
797  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
798  {
799  public:
801  };
802 
803 
804 
805 ///////////////////////////////////////////////////////////////////////
806 ///////////////////////////////////////////////////////////////////////
807 ///////////////////////////////////////////////////////////////////////
808 
809 
810 
811 //=========================================================================
812 /// A class for elements that solve the equations of solid mechanics,
813 /// based on the principle of virtual displacements, with a
814 /// contitutive equation that involves a pressure. This
815 /// formulation is required in the case of incompressible materials, in
816 /// which the additional constraint that volume must be conserved is applied.
817 /// In this case, the Incompressible flag must be set to true. If the
818 /// Incompressible flag is not set to true, we use the nearly-incompressible
819 /// formulation of the constitutive equations.
820 //============================================================================
821  template <unsigned DIM>
822  class PVDEquationsWithPressure : public virtual PVDEquationsBase<DIM>,
824  {
825  public:
826 
827  /// Constructor, by default the element is NOT incompressible.
829  Incompressible(false) {}
830 
831  /// \short Return the 2nd Piola Kirchoff stress tensor, as calculated
832  /// from the constitutive law at specified local coordinate
833  void get_stress(const Vector<double> &s, DenseMatrix<double> &sigma);
834 
835  /// Return whether the material is incompressible
836  bool is_incompressible() const {return Incompressible;}
837 
838  /// Set the material to be incompressible
840 
841  /// Set the material to be compressible
843 
844  /// Return the lth solid pressure
845  virtual double solid_p(const unsigned &l)=0;
846 
847  /// Set the lth solid pressure to p_value
848  virtual void set_solid_p(const unsigned &l, const double &p_value)=0;
849 
850  /// Fill in the residuals
852  {
853  //Call the generic residuals function with flag set to 0
854  //using a dummy matrix argument
858  }
859 
860  /// \short Fill in contribution to Jacobian (either by FD or analytically,
861  /// for the positional variables; control this via
862  /// evaluate_jacobian_by_fd(). Note: Jacobian entries arising from
863  /// derivatives w.r.t. pressure terms are always computed analytically.
865  DenseMatrix<double> &jacobian)
866  {
867 
868  //Solve for the consistent acceleration in the Newmark scheme
869  //Note that this replaces solid entries only
871  (this->Solid_ic_pt!=0))
872  {
873  std::string error_message ="Can't assign consistent Newmark history\n";
874  error_message += " values for solid element with pressure dofs\n";
875 
876  throw OomphLibError(
877  error_message,
878  OOMPH_CURRENT_FUNCTION,
879  OOMPH_EXCEPTION_LOCATION);
880  }
881 
882  // FD
883  if (this->Evaluate_jacobian_by_fd)
884  {
885  // Call the generic routine with the flag set to 2: Computes residuals
886  // and derivatives w.r.t. to pressure variables
888  residuals,jacobian,GeneralisedElement::Dummy_matrix,2);
889 
890  // Call the finite difference routine for the deriatives w.r.t.
891  // the positional variables
893 
894  }
895  // Do it fully analytically
896  else
897  {
898  //Call the generic routine with the flag set to 1: Get residual
899  // and fully analytical Jacobian
901  residuals,jacobian,GeneralisedElement::Dummy_matrix,1);
902  }
903  }
904 
905  /// \short Fill in contribution to Mass matrix and
906  /// Jacobian (either by FD or analytically,
907  /// for the positional variables; control this via
908  /// evaluate_jacobian_by_fd(). Note: Jacobian entries arising from
909  /// derivatives w.r.t. pressure terms are always computed analytically.
910  /// Note that the Jacobian is multiplied by minus one to
911  /// ensure that the mass matrix is positive semi-definite.
913  Vector<double> &residuals,
914  DenseMatrix<double> &jacobian,
915  DenseMatrix<double> &mass_matrix)
916  {
917 
918  //Solve for the consistent acceleration in the Newmark scheme
919  //Note that this replaces solid entries only
921  (this->Solid_ic_pt!=0))
922  {
923  std::string error_message ="Can't assign consistent Newmark history\n";
924  error_message += " values for solid element with pressure dofs\n";
925 
926  throw OomphLibError(
927  error_message,
928  OOMPH_CURRENT_FUNCTION,
929  OOMPH_EXCEPTION_LOCATION);
930  }
931 
932  // FD
933  if (this->Evaluate_jacobian_by_fd)
934  {
935  // Call the generic routine with the flag set to 4: Computes residuals
936  // and derivatives w.r.t. to pressure variables
937  // and the mass matrix
939  residuals,jacobian,mass_matrix,4);
940 
941  // Call the finite difference routine for the deriatives w.r.t.
942  // the positional variables
944  }
945  // Do it fully analytically
946  else
947  {
948  //Call the generic routine with the flag set to 3: Get residual
949  // and fully analytical Jacobian
951  residuals,jacobian,mass_matrix,3);
952  }
953 
954  //Multiply the residuals and jacobian by minus one
955  const unsigned n_dof = this->ndof();
956  for(unsigned i=0;i<n_dof;i++)
957  {
958  residuals[i] *= -1.0;
959  for(unsigned j=0;j<n_dof;j++)
960  {
961  jacobian(i,j) *= -1.0;
962  }
963  }
964  }
965 
966 
967 
968  /// Return the interpolated_solid_pressure
970  {
971  //Find number of nodes
972  unsigned n_solid_pres = this->npres_solid();
973  //Local shape function
974  Shape psisp(n_solid_pres);
975  //Find values of shape function
976  solid_pshape(s,psisp);
977 
978  //Initialise value of solid_p
979  double interpolated_solid_p = 0.0;
980  //Loop over the local nodes and sum
981  for(unsigned l=0;l<n_solid_pres;l++)
982  {interpolated_solid_p += solid_p(l)*psisp[l];}
983 
984  return(interpolated_solid_p);
985  }
986 
987 
988  /// Output: x,y,[z],xi0,xi1,[xi2],p,gamma
989  void output(std::ostream &outfile)
990  {
991  unsigned n_plot=5;
992  output(outfile,n_plot);
993  }
994 
995  /// Output: x,y,[z],xi0,xi1,[xi2],p,gamma
996  void output(std::ostream &outfile, const unsigned &n_plot);
997 
998 
999  /// C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma
1000  void output(FILE* file_pt)
1001  {
1002  unsigned n_plot=5;
1003  output(file_pt,n_plot);
1004  }
1005 
1006  /// C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma
1007  void output(FILE* file_pt, const unsigned &n_plot);
1008 
1009  /// \short Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress
1010  /// components
1011  void extended_output(std::ostream &outfile, const unsigned &n_plot);
1012 
1013 
1014  /// \short Compute the diagonal of the displacement mass matrix for
1015  /// LSC preconditioner
1016  void get_mass_matrix_diagonal(Vector<double> &mass_diag);
1017 
1018  /// \short returns the number of DOF types associated with this element:
1019  /// displacement components and pressure
1020  unsigned ndof_types() const
1021  {
1022  return DIM+1;
1023  }
1024 
1025  /// \short Create a list of pairs for all unknowns in this element,
1026  /// so that the first entry in each pair contains the global equation
1027  /// number of the unknown, while the second one contains the number
1028  /// of the "DOF" that this unknown is associated with.
1029  /// (Function can obviously only be called if the equation numbering
1030  /// scheme has been set up.)
1031  /// There are DIM+1 types of DOF: displacement compnents and
1032  /// pressure
1034  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
1035  {
1036  // temporary pair (used to store dof lookup prior to being added to list
1037  std::pair<unsigned,unsigned> dof_lookup;
1038 
1039  // number of nodes
1040  const unsigned n_node = this->nnode();
1041 
1042  //Get the number of position dofs and dimensions at the node
1043  const unsigned n_position_type = this->nnodal_position_type();
1044  const unsigned nodal_dim = this->nodal_dimension();
1045 
1046  //Integer storage for local unknown
1047  int local_unknown=0;
1048 
1049  //Loop over the nodes
1050  for(unsigned n=0;n<n_node;n++)
1051  {
1052 
1053  //Loop over position dofs
1054  for(unsigned k=0;k<n_position_type;k++)
1055  {
1056  //Loop over dimension
1057  for(unsigned i=0;i<nodal_dim;i++)
1058  {
1059  //If the variable is free
1060  local_unknown = this->position_local_eqn(n,k,i);
1061 
1062  // ignore pinned values
1063  if (local_unknown >= 0)
1064  {
1065  // store dof lookup in temporary pair: First entry in pair
1066  // is global equation number; second entry is dof type
1067  dof_lookup.first = this->eqn_number(local_unknown);
1068  dof_lookup.second = i;
1069 
1070  // add to list
1071  dof_lookup_list.push_front(dof_lookup);
1072  }
1073  }
1074  }
1075  }
1076 
1077  //Do solid pressure degrees of freedom
1078  unsigned np=this->npres_solid();
1079  for (unsigned j=0;j<np;j++)
1080  {
1081  int local_unknown=this->solid_p_local_eqn(j);
1082  // ignore pinned values
1083  if (local_unknown >= 0)
1084  {
1085  // store dof lookup in temporary pair: First entry in pair
1086  // is global equation number; second entry is dof type
1087  dof_lookup.first = this->eqn_number(local_unknown);
1088  dof_lookup.second = DIM;
1089 
1090  // add to list
1091  dof_lookup_list.push_front(dof_lookup);
1092  }
1093  }
1094 
1095  }
1096 
1097  protected:
1098 
1099  /// \short Return the deviatoric part of the 2nd Piola Kirchhoff stress
1100  /// tensor, as calculated from the constitutive law in the nearly
1101  /// incompresible formulation. Also return the contravariant
1102  /// deformed metric tensor, the generalised dilatation, and the
1103  /// inverse of the bulk modulus.
1104  inline void get_stress(const DenseMatrix<double> &g,
1105  const DenseMatrix<double> &G,
1106  DenseMatrix<double> &sigma_dev,
1107  DenseMatrix<double> &Gcontra,
1108  double &gen_dil, double &inv_kappa)
1109  {
1110 #ifdef PARANOID
1111  //If the pointer to the constitutive law hasn't been set, issue an error
1112  if(this->Constitutive_law_pt == 0)
1113  {
1114  //Write an error message
1115  std::string error_message =
1116  "Elements derived from PVDEquationsWithPressure \n";
1117  error_message += "must have a constitutive law:\n";
1118  error_message +=
1119  "set one using the constitutive_law_pt() member function";
1120  //Throw the error
1121  throw OomphLibError(error_message,
1122  OOMPH_CURRENT_FUNCTION,
1123  OOMPH_EXCEPTION_LOCATION);
1124  }
1125 #endif
1126  this->Constitutive_law_pt->
1127  calculate_second_piola_kirchhoff_stress(g,G,sigma_dev,Gcontra,
1128  gen_dil,inv_kappa);
1129  }
1130 
1131 
1132  /// \short Return the derivative of the
1133  /// deviatoric part of the 2nd Piola Kirchhoff stress
1134  /// tensor, as calculated from the constitutive law in the nearly
1135  /// incompresible formulation. Also return the derivative of the
1136  /// generalised dilatation.
1138  const DenseMatrix<double> &G,
1139  const DenseMatrix<double> &sigma,
1140  const double &gen_dil,
1141  const double &inv_kappa,
1142  const double &interpolated_solid_p,
1143  RankFourTensor<double> &d_sigma_dG,
1144  DenseMatrix<double> &d_gen_dil_dG)
1145 
1146  {
1147 #ifdef PARANOID
1148  //If the pointer to the constitutive law hasn't been set, issue an error
1149  if(this->Constitutive_law_pt == 0)
1150  {
1151  //Write an error message
1152  std::string error_message =
1153  "Elements derived from PVDEquationsWithPressure \n";
1154  error_message += "must have a constitutive law:\n";
1155  error_message +=
1156  "set one using the constitutive_law_pt() member function";
1157  //Throw the error
1158  throw OomphLibError(error_message,
1159  OOMPH_CURRENT_FUNCTION,
1160  OOMPH_EXCEPTION_LOCATION);
1161  }
1162 #endif
1163  //Only bother with the symmetric part by passing false as last entry
1164  this->Constitutive_law_pt->
1165  calculate_d_second_piola_kirchhoff_stress_dG(
1166  g,G,sigma,gen_dil,inv_kappa,interpolated_solid_p,
1167  d_sigma_dG,d_gen_dil_dG,false);
1168  }
1169 
1170 
1171  /// Return the solid pressure shape functions
1172  virtual void solid_pshape(const Vector<double> &s, Shape &psi) const=0;
1173 
1174  /// Return the stored solid shape functions at the knots
1175  void solid_pshape_at_knot(const unsigned &ipt, Shape &psi) const
1176  {
1177  //Find the dimension of the element
1178  unsigned Dim = this->dim();
1179  //Storage for local coordinates of the integration point
1180  Vector<double> s(Dim);
1181  //Set the local coordinates
1182  for(unsigned i=0;i<Dim;i++) {s[i] = this->integral_pt()->knot(ipt,i);}
1183  //Get the shape function
1184  solid_pshape(s,psi);
1185  }
1186 
1187  protected:
1188 
1189  /// Boolean to determine whether the solid is incompressible or not
1191 
1192  /// \short Returns the residuals for the discretised principle of
1193  /// virtual displacements, formulated in the incompressible/
1194  /// near-incompressible case.
1195  /// - If flag==0, compute only the residual vector.
1196  /// - If flag==1, compute residual vector and fully analytical Jacobian
1197  /// - If flag==2, also compute the pressure-related entries
1198  /// in the Jacobian (all others need to be done by finite differencing.
1200  Vector<double> &residuals, DenseMatrix<double> &jacobian,
1201  DenseMatrix<double> &mass_matrix,
1202  const unsigned& flag);
1203 
1204  /// \short Return the deviatoric part of the 2nd Piola Kirchhoff stress
1205  /// tensor, as calculated from the constitutive law in the
1206  /// incompresible formulation. Also return the contravariant
1207  /// deformed metric tensor, and the
1208  /// determinant of the deformed covariant metric tensor
1209  /// (likely to be needed in the incompressibility constraint)
1210  inline void get_stress(const DenseMatrix<double> &g,
1211  const DenseMatrix<double> &G,
1212  DenseMatrix<double> &sigma_dev,
1213  DenseMatrix<double> &Gcontra,
1214  double &detG)
1215  {
1216 #ifdef PARANOID
1217  //If the pointer to the constitutive law hasn't been set, issue an error
1218  if(this->Constitutive_law_pt == 0)
1219  {
1220  //Write an error message
1221  std::string error_message =
1222  "Elements derived from PVDEquationsWithPressure \n";
1223  error_message += "must have a constitutive law:\n";
1224  error_message +=
1225  "set one using the constitutive_law_pt() member function";
1226  //Throw the error
1227  throw OomphLibError(error_message,
1228  OOMPH_CURRENT_FUNCTION,
1229  OOMPH_EXCEPTION_LOCATION);
1230  }
1231 #endif
1232  this->Constitutive_law_pt->
1233  calculate_second_piola_kirchhoff_stress(g,G,sigma_dev,Gcontra,detG);
1234  }
1235 
1236  /// \short Return the derivative of the 2nd Piola Kirchhoff stress
1237  /// tensor, as calculated from the constitutive law in the
1238  /// incompresible formulation. Also return
1239  /// derivative of the determinant of the deformed covariant metric tensor
1240  /// (likely to be needed in the incompressibility constraint)
1242  const DenseMatrix<double> &G,
1243  const DenseMatrix<double> &sigma,
1244  const double &detG,
1245  const double &interpolated_solid_p,
1246  RankFourTensor<double> &d_sigma_dG,
1247  DenseMatrix<double> &d_detG_dG)
1248  {
1249 #ifdef PARANOID
1250  //If the pointer to the constitutive law hasn't been set, issue an error
1251  if(this->Constitutive_law_pt == 0)
1252  {
1253  //Write an error message
1254  std::string error_message =
1255  "Elements derived from PVDEquationsWithPressure \n";
1256  error_message += "must have a constitutive law:\n";
1257  error_message +=
1258  "set one using the constitutive_law_pt() member function";
1259  //Throw the error
1260  throw OomphLibError(error_message,
1261  OOMPH_CURRENT_FUNCTION,
1262  OOMPH_EXCEPTION_LOCATION);
1263  }
1264 #endif
1265  //Only bother with the symmetric part by passing false as last entry
1266  this->Constitutive_law_pt->
1267  calculate_d_second_piola_kirchhoff_stress_dG(
1268  g,G,sigma,detG,interpolated_solid_p,d_sigma_dG,d_detG_dG,false);
1269  }
1270 
1271 
1272  };
1273 
1274 
1275 
1276 //////////////////////////////////////////////////////////////////////////////
1277 //////////////////////////////////////////////////////////////////////////////
1278 //////////////////////////////////////////////////////////////////////////////
1279 
1280 
1281 
1282 
1283 //===========================================================================
1284 /// An Element that solves the equations of solid mechanics, using the
1285 /// principle of virtual displacements, with quadratic interpolation
1286 /// for the positions and a discontinuous linear solid pressure. This is
1287 /// analogous to the QCrouzeixRaviartElement element for fluids.
1288 //============================================================================
1289 template<unsigned DIM>
1290 class QPVDElementWithPressure : public virtual SolidQElement<DIM,3>,
1291  public virtual PVDEquationsWithPressure<DIM>
1292 {
1293 
1294  /// \short Unpin all solid pressure dofs in the element
1296  {
1297  unsigned n_pres = this->npres_solid();
1298  // loop over pressure dofs and unpin them
1299  for(unsigned l=0;l<n_pres;l++)
1300  {this->internal_data_pt(this->P_solid_internal_index)->unpin(l);}
1301  }
1302 
1303  protected:
1304 
1305  /// \short Internal index that indicates at which internal data value the
1306  /// solid presure is stored
1308 
1309  /// \short Overload the access function
1310  /// that is used to return local equation corresponding to the i-th
1311  /// solid pressure value
1312  inline int solid_p_local_eqn(const unsigned &i) const
1313  {return this->internal_local_eqn(P_solid_internal_index,i);}
1314 
1315  /// Return the pressure shape functions
1316  inline void solid_pshape(const Vector<double> &s, Shape &psi) const;
1317 
1318  public:
1319 
1320  /// \short There is internal solid data so we can't use the automatic
1321  /// assignment of consistent initial conditions for time-dependent problems.
1322  bool has_internal_solid_data() {return true;}
1323 
1324  /// Constructor, there are DIM+1 internal data points
1326  PVDEquationsWithPressure<DIM>()
1327  {
1328  //Allocate and add one Internal data object that stores DIM+1 pressure
1329  //values
1330  P_solid_internal_index = this->add_internal_data(new Data(DIM+1));
1331  }
1332 
1333  /// Return the lth pressure value
1334  double solid_p(const unsigned &l)
1335  {return this->internal_data_pt(P_solid_internal_index)->value(l);}
1336 
1337  /// Set the l-th pressure value to p_value
1338  void set_solid_p(const unsigned &l, const double &p_value)
1339  {this->internal_data_pt(P_solid_internal_index)->set_value(l,p_value);}
1340 
1341  /// Return number of pressure values
1342  unsigned npres_solid() const {return DIM+1;}
1343 
1344  /// Fix the pressure dof l to be the value pvalue
1345  void fix_solid_pressure(const unsigned &l, const double &pvalue)
1346  {
1349  }
1350 
1351  /// Generic FiniteElement output function
1352  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
1353 
1354  /// PVDEquationsWithPressure output function
1355  void output(std::ostream &outfile, const unsigned &n_plot)
1356  {PVDEquationsWithPressure<DIM>::output(outfile,n_plot);}
1357 
1358 
1359  /// C-style Generic FiniteElement output function
1360  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
1361 
1362  /// C-style PVDEquationsWithPressure output function
1363  void output(FILE* file_pt, const unsigned &n_plot)
1364  {PVDEquationsWithPressure<DIM>::output(file_pt,n_plot);}
1365 
1366 };
1367 
1368 
1369 //=====================================================================
1370 /// Pressure shape functions for 2D QPVDElementWithPressure elements
1371 //=====================================================================
1372 template<>
1374  Shape &psi) const
1375 {
1376  psi[0] = 1.0;
1377  psi[1] = s[0];
1378  psi[2] = s[1];
1379 }
1380 
1381 
1382 //=====================================================================
1383 /// Pressure shape functions for 3D QPVDElementWithPressure elements
1384 //=====================================================================
1385 template<>
1387  Shape &psi) const
1388 {
1389  psi[0] = 1.0;
1390  psi[1] = s[0];
1391  psi[2] = s[1];
1392  psi[3] = s[2];
1393 }
1394 
1395 
1396 
1397 //======================================================================
1398 /// FaceGeometry of 2D QPVDElementWithPressure
1399 //======================================================================
1400 template<>
1402 public virtual SolidQElement<1,3>
1403 {
1404  public:
1405  /// Constructor must call constructor of underlying solid element
1407 };
1408 
1409 
1410 //======================================================================
1411 /// FaceGeometry of FaceGeometry of 2D QPVDElementWithPressure
1412 //======================================================================
1413 template<>
1415 public virtual PointElement
1416 {
1417  public:
1418  /// Constructor must call constructor of underlying solid element
1420 };
1421 
1422 //======================================================================
1423 /// FaceGeometry of 3D QPVDElementWithPressure
1424 //======================================================================
1425 template<>
1427 public virtual SolidQElement<2,3>
1428 {
1429  public:
1430  /// Constructor must call constructor of underlying solid element
1432 };
1433 
1434 
1435 //======================================================================
1436 /// FaceGeometry of FaceGeometry of 3D QPVDElementWithPressure
1437 //======================================================================
1438 template<>
1440 public virtual SolidQElement<1,3>
1441 {
1442  public:
1443  /// Constructor must call constructor of underlying solid element
1445 };
1446 
1447 
1448 ///////////////////////////////////////////////////////////////////////////
1449 ///////////////////////////////////////////////////////////////////////////
1450 ///////////////////////////////////////////////////////////////////////////
1451 
1452 
1453 
1454 
1455 //===========================================================================
1456 /// An Element that solves the equations of solid mechanics, based on
1457 /// the discretised principle of virtual displacements, using quadratic
1458 /// interpolation
1459 /// for the positions and continuous linear solid pressure. This is analagous
1460 /// to the QTaylorHoodElement fluid element.
1461 //============================================================================
1462 template<unsigned DIM>
1463 class QPVDElementWithContinuousPressure : public virtual SolidQElement<DIM,3>,
1464  public virtual PVDEquationsWithPressure<DIM>
1465 {
1466  private:
1467 
1468  /// Static array of ints to hold number of solid pressure values at each node
1469  static const unsigned Initial_Nvalue[];
1470 
1471  /// \short Unpin all solid pressure dofs in the element
1473  {
1474  //find the index at which the pressure is stored
1475  int p_index = this->solid_p_nodal_index();
1476  unsigned n_node = this->nnode();
1477  // loop over nodes
1478  for(unsigned n=0;n<n_node;n++)
1479  {this->node_pt(n)->unpin(p_index);}
1480  }
1481 
1482 protected:
1483 
1484  /// \short Static array of ints to hold conversion from pressure node
1485  /// numbers to actual node numbers
1486  static const unsigned Pconv[];
1487 
1488  /// \short Overload the access function
1489  /// that is used to return local equation corresponding to the i-th
1490  /// solid pressure value
1491  inline int solid_p_local_eqn(const unsigned &i) const
1492  {return this->nodal_local_eqn(Pconv[i],this->solid_p_nodal_index());}
1493 
1494  /// Return the pressure shape functions
1495  inline void solid_pshape(const Vector<double> &s, Shape &psi) const;
1496 
1497 public:
1498 
1499  /// Constructor
1501  PVDEquationsWithPressure<DIM>()
1502  { }
1503 
1504  /// \short Set the value at which the solid pressure is stored in the nodes
1505  inline int solid_p_nodal_index() const {return 0;}
1506 
1507  /// \short Number of values (pinned or dofs) required at node n. Can
1508  /// be overwritten for hanging node version
1509  inline virtual unsigned required_nvalue(const unsigned &n) const
1510  {return Initial_Nvalue[n];}
1511 
1512  /// Return the l-th pressure value, make sure to use the hanging
1513  /// representation if there is one!
1514  double solid_p(const unsigned &l)
1515  {return this->nodal_value(Pconv[l],this->solid_p_nodal_index());}
1516 
1517  /// Set the l-th solid pressure value to p_value
1518  void set_solid_p(const unsigned &l, const double &p_value)
1519  {this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(),p_value);}
1520 
1521  /// Return number of pressure values
1522  unsigned npres_solid() const
1523  {return static_cast<unsigned>(pow(2.0,static_cast<int>(DIM)));}
1524 
1525  /// Fix the pressure dof l to be the value pvalue
1526  void fix_solid_pressure(const unsigned &l, const double &pvalue)
1527  {
1528  this->node_pt(Pconv[l])->pin(this->solid_p_nodal_index());
1529  this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(),pvalue);
1530  }
1531 
1532  /// Generic FiniteElement output function
1533  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
1534 
1535  /// PVDEquationsWithPressure output function
1536  void output(std::ostream &outfile, const unsigned &n_plot)
1537  {PVDEquationsWithPressure<DIM>::output(outfile,n_plot);}
1538 
1539 
1540  /// C-style generic FiniteElement output function
1541  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
1542 
1543  /// C-style PVDEquationsWithPressure output function
1544  void output(FILE* file_pt, const unsigned &n_plot)
1545  {PVDEquationsWithPressure<DIM>::output(file_pt,n_plot);}
1546 
1547 };
1548 
1549 
1550 //===============================================================
1551 /// Pressure shape functions for 2D QPVDElementWithContinuousPressure
1552 /// elements
1553 //===============================================================
1554 template<>
1556  const Vector<double> &s, Shape &psi) const
1557 {
1558  //Local storage
1559  double psi1[2], psi2[2];
1560  //Call the OneDimensional Shape functions
1561  OneDimLagrange::shape<2>(s[0],psi1);
1562  OneDimLagrange::shape<2>(s[1],psi2);
1563 
1564  //Now let's loop over the nodal points in the element
1565  //s1 is the "x" coordinate, s2 the "y"
1566  for(unsigned i=0;i<2;i++)
1567  {
1568  for(unsigned j=0;j<2;j++)
1569  {
1570  /*Multiply the two 1D functions together to get the 2D function*/
1571  psi[2*i + j] = psi2[i]*psi1[j];
1572  }
1573  }
1574 }
1575 
1576 //===============================================================
1577 /// Pressure shape functions for 3D QPVDElementWithContinuousPressure
1578 /// elements
1579 //===============================================================
1580 template<>
1582  const Vector<double> &s, Shape &psi) const
1583 {
1584  //Local storage
1585  double psi1[2], psi2[2], psi3[2];
1586  //Call the OneDimensional Shape functions
1587  OneDimLagrange::shape<2>(s[0],psi1);
1588  OneDimLagrange::shape<2>(s[1],psi2);
1589  OneDimLagrange::shape<2>(s[2],psi3);
1590 
1591  //Now let's loop over the nodal points in the element
1592  //s1 is the "x" coordinate, s2 the "y"
1593  for(unsigned i=0;i<2;i++)
1594  {
1595  for(unsigned j=0;j<2;j++)
1596  {
1597  for(unsigned k=0;k<2;k++)
1598  {
1599  /*Multiply the two 1D functions together to get the 3D function*/
1600  psi[4*i + 2*j + k] = psi3[i]*psi2[j]*psi1[k];
1601  }
1602  }
1603  }
1604 }
1605 
1606 
1607 
1608 
1609 //===============================================================
1610 /// FaceGeometry for 2D QPVDElementWithContinuousPressure element
1611 //===============================================================
1612 template<>
1614 public virtual SolidQElement<1,3>
1615 {
1616  public:
1617  /// Constructor must call constructor of the underlying Solid element
1619 };
1620 
1621 
1622 
1623 //===============================================================
1624 /// FaceGeometry of FaceGeometry
1625 /// for 2D QPVDElementWithContinuousPressure element
1626 //===============================================================
1627 template<>
1629 public virtual PointElement
1630 {
1631  public:
1632  /// Constructor must call constructor of the underlying Point element
1634 };
1635 
1636 
1637 //===============================================================
1638 /// FaceGeometry for 3D QPVDElementWithContinuousPressure element
1639 //===============================================================
1640 template<>
1642 public virtual SolidQElement<2,3>
1643 {
1644  public:
1645  /// Constructor must call constructor of the underlying Solid element
1647 };
1648 
1649 
1650 
1651 //===============================================================
1652 /// FaceGeometry of FaceGeometry
1653 /// for 3D QPVDElementWithContinuousPressure element
1654 //===============================================================
1655 template<>
1657 public virtual SolidQElement<1,3>
1658 {
1659  public:
1660  /// Constructor must call constructor of the underlying element
1662 };
1663 
1664 
1665 
1666 
1667 //////////////////////////////////////////////////////////////////////
1668 //////////////////////////////////////////////////////////////////////
1669 //////////////////////////////////////////////////////////////////////
1670 
1671 
1672 
1673 //===========================================================================
1674 /// An Element that solves the solid mechanics equations, based on
1675 /// the principle of virtual displacements in Cartesian coordinates,
1676 /// using SolidTElements for the interpolation of the variable positions.
1677 //============================================================================
1678 template<unsigned DIM, unsigned NNODE_1D>
1679  class TPVDElement : public virtual SolidTElement<DIM,NNODE_1D>,
1680  public virtual PVDEquations<DIM>,
1681  public virtual ElementWithZ2ErrorEstimator
1682 {
1683  public:
1684 
1685  /// Constructor, there are no internal data points
1686  TPVDElement() : SolidTElement<DIM,NNODE_1D>(), PVDEquations<DIM>() { }
1687 
1688  /// Output function
1689  void output(std::ostream &outfile) {PVDEquations<DIM>::output(outfile);}
1690 
1691  /// Output function
1692  void output(std::ostream &outfile, const unsigned &n_plot)
1693  {PVDEquations<DIM>::output(outfile,n_plot);}
1694 
1695 
1696  /// C-style output function
1697  void output(FILE* file_pt) {PVDEquations<DIM>::output(file_pt);}
1698 
1699  /// C-style output function
1700  void output(FILE* file_pt, const unsigned &n_plot)
1701  {PVDEquations<DIM>::output(file_pt,n_plot);}
1702 
1703  /// \short Order of recovery shape functions for Z2 error estimation:
1704  /// Same order as shape functions.
1705  unsigned nrecovery_order() {return (NNODE_1D-1);}
1706 
1707  /// \short Number of vertex nodes in the element
1708  unsigned nvertex_node() const
1710 
1711  /// \short Pointer to the j-th vertex node in the element
1712  Node* vertex_node_pt(const unsigned& j) const
1714 
1715  /// \short Function to describe the local dofs of the element. The ostream
1716  /// specifies the output stream to which the description
1717  /// is written; the string stores the currently
1718  /// assembled output that is ultimately written to the
1719  /// output stream by Data::desribe_dofs(...); it is typically
1720  /// built up incrementally as we descend through the
1721  /// call hierarchy of this function when called from
1722  /// Problem::describe_dofs(...)
1724 
1725  /// Number of 'flux' terms for Z2 error estimation
1727  {
1728  // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
1729  return DIM + DIM*(DIM-1)/2;
1730  }
1731 
1732  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
1733  /// in strain tensor.
1735  {
1736 #ifdef PARANOID
1737  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
1738  if (flux.size()!=num_entries)
1739  {
1740  std::ostringstream error_message;
1741  error_message << "The flux vector has the wrong number of entries, "
1742  << flux.size() << ", whereas it should be "
1743  << num_entries << std::endl;
1744  throw OomphLibError(error_message.str(),
1745  OOMPH_CURRENT_FUNCTION,
1746  OOMPH_EXCEPTION_LOCATION);
1747  }
1748 #endif
1749 
1750  // Get strain matrix
1751  DenseMatrix<double> strain(DIM);
1752  this->get_strain(s,strain);
1753 
1754  // Pack into flux Vector
1755  unsigned icount=0;
1756 
1757  // Start with diagonal terms
1758  for(unsigned i=0;i<DIM;i++)
1759  {
1760  flux[icount]=strain(i,i);
1761  icount++;
1762  }
1763  //Off diagonals row by row
1764  for(unsigned i=0;i<DIM;i++)
1765  {
1766  for(unsigned j=i+1;j<DIM;j++)
1767  {
1768  flux[icount]=strain(i,j);
1769  icount++;
1770  }
1771  }
1772  }
1773 
1774 
1775 };
1776 
1777 
1778 //============================================================================
1779 /// FaceGeometry of a 2D TPVDElement element
1780 //============================================================================
1781 template<unsigned NNODE_1D>
1782 class FaceGeometry<TPVDElement<2,NNODE_1D> > :
1783  public virtual SolidTElement<1,NNODE_1D>
1784 {
1785  public:
1786  /// Constructor must call the constructor of the underlying solid element
1787  FaceGeometry() : SolidTElement<1,NNODE_1D>() {}
1788 };
1789 
1790 
1791 //==============================================================
1792 /// FaceGeometry of the FaceGeometry of the 2D TPVDElement
1793 //==============================================================
1794 template<unsigned NNODE_1D>
1795 class FaceGeometry<FaceGeometry<TPVDElement<2,NNODE_1D> > >:
1796  public virtual PointElement
1797 {
1798  public:
1799  //Make sure that we call the constructor of the SolidQElement
1800  //Only the Intel compiler seems to need this!
1802 };
1803 
1804 
1805 //============================================================================
1806 /// FaceGeometry of a 3D TPVDElement element
1807 //============================================================================
1808 template<unsigned NNODE_1D>
1809 class FaceGeometry<TPVDElement<3,NNODE_1D> > :
1810  public virtual SolidTElement<2,NNODE_1D>
1811 {
1812  public:
1813  /// Constructor must call the constructor of the underlying solid element
1814  FaceGeometry() : SolidTElement<2,NNODE_1D>() {}
1815 };
1816 
1817 //============================================================================
1818 /// FaceGeometry of FaceGeometry of a 3D TPVDElement element
1819 //============================================================================
1820 template<unsigned NNODE_1D>
1821 class FaceGeometry<FaceGeometry<TPVDElement<3,NNODE_1D> > > :
1822  public virtual SolidTElement<1,NNODE_1D>
1823 {
1824  public:
1825  /// Constructor must call the constructor of the underlying solid element
1826  FaceGeometry() : SolidTElement<1,NNODE_1D>() {}
1827 };
1828 
1829 
1830 ////////////////////////////////////////////////////////////////////////
1831 ////////////////////////////////////////////////////////////////////////
1832 ////////////////////////////////////////////////////////////////////////
1833 
1834 //===========================================================================
1835 /// An Element that solves the solid mechanics equations, based on
1836 /// the principle of virtual displacements in Cartesian coordinates,
1837 /// using SolidTBubbleEnrichedElements for the interpolation of the
1838 /// variable positions. These elements are typically required when using
1839 /// pseudo-elasticity to move internal mesh nodes and TCrouzeixRaviartFluid
1840 /// elements.
1841 //============================================================================
1842 template<unsigned DIM, unsigned NNODE_1D>
1844 public virtual SolidTBubbleEnrichedElement<DIM,NNODE_1D>,
1845  public virtual PVDEquations<DIM>,
1846  public virtual ElementWithZ2ErrorEstimator
1847 {
1848  public:
1849 
1850  /// Constructor, there are no internal data points
1852  SolidTBubbleEnrichedElement<DIM,NNODE_1D>(), PVDEquations<DIM>() { }
1853 
1854  /// Output function
1855  void output(std::ostream &outfile) {PVDEquations<DIM>::output(outfile);}
1856 
1857  /// Output function
1858  void output(std::ostream &outfile, const unsigned &n_plot)
1859  {PVDEquations<DIM>::output(outfile,n_plot);}
1860 
1861 
1862  /// C-style output function
1863  void output(FILE* file_pt) {PVDEquations<DIM>::output(file_pt);}
1864 
1865  /// C-style output function
1866  void output(FILE* file_pt, const unsigned &n_plot)
1867  {PVDEquations<DIM>::output(file_pt,n_plot);}
1868 
1869 
1870  /// \short Order of recovery shape functions for Z2 error estimation:
1871  /// Same order as shape functions.
1872  unsigned nrecovery_order() {return (NNODE_1D-1);}
1873 
1874  /// \short Number of vertex nodes in the element
1875  unsigned nvertex_node() const
1877 
1878  /// \short Pointer to the j-th vertex node in the element
1879  Node* vertex_node_pt(const unsigned& j) const
1881 
1882  /// Number of 'flux' terms for Z2 error estimation
1884  {
1885  // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
1886  return DIM + DIM*(DIM-1)/2;
1887  }
1888 
1889  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
1890  /// in strain tensor.
1892  {
1893 #ifdef PARANOID
1894  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
1895  if (flux.size()!=num_entries)
1896  {
1897  std::ostringstream error_message;
1898  error_message << "The flux vector has the wrong number of entries, "
1899  << flux.size() << ", whereas it should be "
1900  << num_entries << std::endl;
1901  throw OomphLibError(error_message.str(),
1902  OOMPH_CURRENT_FUNCTION,
1903  OOMPH_EXCEPTION_LOCATION);
1904  }
1905 #endif
1906 
1907  // Get strain matrix
1908  DenseMatrix<double> strain(DIM);
1909  this->get_strain(s,strain);
1910 
1911  // Pack into flux Vector
1912  unsigned icount=0;
1913 
1914  // Start with diagonal terms
1915  for(unsigned i=0;i<DIM;i++)
1916  {
1917  flux[icount]=strain(i,i);
1918  icount++;
1919  }
1920  //Off diagonals row by row
1921  for(unsigned i=0;i<DIM;i++)
1922  {
1923  for(unsigned j=i+1;j<DIM;j++)
1924  {
1925  flux[icount]=strain(i,j);
1926  icount++;
1927  }
1928  }
1929  }
1930 
1931 
1932 };
1933 
1934 
1935 //============================================================================
1936 /// FaceGeometry of a 2D TPVDBubbleEnrichedElement element
1937 //============================================================================
1938 template<unsigned NNODE_1D>
1940  public virtual SolidTElement<1,NNODE_1D>
1941 {
1942  public:
1943  /// Constructor must call the constructor of the underlying solid element
1944  FaceGeometry() : SolidTElement<1,NNODE_1D>() {}
1945 };
1946 
1947 
1948 //==============================================================
1949 /// FaceGeometry of the FaceGeometry of the 2D TPVDBubbleEnrichedElement
1950 //==============================================================
1951 template<unsigned NNODE_1D>
1953 public virtual PointElement
1954 {
1955  public:
1956  //Make sure that we call the constructor of the SolidQElement
1957  //Only the Intel compiler seems to need this!
1959 };
1960 
1961 
1962 //============================================================================
1963 /// FaceGeometry of a 3D TPVDBubbleEnrichedElement element
1964 //============================================================================
1965 template<unsigned NNODE_1D>
1967  public virtual SolidTBubbleEnrichedElement<2,NNODE_1D>
1968 {
1969  public:
1970  /// Constructor must call the constructor of the underlying solid element
1972 };
1973 
1974 //============================================================================
1975 /// FaceGeometry of FaceGeometry of a 3D TPVDElement element
1976 //============================================================================
1977 template<unsigned NNODE_1D>
1979  public virtual SolidTElement<1,NNODE_1D>
1980 {
1981  public:
1982  /// Constructor must call the constructor of the underlying solid element
1983  FaceGeometry() : SolidTElement<1,NNODE_1D>() {}
1984 };
1985 
1986 
1987 
1988 
1989 //=======================================================================
1990 /// An Element that solves the solid mechanics equations in an
1991 /// (near) incompressible formulation
1992 /// with quadratic interpolation for velocities and positions and
1993 /// continous linear pressure interpolation. This is equivalent to the
1994 /// TTaylorHoodElement element for fluids.
1995 //=======================================================================
1996 template <unsigned DIM>
1997 class TPVDElementWithContinuousPressure : public virtual SolidTElement<DIM,3>,
1998  public virtual PVDEquationsWithPressure<DIM>,
1999  public virtual ElementWithZ2ErrorEstimator
2000 {
2001  private:
2002 
2003  /// Static array of ints to hold number of variables at node
2004  static const unsigned Initial_Nvalue[];
2005 
2006  /// \short Unpin all solid pressure dofs in the element
2008  {
2009  //find the index at which the pressure is stored
2010  int p_index = this->solid_p_nodal_index();
2011  unsigned n_node = this->nnode();
2012  // loop over nodes
2013  for(unsigned n=0;n<n_node;n++)
2014  {this->node_pt(n)->unpin(p_index);}
2015  }
2016 
2017  protected:
2018 
2019  /// \short Static array of ints to hold conversion from pressure
2020  /// node numbers to actual node numbers
2021  static const unsigned Pconv[];
2022 
2023  /// \short Overload the access function
2024  /// that is used to return local equation corresponding to the i-th
2025  /// solid pressure value
2026  inline int solid_p_local_eqn(const unsigned &i) const
2027  {return this->nodal_local_eqn(Pconv[i],this->solid_p_nodal_index());}
2028 
2029  /// Pressure shape functions at local coordinate s
2030  inline void solid_pshape(const Vector<double> &s, Shape &psi) const;
2031 
2032 public:
2033 
2034  /// Constructor
2036  PVDEquationsWithPressure<DIM>()
2037  { }
2038 
2039 
2040  /// Broken copy constructor
2043  {
2044  BrokenCopy::broken_copy("TPVDElementWithContinuousPressu");
2045  }
2046 
2047  /// Broken assignment operator
2048 //Commented out broken assignment operator because this can lead to a conflict warning
2049 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
2050 //realise that two separate implementations of the broken function are the same and so,
2051 //quite rightly, it shouts.
2052  /*void operator=(const TPVDElementWithContinuousPressure<DIM>&)
2053  {
2054  BrokenCopy::broken_assign("TPVDElementWithContinuousPressure");
2055  }*/
2056 
2057  /// \short Set the value at which the solid pressure is stored in the nodes
2058  inline int solid_p_nodal_index() const {return 0;}
2059 
2060  /// \short Number of values (pinned or dofs) required at node n. Can
2061  /// be overwritten for hanging node version
2062  inline virtual unsigned required_nvalue(const unsigned &n) const
2063  {return Initial_Nvalue[n];}
2064 
2065  /// Return the l-th pressure value, make sure to use the hanging
2066  /// representation if there is one!
2067  double solid_p(const unsigned &l)
2068  {return this->nodal_value(Pconv[l],this->solid_p_nodal_index());}
2069 
2070  /// Set the l-th solid pressure value to p_value
2071  void set_solid_p(const unsigned &l, const double &p_value)
2072  {this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(),p_value);}
2073 
2074  /// Return number of pressure values
2075  unsigned npres_solid() const;
2076 
2077  /// Fix the pressure dof l to be the value pvalue
2078  void fix_solid_pressure(const unsigned &l, const double &pvalue)
2079  {
2080  this->node_pt(Pconv[l])->pin(this->solid_p_nodal_index());
2081  this->node_pt(Pconv[l])->set_value(this->solid_p_nodal_index(),pvalue);
2082  }
2083 
2084  /// Generic FiniteElement output function
2085  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
2086 
2087  /// PVDEquationsWithPressure output function
2088  void output(std::ostream &outfile, const unsigned &n_plot)
2089  {PVDEquationsWithPressure<DIM>::output(outfile,n_plot);}
2090 
2091 
2092  /// C-style generic FiniteElement output function
2093  void output(FILE* file_pt) {FiniteElement::output(file_pt);}
2094 
2095  /// C-style PVDEquationsWithPressure output function
2096  void output(FILE* file_pt, const unsigned &n_plot)
2097  {PVDEquationsWithPressure<DIM>::output(file_pt,n_plot);}
2098 
2099  /// \short Order of recovery shape functions for Z2 error estimation:
2100  /// Same order as shape functions.
2101  unsigned nrecovery_order() {return 2;}
2102 
2103  /// \short Number of vertex nodes in the element
2104  unsigned nvertex_node() const
2105  {return TElement<DIM,3>::nvertex_node();}
2106 
2107  /// \short Pointer to the j-th vertex node in the element
2108  Node* vertex_node_pt(const unsigned& j) const
2109  {return TElement<DIM,3>::vertex_node_pt(j);}
2110 
2111  /// Number of 'flux' terms for Z2 error estimation
2113  {
2114  // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
2115  return DIM + DIM*(DIM-1)/2;
2116  }
2117 
2118  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
2119  /// in strain tensor.
2121  {
2122 #ifdef PARANOID
2123  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
2124  if (flux.size()!=num_entries)
2125  {
2126  std::ostringstream error_message;
2127  error_message << "The flux vector has the wrong number of entries, "
2128  << flux.size() << ", whereas it should be "
2129  << num_entries << std::endl;
2130  throw OomphLibError(error_message.str(),
2131  OOMPH_CURRENT_FUNCTION,
2132  OOMPH_EXCEPTION_LOCATION);
2133  }
2134 #endif
2135 
2136  // Get strain matrix
2137  DenseMatrix<double> strain(DIM);
2138  this->get_strain(s,strain);
2139 
2140  // Pack into flux Vector
2141  unsigned icount=0;
2142 
2143  // Start with diagonal terms
2144  for(unsigned i=0;i<DIM;i++)
2145  {
2146  flux[icount]=strain(i,i);
2147  icount++;
2148  }
2149  //Off diagonals row by row
2150  for(unsigned i=0;i<DIM;i++)
2151  {
2152  for(unsigned j=i+1;j<DIM;j++)
2153  {
2154  flux[icount]=strain(i,j);
2155  icount++;
2156  }
2157  }
2158  }
2159 
2160 
2161 };
2162 
2163 //Inline functions
2164 
2165 
2166 
2167 //==========================================================================
2168 /// 2D :
2169 /// Number of pressure values
2170 //==========================================================================
2171 template<>
2173 {return 3;}
2174 
2175 //==========================================================================
2176 /// 3D :
2177 /// Number of pressure values
2178 //==========================================================================
2179 template<>
2181 {return 4;}
2182 
2183 
2184 //==========================================================================
2185 /// 2D :
2186 /// Pressure shape functions
2187 //==========================================================================
2188 template<>
2190  const Vector<double> &s, Shape &psi) const
2191 {
2192  psi[0] = s[0];
2193  psi[1] = s[1];
2194  psi[2] = 1.0-s[0]-s[1];
2195 }
2196 
2197 //==========================================================================
2198 /// 3D :
2199 /// Pressure shape functions
2200 //==========================================================================
2201 template<>
2203  const Vector<double> &s, Shape &psi) const
2204 {
2205  psi[0] = s[0];
2206  psi[1] = s[1];
2207  psi[2] = s[2];
2208  psi[3] = 1.0-s[0]-s[1]-s[2];
2209 }
2210 
2211 
2212 //=======================================================================
2213 /// Face geometry of the 2D Taylor_Hood elements
2214 //=======================================================================
2215 template<>
2217 public virtual SolidTElement<1,3>
2218 {
2219  public:
2220 
2221  /// Constructor: Call constructor of base
2223 };
2224 
2225 //=======================================================================
2226 /// Face geometry of the 3D Taylor_Hood elements
2227 //=======================================================================
2228 template<>
2230 public virtual SolidTElement<2,3>
2231 {
2232 
2233  public:
2234 
2235  /// Constructor: Call constructor of base
2237 };
2238 
2239 
2240 
2241 //////////////////////////////////////////////////////////////////////
2242 //////////////////////////////////////////////////////////////////////
2243 //////////////////////////////////////////////////////////////////////
2244 
2245 
2246 //====================================================================
2247 /// Namespace for solid mechanics helper functions
2248 //====================================================================
2249 namespace SolidHelpers
2250 {
2251 
2252  /// \short Document the principal stresses in a 2D SolidMesh
2253  /// pointed to by \c mesh_pt, in the directory specified
2254  /// by the DocInfo object, in a format that can be processed with
2255  /// tecplot macro.
2256  template<class ELEMENT>
2257  void doc_2D_principal_stress(DocInfo& doc_info, SolidMesh* mesh_pt)
2258  {
2259  // Output principal stress vectors at the centre of all elements
2260  std::ofstream pos_file;
2261  std::ofstream neg_file;
2262  std::ostringstream filename;
2263  filename << doc_info.directory() << "/pos_principal_stress"
2264  << doc_info.number() << ".dat";
2265  pos_file.open(filename.str().c_str());
2266  filename.str("");
2267  filename << doc_info.directory() << "/neg_principal_stress"
2268  << doc_info.number() << ".dat";
2269  neg_file.open(filename.str().c_str());
2270 
2271  // Write dummy data in both so there's at lest one zone in each
2272  pos_file << 0.0 << " " << 0.0 << " " << 0.0 << " " << 0.0 << " " << std::endl;
2273  neg_file << 0.0 << " " << 0.0 << " " << 0.0 << " " << 0.0 << " " << std::endl;
2274 
2275 
2276  Vector<double> s(2);
2277  Vector<double> x(2);
2278  s[0]=0.0;
2279  s[1]=0.0;
2280  unsigned n_solid_element=mesh_pt->nelement();
2281  for (unsigned e=0;e<n_solid_element;e++)
2282  {
2283  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
2284 
2285  // Get principal stress
2286  DenseMatrix<double> principal_stress_vector(2);
2287  Vector<double> principal_stress(2);
2288  el_pt->get_principal_stress(s,principal_stress_vector,principal_stress);
2289 
2290  // Get position of centre of element
2291  el_pt->interpolated_x(s,x);
2292 
2293  // compute vectors at 45 degree for nearly hydrostatic pressure state
2294  DenseMatrix<double> rot(2);
2295 
2296  bool hydrostat=false;
2297 
2298  // Max. relative difference between principal stresses
2299  // required to classify stress state as non-hydrostatic: 1%
2300  double dev_max=1.0e-2;
2301  if (principal_stress[0]!=0.0)
2302  {
2303  if (std::fabs((principal_stress[0]-principal_stress[1])/
2304  principal_stress[0])<dev_max)
2305  {
2306  hydrostat=true;
2307  double Cos=cos(0.25*3.14159);
2308  double Sin=sin(0.25*3.14159);
2309  rot(0,0) =
2310  Cos*principal_stress_vector(0,0) - Sin*principal_stress_vector(0,1);
2311  rot(0,1) =
2312  Sin*principal_stress_vector(0,0) + Cos*principal_stress_vector(0,1);
2313  rot(1,0) =
2314  Cos*principal_stress_vector(1,0) - Sin*principal_stress_vector(1,1);
2315  rot(1,1) =
2316  Sin*principal_stress_vector(1,0) + Cos*principal_stress_vector(1,1);
2317  }
2318  }
2319 
2320  // Loop over two principal stresses:
2321  for (unsigned i=0;i<2;i++)
2322  {
2323  if (principal_stress[i]>0.0)
2324  {
2325  pos_file << x[0] << " " << x[1] << " "
2326  << principal_stress_vector(i,0) << " "
2327  << principal_stress_vector(i,1) << std::endl;
2328  pos_file << x[0] << " " << x[1] << " "
2329  << -principal_stress_vector(i,0) << " "
2330  << -principal_stress_vector(i,1) << std::endl;
2331  if (hydrostat)
2332  {
2333  pos_file << x[0] << " " << x[1] << " "
2334  << rot(i,0) << " "
2335  << rot(i,1) << std::endl;
2336  pos_file << x[0] << " " << x[1] << " "
2337  << -rot(i,0) << " "
2338  << -rot(i,1) << std::endl;
2339  }
2340  }
2341  else
2342  {
2343  neg_file << x[0] << " " << x[1] << " "
2344  << principal_stress_vector(i,0) << " "
2345  << principal_stress_vector(i,1) << std::endl;
2346  neg_file << x[0] << " " << x[1] << " "
2347  << -principal_stress_vector(i,0) << " "
2348  << -principal_stress_vector(i,1) << std::endl;
2349  if (hydrostat)
2350  {
2351  neg_file << x[0] << " " << x[1] << " "
2352  << rot(i,0) << " "
2353  << rot(i,1) << std::endl;
2354  neg_file << x[0] << " " << x[1] << " "
2355  << -rot(i,0) << " "
2356  << -rot(i,1) << std::endl;
2357  }
2358  }
2359  }
2360  }
2361 
2362  pos_file.close();
2363  neg_file.close();
2364 
2365  }
2366 
2367 };
2368 
2369 
2370 
2371 //==========================================================
2372 /// PVDElementWithContinuousPressure upgraded to become projectable
2373 //==========================================================
2374  template<class PVD_ELEMENT>
2376  public virtual ProjectableElement<PVD_ELEMENT>
2377  {
2378 
2379  public:
2380 
2381  /// \short Constructor [this was only required explicitly
2382  /// from gcc 4.5.2 onwards...]
2384 
2385 
2386  /// \short Specify the values associated with field fld.
2387  /// The information is returned in a vector of pairs which comprise
2388  /// the Data object and the value within it, that correspond to field fld.
2389  /// In the underlying PVD elements the pressures (the first
2390  /// field) are the first values at the vertex nodes etc.
2392  {
2393  // Create the vector
2394  Vector<std::pair<Data*,unsigned> > data_values;
2395 
2396  // Loop over all vertex nodes
2397  const unsigned n_solid_pres=this->npres_solid();
2398  for(unsigned j=0;j<n_solid_pres;j++)
2399  {
2400  // Add the data value associated with the pressure components
2401  unsigned vertex_index=this->Pconv[j];
2402  data_values.push_back(std::make_pair(this->node_pt(vertex_index),0));
2403  }
2404 
2405  // Return the vector
2406  return data_values;
2407  }
2408 
2409  /// \short Number of fields to be projected: 1, corresponding to
2410  /// the pressure only
2411  unsigned nfields_for_projection() {return 1;}
2412 
2413  /// \short Number of history values to be stored for fld-th field.
2414  /// (Includes the current value!)
2415  unsigned nhistory_values_for_projection(const unsigned &fld)
2416  {
2417  //pressure doesn't have history values as such so only one value
2418  // representing the current value
2419  return 1;
2420  }
2421 
2422  ///\short Number of positional history values (Includes the current value!)
2424  {
2425  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
2426  }
2427 
2428  /// \short Return Jacobian of mapping and shape functions of field fld
2429  /// at local coordinate s
2430  double jacobian_and_shape_of_field(const unsigned &fld,
2431  const Vector<double> &s,
2432  Shape &psi)
2433  {
2434  //Get the solid pressure shape function
2435  this->solid_pshape(s,psi);
2436  //Return the Jacobian of the eulerian mapping
2437  return this->J_eulerian(s);
2438  }
2439 
2440  /// \short Return interpolated field fld at local coordinate s, at time level
2441  /// t (t=0: present; t>0: history values)
2442  double get_field(const unsigned &t,
2443  const unsigned &fld,
2444  const Vector<double>& s)
2445  {
2446  return this->interpolated_solid_p(s);
2447  }
2448 
2449 
2450 
2451  ///Return number of values in field fld
2452  unsigned nvalue_of_field(const unsigned &fld)
2453  {
2454  return this->npres_solid();
2455  }
2456 
2457 
2458  ///Return local equation number of value j in field fld.
2459  int local_equation(const unsigned &fld,
2460  const unsigned &j)
2461  {
2462  return this->solid_p_local_eqn(j);
2463  }
2464 
2465  };
2466 
2467 
2468 //=======================================================================
2469 /// Face geometry for element is the same as that for the underlying
2470 /// wrapped element
2471 //=======================================================================
2472  template<class ELEMENT>
2474  : public virtual FaceGeometry<ELEMENT>
2475  {
2476  public:
2477  FaceGeometry() : FaceGeometry<ELEMENT>() {}
2478  };
2479 
2480 
2481 //=======================================================================
2482 /// Face geometry of the Face Geometry for element is the same as
2483 /// that for the underlying wrapped element
2484 //=======================================================================
2485  template<class ELEMENT>
2488  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
2489  {
2490  public:
2492  };
2493 
2494 
2495 
2496 }
2497 
2498 #endif
2499 
2500 
2501 
2502 
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:386
unsigned ndof_types() const
returns the number of DOF types associated with this element.
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
HermitePVDElement()
Constructor, there are no internal data points.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void output(std::ostream &outfile)
Generic FiniteElement output function.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
void output(std::ostream &outfile)
Output function.
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...
unsigned nvertex_node() const
Number of vertex nodes in the element.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Fill in contribution to Mass matrix and Jacobian (either by FD or analytically, for the positional va...
void get_energy(double &pot_en, double &kin_en)
Get potential (strain) and kinetic energy.
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n. Can be overwritten for hanging node version...
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
virtual void calculate_d_second_piola_kirchhoff_stress_dG(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, RankFourTensor< double > &d_sigma_dG, const bool &symmetrize_tensor=true)
Calculate the derivatives of the contravariant 2nd Piola Kirchhoff stress tensor with respect to the ...
virtual void fill_in_generic_contribution_to_residuals_pvd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
void disable_evaluate_jacobian_by_fd()
Set Jacobian to be evaluated analytically Else: by FD.
void unpin_elemental_solid_pressure_dofs()
Unpin all solid pressure dofs in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
virtual void pin_elemental_redundant_nodal_solid_pressures()
Pin the element's redundant solid pressures (needed for refinement)
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void output(std::ostream &outfile)
Generic FiniteElement output function.
bool is_jacobian_evaluated_by_fd() const
Return the flag indicating whether the jacobian is evaluated by fd.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
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
static double Default_lambda_sq_value
Static default value for timescale ratio (1.0 – for natural scaling)
void output(FILE *file_pt, const unsigned &n_plot)
C-style PVDEquationsWithPressure output function.
void set_solid_p(const unsigned &l, const double &p_value)
Set the l-th solid pressure value to p_value.
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
Definition: elements.cc:6197
PVDElementWithContinuousPressure upgraded to become projectable.
static int Solid_pressure_not_stored_at_node
Static "magic" number that indicates that the solid pressure is not stored at a node.
Information for documentation of results: Directory and file number to enable output in the form RESL...
PrestressFctPt & prestress_fct_pt()
Access function: Pointer to pre-stress function.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
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...
cstr elem_len * i
Definition: cfortran.h:607
int solid_p_nodal_index() const
Set the value at which the solid pressure is stored in the nodes.
void output(FILE *file_pt)
C-style output function.
unsigned npres_solid() const
Return number of pressure values.
void output(FILE *file_pt)
C-style output function.
void output(FILE *file_pt)
C-style generic FiniteElement output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
virtual double solid_p(const unsigned &l)=0
Return the lth solid pressure.
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
ConstitutiveLaw * Constitutive_law_pt
Pointer to the constitutive law.
virtual unsigned npres_solid() const
Return the number of solid pressure degrees of freedom Default is that there are no solid pressures...
double solid_p(const unsigned &l)
Return the lth pressure value.
virtual void get_isotropic_growth(const unsigned &ipt, const Vector< double > &s, const Vector< double > &xi, double &gamma) const
Evaluate isotropic growth function at Lagrangian coordinate xi and/or local coordinate s...
void get_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &Gcontra, double &detG)
Return the deviatoric part of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitut...
void enable_inertia()
Switch on solid inertia.
void disable_inertia()
Switch off solid inertia.
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
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain tensor.
unsigned nfields_for_projection()
Number of fields to be projected: 1, corresponding to the pressure only.
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
void get_d_stress_dG_upper(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, RankFourTensor< double > &d_sigma_dG)
Return the derivatives of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive ...
void output(std::ostream &outfile, const unsigned &n_plot)
PVDEquationsWithPressure output function.
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 ...
char t
Definition: cfortran.h:572
IsotropicGrowthFctPt isotropic_growth_fct_pt() const
Access function: Pointer to isotropic growth function (const version)
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Includes the current value!)
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
void set_solid_p(const unsigned &l, const double &p_value)
Set the l-th solid pressure value to p_value.
FaceGeometry()
Constructor must call constructor of the underlying Solid element.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
QPVDElement()
Constructor, there are no internal data points.
void output(std::ostream &outfile)
Generic FiniteElement output function.
virtual int solid_p_local_eqn(const unsigned &i) const
Return the local degree of freedom associated with the i-th solid pressure. Default is that there are...
int solid_p_local_eqn(const unsigned &i) const
Overload the access function that is used to return local equation corresponding to the i-th solid pr...
void solid_pshape(const Vector< double > &s, Shape &psi) const
Return the pressure shape functions.
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
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...
PVDElementWithContinuousPressure upgraded to become projectable.
FaceGeometry()
Constructor must call constructor of the underlying element.
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
e
Definition: cfortran.h:575
PrestressFctPt Prestress_fct_pt
Pointer to prestress function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to Jacobian (either by FD or analytically, control this via evaluate_jacobian_by...
A Rank 4 Tensor class.
Definition: matrices.h:1625
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
void output(FILE *file_pt)
C-style output: x,y,[z],xi0,xi1,[xi2],p,gamma.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Includes the current value!) ...
ProjectablePVDElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
ProjectablePVDElementWithContinuousPressure()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
unsigned npres_solid() const
Return number of pressure values.
void(* BodyForceFctPt)(const double &t, const Vector< double > &xi, Vector< double > &b)
Function pointer to function that specifies the body force as a function of the Lagrangian coordinate...
unsigned & number()
Number used (e.g.) for labeling output files.
FaceGeometry()
Constructor must call constructor of underlying solid element.
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
Definition: multi_domain.cc:66
virtual void fill_in_generic_residual_contribution_pvd_with_pressure(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
Returns the residuals for the discretised principle of virtual displacements, formulated in the incom...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals.
void fill_in_jacobian_for_solid_ic(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the residuals and Jacobian for the setup of an initial condition. The global equations are: ...
Definition: elements.h:3777
void output(std::ostream &outfile, const unsigned &n_plot)
SolidQHermiteElement output function.
PVDEquationsWithPressure()
Constructor, by default the element is NOT incompressible.
bool Unsteady
Flag that switches inertia on/off.
PVDEquations()
Constructor.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
virtual void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma)=0
Return the 2nd Piola Kirchoff stress tensor, as calculated from the constitutive law at specified loc...
void output(FILE *file_pt, const unsigned &n_plot)
C-style PVDEquationsWithPressure output function.
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 doc_2D_principal_stress(DocInfo &doc_info, SolidMesh *mesh_pt)
Document the principal stresses in a 2D SolidMesh pointed to by mesh_pt, in the directory specified b...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain tensor.
virtual void set_solid_p(const unsigned &l, const double &p_value)=0
Set the lth solid pressure to p_value.
void fix_solid_pressure(const unsigned &l, const double &pvalue)
Fix the pressure dof l to be the value pvalue.
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.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Includes the current value!)
unsigned npres_solid() const
Return number of pressure values.
void get_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Return the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive law: Pass metric te...
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
void set_solid_p(const unsigned &l, const double &p_value)
Set the l-th pressure value to p_value.
TPVDElement()
Constructor, there are no internal data points.
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
virtual void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)=0
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
FaceGeometry()
Constructor must call constructor of the underlying Solid element.
void output(std::ostream &outfile)
Overload the output function.
virtual void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Use finite differences to calculate the Jacobian entries corresponding to the solid positions...
Definition: elements.cc:6658
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
int solid_p_local_eqn(const unsigned &i) const
Overload the access function that is used to return local equation corresponding to the i-th solid pr...
virtual void solid_pshape(const Vector< double > &s, Shape &psi) const =0
Return the solid pressure shape functions.
static char t char * s
Definition: cfortran.h:572
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
FaceGeometry()
Constructor: Call constructor of base.
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma)
Return the 2nd Piola Kirchoff stress tensor, as calculated from the constitutive law at specified loc...
void extended_output(std::ostream &outfile, const unsigned &n_plot)
Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress components.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
void get_d_stress_dG_upper(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, const double &detG, const double &interpolated_solid_p, RankFourTensor< double > &d_sigma_dG, DenseMatrix< double > &d_detG_dG)
Return the derivative of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive l...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
IsotropicGrowthFctPt & isotropic_growth_fct_pt()
Access function: Pointer to isotropic growth function.
bool is_inertia_enabled() const
Access function to flag that switches inertia on/off (const version)
void fix_solid_pressure(const unsigned &l, const double &pvalue)
Fix the pressure dof l to be the value pvalue.
void unpin_elemental_solid_pressure_dofs()
Unpin all solid pressure dofs in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
void get_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma_dev, DenseMatrix< double > &Gcontra, double &gen_dil, double &inv_kappa)
Return the deviatoric part of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitut...
FaceGeometry()
Constructor must call constructor of underlying solid element.
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 fix_solid_pressure(const unsigned &l, const double &pvalue)
Fix the pressure dof l to be the value pvalue.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals for the solid equations (the discretised principle of virtual displacements) ...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
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...
FaceGeometry()
Constructor must call constructor of underlying solid element.
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
void set_compressible()
Set the material to be compressible.
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
void output(std::ostream &outfile, const unsigned &n_plot)
PVDEquationsWithPressure output function.
double(* PrestressFctPt)(const unsigned &i, const unsigned &j, const Vector< double > &xi)
Function pointer to function that specifies the pre-stress sigma_0(i,j) as a function of the Lagrangi...
bool Incompressible
Boolean to determine whether the solid is incompressible or not.
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
void solid_pshape_at_knot(const unsigned &ipt, Shape &psi) const
Return the stored solid shape functions at the knots.
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
unsigned nfields_for_projection()
Number of fields to be projected: 0.
double prestress(const unsigned &i, const unsigned &j, const Vector< double > xi)
Return (i,j)-th component of second Piola Kirchhoff membrane prestress at Lagrangian coordinate xi...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field (Includes the current value!). No nodal data.
QPVDElementWithPressure()
Constructor, there are DIM+1 internal data points.
void get_d_stress_dG_upper(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, const double &gen_dil, const double &inv_kappa, const double &interpolated_solid_p, RankFourTensor< double > &d_sigma_dG, DenseMatrix< double > &d_gen_dil_dG)
Return the derivative of the deviatoric part of the 2nd Piola Kirchhoff stress tensor, as calculated from the constitutive law in the nearly incompresible formulation. Also return the derivative of the generalised dilatation.
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.
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
void enable_evaluate_jacobian_by_fd()
Set Jacobian to be evaluated by FD? Else: Analytically.
void get_deformed_covariant_basis_vectors(const Vector< double > &s, DenseMatrix< double > &def_covariant_basis)
Return the deformed covariant basis vectors at specified local coordinate: def_covariant_basis(i,j) is the j-th component of the i-th basis vector.
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
bool Evaluate_jacobian_by_fd
Use FD to evaluate Jacobian.
virtual int solid_p_nodal_index() const
Return the index at which the solid pressure is stored if it is stored at the nodes. If not stored at the nodes this will return a negative number.
void extended_output(std::ostream &outfile, const unsigned &n_plot)
Output: x,y,[z],xi0,xi1,[xi2],gamma and the strain and stress.
void unpin_elemental_solid_pressure_dofs()
Unpin all solid pressure dofs in the element.
void fill_in_jacobian_for_newmark_accel(DenseMatrix< double > &jacobian)
Fill in the contributions of the Jacobian matrix for the consistent assignment of the initial "accele...
Definition: elements.cc:6900
bool is_incompressible() const
Return whether the material is incompressible.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
General SolidMesh class.
Definition: mesh.h:2252
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
TPVDElementWithContinuousPressure(const TPVDElementWithContinuousPressure< DIM > &dummy)
Broken copy constructor.
IsotropicGrowthFctPt Isotropic_growth_fct_pt
Pointer to isotropic growth function.
PVDEquationsBase()
Constructor: Set null pointers for constitutive law and for isotropic growth function. Set physical parameter values to default values, enable inertia and set body force to zero. Default evaluation of Jacobian: analytically rather than by FD.
bool Solve_for_consistent_newmark_accel_flag
Flag to indicate which system of equations to solve when assigning initial conditions for time-depend...
Definition: elements.h:4049
std::string directory() const
Output directory.
FaceGeometry()
Constructor must call constructor of underlying solid element.
void output(std::ostream &outfile)
SolidQHermiteElement output function.
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
void get_stress(const Vector< double > &s, DenseMatrix< double > &sigma)
Return the 2nd Piola Kirchoff stress tensor, as calculated from the constitutive law at specified loc...
unsigned ndof_types() const
returns the number of DOF types associated with this element: displacement components and pressure ...
BodyForceFctPt body_force_fct_pt() const
Access function: Pointer to body force function (const version)
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
void get_mass_matrix_diagonal(Vector< double > &mass_diag)
Compute the diagonal of the displacement mass matrix for LSC preconditioner.
static void unpin_all_solid_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void output(FILE *file_pt)
C-style output function.
void output(std::ostream &outfile)
Output: x,y,[z],xi0,xi1,[xi2],gamma.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
bool has_internal_solid_data()
There is internal solid data so we can't use the automatic assignment of consistent initial condition...
void output(std::ostream &outfile)
Output function.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void output(FILE *file_pt)
C-style output: x,y,[z],xi0,xi1,[xi2],gamma.
SolidInitialCondition * Solid_ic_pt
Pointer to object that specifies the initial condition.
Definition: elements.h:3874
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void body_force(const Vector< double > &xi, Vector< double > &b) const
Evaluate body force at Lagrangian coordinate xi at present time (returns zero vector if no body force...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:66
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n. Can be overwritten for hanging node version...
TPVDBubbleEnrichedElement()
Constructor, there are no internal data points.
ConstitutiveLaw *& constitutive_law_pt()
Return the constitutive law pointer.
BodyForceFctPt Body_force_fct_pt
Pointer to body force function.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain tensor.
void solid_pshape(const Vector< double > &s, Shape &psi) const
Return the pressure shape functions.
void unpin_elemental_solid_pressure_dofs()
Unpin all solid pressure dofs – empty as there are no pressures.
FaceGeometry()
Constructor must call constructor of the underlying Point element.
FaceGeometry()
Constructor: Call constructor of base.
unsigned P_solid_internal_index
Internal index that indicates at which internal data value the solid presure is stored.
void set_incompressible()
Set the material to be incompressible.
void solid_pshape(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void output(std::ostream &outfile)
Output: x,y,[z],xi0,xi1,[xi2],p,gamma.
void output(FILE *file_pt, const unsigned &n_plot)
C-style SolidQHermiteElement output function.
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
double interpolated_solid_p(const Vector< double > &s)
Return the interpolated_solid_pressure.
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
Definition: elements.h:3881
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
void output(std::ostream &outfile, const unsigned &n_plot)
PVDEquationsWithPressure output function.
void output(FILE *file_pt)
C-style SolidQHermiteElement output function.
SolidFiniteElement class.
Definition: elements.h:3320
virtual void unpin_elemental_solid_pressure_dofs()=0
Unpin all solid pressure dofs in the element.
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
void get_principal_stress(const Vector< double > &s, DenseMatrix< double > &principal_stress_vector, Vector< double > &principal_stress)
Compute principal stress vectors and (scalar) principal stresses at specified local coordinate...
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
int solid_p_local_eqn(const unsigned &i) const
Overload the access function that is used to return local equation corresponding to the i-th solid pr...
void(* IsotropicGrowthFctPt)(const Vector< double > &xi, double &gamma)
Function pointer to function that specifies the isotropic growth as a function of the Lagrangian coor...
void output(FILE *file_pt)
C-style generic FiniteElement output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style PVDEquationsWithPressure output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to Jacobian (either by FD or analytically, for the positional variables; control...
FaceGeometry()
Constructor must call the constructor of the underlying solid element.
int solid_p_nodal_index() const
Broken assignment operator.
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:596
void output(FILE *file_pt)
C-style Generic FiniteElement output function.
void output(std::ostream &outfile)
Output function.