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