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-2022 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
48namespace 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,
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)
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
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 {
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)
419
420 /// Flag that switches inertia on/off
422
423 /// Pointer to body force function
425
426 /// Static default value for timescale ratio (1.0 -- for natural scaling)
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 {
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
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
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 {
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 {
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>
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>
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
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
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 {
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 {
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 {
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 {
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 {
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 {
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.
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 {
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 {
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.
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 {
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 {
2216 }
2217
2218 /// Order of recovery shape functions for Z2 error estimation:
2219 /// Same order as shape functions.
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>
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
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
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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.
ConstitutiveLaw *& constitutive_law_pt()
Return the constitutive law pointer.
static void unpin_all_solid_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
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)
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.
const double & lambda_sq() const
Access function for timescale ratio (nondim density)
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(* IsotropicGrowthFctPt)(const Vector< double > &xi, double &gamma)
Function pointer to function that specifies the isotropic growth as a function of the Lagrangian coor...
PrestressFctPt & prestress_fct_pt()
Access function: Pointer to pre-stress function.
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....
IsotropicGrowthFctPt & isotropic_growth_fct_pt()
Access function: Pointer to isotropic growth function.
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)
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.
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
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.
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 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 enable_inertia()
Switch on solid inertia.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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_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 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...
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.
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 ...
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.
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.
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 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!)
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.
unsigned nvertex_node() const
Number of vertex nodes in the element.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node 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...
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.
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 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.
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.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node 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...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...