refineable_solid_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for refineable solid mechanics elements
31 
32 //Include guard to prevent multiple inclusions of this header
33 #ifndef OOMPH_REFINEABLE_ELASTICITY_ELEMENTS_HEADER
34 #define OOMPH_REFINEABLE_ELASTICITY_ELEMENTS_HEADER
35 
36 //oomph-lib headers
37 #include "solid_elements.h"
38 #include "../generic/refineable_quad_element.h"
39 #include "../generic/refineable_brick_element.h"
40 #include "../generic/error_estimator.h"
41 
42 namespace oomph
43 {
44 
45 //========================================================================
46 /// Class for Refineable PVD equations
47 //========================================================================
48 template<unsigned DIM>
49 class RefineablePVDEquations: public virtual PVDEquations<DIM>,
50 public virtual RefineableSolidElement,
51 public virtual ElementWithZ2ErrorEstimator
52 {
53 public:
54 
55  /// \short Constructor
57  PVDEquations<DIM>(),
61  {}
62 
63  /// Call the residuals including hanging node cases
65  Vector<double> &residuals,
66  DenseMatrix<double> &jacobian,
67  const unsigned& flag);
68 
69  /// No values are interpolated in this element (pure solid)
70  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
71  Vector<double>& values) {values.clear();}
72 
73  /// No values are interpolated in this element (pure solid)
75  {values.clear();}
76 
77  /// Number of 'flux' terms for Z2 error estimation
78  unsigned num_Z2_flux_terms()
79  {
80  // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
81  return DIM + DIM*(DIM-1)/2;
82  }
83 
84  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
85  /// in strain tensor.
87 {
88 #ifdef PARANOID
89  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
90  if (flux.size()!=num_entries)
91  {
92  std::ostringstream error_message;
93  error_message << "The flux vector has the wrong number of entries, "
94  << flux.size() << ", whereas it should be "
95  << num_entries << std::endl;
96  throw OomphLibError(error_message.str(),
97  OOMPH_CURRENT_FUNCTION,
98  OOMPH_EXCEPTION_LOCATION);
99  }
100 #endif
101 
102  // Get strain matrix
103  DenseMatrix<double> strain(DIM);
104  this->get_strain(s,strain);
105 
106  // Pack into flux Vector
107  unsigned icount=0;
108 
109  // Start with diagonal terms
110  for(unsigned i=0;i<DIM;i++)
111  {
112  flux[icount]=strain(i,i);
113  icount++;
114  }
115 
116  //Off diagonals row by row
117  for(unsigned i=0;i<DIM;i++)
118  {
119  for(unsigned j=i+1;j<DIM;j++)
120  {
121  flux[icount]=strain(i,j);
122  icount++;
123  }
124  }
125 }
126 
127  /// Number of continuously interpolated values: 0 (pure solid problem)
128  unsigned ncont_interpolated_values() const {return 0;}
129 
130  //Return a pointer to the solid node at which pressure dof l2 is stored
131  //This is only required so that the generic templating in
132  //PseudoSolidNodeUpdateElements works OK
133  virtual Node* solid_pressure_node_pt(const unsigned &l) {return 0;}
134 
135  /// Further build function, pass the pointers down to the sons
137  {
138  RefineablePVDEquations<DIM>* cast_father_element_pt
139  = dynamic_cast<RefineablePVDEquations<DIM>*>
140  (this->father_element_pt());
141 
142  // Do whatever needs to be done in the base class
144 
145  //Set pointer to isotropic growth function
146  this->Isotropic_growth_fct_pt =
147  cast_father_element_pt->isotropic_growth_fct_pt();
148 
149  // Set pointer to body force function
150  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
151 
152  // Set pointer to the contitutive law
153  this->Constitutive_law_pt = cast_father_element_pt->constitutive_law_pt();
154 
155  // Set the timescale ratio (non-dim. density)
156  this->Lambda_sq_pt = cast_father_element_pt->lambda_sq_pt();
157 
158  // Set the flag that switches inertia on/off
159  this->Unsteady = cast_father_element_pt->is_inertia_enabled();
160 
161  // Evaluation of Jacobian by same method as father
163  cast_father_element_pt->is_jacobian_evaluated_by_fd();
164  }
165 
166 };
167 
168 //========================================================================
169 /// Class for refineable QPVDElement elements
170 //========================================================================
171 template<unsigned DIM, unsigned NNODE_1D>
172 class RefineableQPVDElement : public virtual QPVDElement<DIM,NNODE_1D>,
173  public virtual RefineablePVDEquations<DIM>,
174  public virtual RefineableSolidQElement<DIM>
175 {
176 public:
177 
178  /// \short Constructor:
180  QPVDElement<DIM,NNODE_1D>(),
183  RefineablePVDEquations<DIM>(),
184  RefineableSolidQElement<DIM>() {}
185 
186  /// Empty rebuild from sons, no need to reconstruct anything here
187  void rebuild_from_sons(Mesh* &mesh_pt) {}
188 
189  /// \short Number of vertex nodes in the element
190  unsigned nvertex_node() const
192 
193  /// \short Pointer to the j-th vertex node in the element
194  Node* vertex_node_pt(const unsigned& j) const
196 
197  /// \short Order of recovery shape functions for Z2 error estimation:
198  /// Same order as shape functions.
199  unsigned nrecovery_order() {return NNODE_1D-1;}
200 
201  /// \short No additional hanging node procedures are required
202  /// for the solid elements.
204 
205 };
206 
207 //==============================================================
208 /// FaceGeometry of the 2D RefineableQPVDElement elements
209 //==============================================================
210 template<unsigned NNODE_1D>
212  public virtual SolidQElement<1,NNODE_1D>
213 {
214  public:
215  //Make sure that we call the constructor of the SolidQElement
216  //Only the Intel compiler seems to need this!
217  FaceGeometry() : SolidQElement<1,NNODE_1D>() {}
218 };
219 
220 //==============================================================
221 /// FaceGeometry of the FaceGeometry of the 2D RefineableQPVDElement
222 //==============================================================
223 template<unsigned NNODE_1D>
225  public virtual PointElement
226 {
227  public:
228  //Make sure that we call the constructor of the SolidQElement
229  //Only the Intel compiler seems to need this!
231 };
232 
233 
234 //==============================================================
235 /// FaceGeometry of the 3D RefineableQPVDElement elements
236 //==============================================================
237 template<unsigned NNODE_1D>
239  public virtual SolidQElement<2,NNODE_1D>
240 {
241  public:
242  //Make sure that we call the constructor of the SolidQElement
243  //Only the Intel compiler seems to need this!
244  FaceGeometry() : SolidQElement<2,NNODE_1D>() {}
245 };
246 
247 //==============================================================
248 /// FaceGeometry of the FaceGeometry of the 3D RefineableQPVDElement
249 //==============================================================
250 template<unsigned NNODE_1D>
252 public virtual SolidQElement<1,NNODE_1D>
253 {
254  public:
255  //Make sure that we call the constructor of the SolidQElement
256  //Only the Intel compiler seems to need this!
257  FaceGeometry() : SolidQElement<1,NNODE_1D>() {}
258 };
259 
260 
261 //===========================================================================
262 /// Class for Refineable solid mechanics elements in near-incompressible/
263 /// incompressible formulation, so a pressure is included! In this case,
264 /// the pressure interpolation is discontinuous, a la Crouzeix Raviart
265 //===========================================================================
266 template<unsigned DIM>
268 public virtual PVDEquationsWithPressure<DIM>,
269  public virtual RefineableSolidElement,
270  public virtual ElementWithZ2ErrorEstimator
271 {
272 
273  public:
274 
275  /// \short Constructor:
281  {}
282 
283 /// \short Add element's contribution to elemental residual vector and/or
284 /// Jacobian matrix
285 /// flag=1: compute both
286 /// flag=0: compute only residual vector
288  Vector<double> &residuals, DenseMatrix<double> &jacobian,
289  DenseMatrix<double> &mass_matrix,
290  const unsigned& flag);
291 
292  /// No values are interpolated in this element (pure solid)
293  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
294  Vector<double>& values) {values.clear();}
295 
296  /// No values are interpolated in this element (pure solid)
298  {values.clear();}
299 
300  /// Number of 'flux' terms for Z2 error estimation
301  unsigned num_Z2_flux_terms()
302  {
303  // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
304  return DIM + DIM*(DIM-1)/2;
305  }
306 
307 // Get 'flux' for Z2 error recovery: Upper triangular entries
308 /// in strain tensor.
309 //----------------------------------------------------------------
311  {
312  //Find the dimension of the problem
313 #ifdef PARANOID
314  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
315  if (flux.size()!=num_entries)
316  {
317  std::ostringstream error_message;
318  error_message << "The flux vector has the wrong number of entries, "
319  << flux.size() << ", whereas it should be "
320  << num_entries << std::endl;
321  throw OomphLibError(error_message.str(),
322  OOMPH_CURRENT_FUNCTION,
323  OOMPH_EXCEPTION_LOCATION);
324  }
325 #endif
326 
327  // Get strain matrix
328  DenseMatrix<double> strain(DIM);
329  this->get_strain(s,strain);
330 
331  // Pack into flux Vector
332  unsigned icount=0;
333 
334  // Start with diagonal terms
335  for(unsigned i=0;i<DIM;i++)
336  {
337  flux[icount]=strain(i,i);
338  icount++;
339  }
340 
341  //Off diagonals row by row
342  for(unsigned i=0;i<DIM;i++)
343  {
344  for(unsigned j=i+1;j<DIM;j++)
345  {
346  flux[icount]=strain(i,j);
347  icount++;
348  }
349  }
350  }
351 
352 /// Number of continuously interpolated values: 0 (pure solid problem)
353  unsigned ncont_interpolated_values() const {return 0;}
354 
355  //Return a pointer to the solid node at which pressure dof l2 is stored
356  virtual Node* solid_pressure_node_pt(const unsigned &l) {return 0;}
357 
358 
359  ///Pass the generic stuff down to the sons
361  {
362  RefineablePVDEquationsWithPressure<DIM>* cast_father_element_pt
364  (this->father_element_pt());
365 
366  // Do whatever needs to be done in the base class
368 
369  //Set pointer to isotropic growth function
370  this->Isotropic_growth_fct_pt =
371  cast_father_element_pt->isotropic_growth_fct_pt();
372 
373  // Set pointer to body force function
374  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
375 
376  // Set pointer to the contitutive law
377  this->Constitutive_law_pt = cast_father_element_pt->constitutive_law_pt();
378 
379  // Set the timescale ratio (non-dim. density)
380  this->Lambda_sq_pt = cast_father_element_pt->lambda_sq_pt();
381 
382  // Set the flag that switches inertia on/off
383  this->Unsteady = cast_father_element_pt->is_inertia_enabled();
384 
385  // Set the incompressibility flag
386  this->Incompressible = cast_father_element_pt->is_incompressible();
387 
388  // Evaluation of Jacobian by same method as father
390  cast_father_element_pt->is_jacobian_evaluated_by_fd();
391  }
392 
393 
394  /// \short Compute the diagonal of the displacement mass matrix for
395  /// LSC preconditioner
396  void get_mass_matrix_diagonal(Vector<double> &mass_diag);
397 
398 };
399 
400 //===========================================================================
401 /// Class for refineable solid mechanics elements in near-incompressible/
402 /// incompressible formulation, so a pressure is included! In this case,
403 /// the pressure interpolation is discontinuous, a la Crouzeix Raviart,
404 /// and the displacement is always quadratic.
405 //===========================================================================
406 template<unsigned DIM>
408 public virtual QPVDElementWithPressure<DIM>,
409 public virtual RefineablePVDEquationsWithPressure<DIM>,
410 public virtual RefineableSolidQElement<DIM>
411 {
412 
413  private:
414 
415  /// Unpin all solid pressure dofs
417  {
418  unsigned n_pres = this->npres_solid();
419  // loop over pressure dofs and unpin them
420  for(unsigned l=0;l<n_pres;l++)
421  {this->internal_data_pt(this->P_solid_internal_index)->unpin(l);}
422  }
423 
424  public:
425 
426  /// \short Constructor:
432  RefineableSolidQElement<DIM>() {}
433 
434 
435  /// \short Reconstruct the pressure from the sons
436  /// Must be specialized for each dimension
437  inline void rebuild_from_sons(Mesh* &mesh_pt);
438 
439  /// \short Number of vertex nodes in the element
440  unsigned nvertex_node() const
442 
443  /// \short Pointer to the j-th vertex node in the element
444  Node* vertex_node_pt(const unsigned& j) const
446 
447  /// \short Order of recovery shape functions for Z2 error estimation:
448  /// Same order as shape functions.
449  unsigned nrecovery_order() {return 2;} //NNODE_1D-1;}
450 
451  /// \short No additional hanging node procedures are required for
452  /// discontinuous solid pressures.
454 
455  /// \short Further build: Interpolate the solid pressure values
456  /// Again this must be specialised for each dimension
457  inline void further_build();
458 
459 
460  /// Number of continuously interpolated values: 0 (pure solid problem)
461  unsigned ncont_interpolated_values() const {return 0;}
462 
463 };
464 
465 //======================================================================
466 ///FaceGeometry of the 2D RefineableQPVDElementWithPressure
467 //=======================================================================
468 template<>
470  public virtual SolidQElement<1,3>
471 {
472  public:
473  //Make sure that we call the constructor of the SolidQElement
474  //Only the Intel compiler seems to need this!
476 };
477 
478 
479 //===========================================================================
480 ///FaceGeometry of the FaceGeometry of the 2D RefineableQPVDElementWithPressure
481 //============================================================================
482 template<>
484 public virtual PointElement
485 {
486  public:
487  //Make sure that we call the constructor of the SolidQElement
488  //Only the Intel compiler seems to need this!
490 };
491 
492 //====================================================================
493 /// 2D rebuild from sons: reconstruct the solid pressure
494 //===================================================================
495 template<>
498 {
499  using namespace QuadTreeNames;
500 
501  //Storage for the solid pressure of each son
502  double centre_solid_p[4]={0.0,0.0,0.0,0.0};
503 
504  // Loop over the sons and assign the central solid pressure
505  for (unsigned ison=0;ison<4;ison++)
506  {
507  // Add the sons midnode pressure
508  centre_solid_p[ison] =
510  (quadtree_pt()->son_pt(ison)->object_pt())->solid_p(0);
511  }
512 
513  // Use the average for the central solid pressure
514  double p_value = 0.25*(centre_solid_p[0] + centre_solid_p[1] +
515  centre_solid_p[2] + centre_solid_p[3]);
516  //Actually set the pressure
517  set_solid_p(0,p_value);
518 
519 
520  //Slope in s_0 direction
521  //----------------------
522 
523  // Use average of the 2 FD approximations based on the
524  // elements central pressure values
525  // [Other options: Take average of the four
526  // pressure derivatives]
527  double slope1=centre_solid_p[SE] - centre_solid_p[SW];
528  double slope2=centre_solid_p[NE] - centre_solid_p[NW];
529 
530 
531  // Use the average value
532  p_value = 0.5*(slope1 + slope2);
533  set_solid_p(1,p_value);
534 
535  //Slope in s_1 direction
536  //----------------------
537 
538  // Use average of the 2 FD approximations based on the
539  // elements central pressure values
540  // [Other options: Take average of the four
541  // pressure derivatives]
542 
543  slope1= centre_solid_p[NE] - centre_solid_p[SE];
544  slope2= centre_solid_p[NW] - centre_solid_p[SW];
545 
546  // Use the average
547  p_value = 0.5*(slope1 + slope2);
548  set_solid_p(2,p_value);
549 }
550 
551 
552 //==============================================================
553 /// 2D further build interpolates the internal solid pressure
554 /// from the father.
555 //=============================================================
556 template<>
558 {
560 
561  using namespace QuadTreeNames;
562 
563  // What type of son am I? Ask my quadtree representation...
564  int son_type=quadtree_pt()->son_type();
565 
566  // Pointer to my father (in element impersonation)
567  RefineableElement* father_el_pt=
568  quadtree_pt()->father_pt()->object_pt();
569 
570  Vector<double> s_father(2);
571 
572  // Son midpoint is located at the following coordinates in father element:
573 
574  // South west son
575  if (son_type==SW)
576  {
577  s_father[0]=-0.5;
578  s_father[1]=-0.5;
579  }
580  // South east son
581  else if (son_type==SE)
582  {
583  s_father[0]= 0.5;
584  s_father[1]=-0.5;
585  }
586  // North east son
587  else if (son_type==NE)
588  {
589  s_father[0]=0.5;
590  s_father[1]=0.5;
591  }
592 
593  // North west son
594  else if (son_type==NW)
595  {
596  s_father[0]=-0.5;
597  s_father[1]= 0.5;
598  }
599 
600  // Pressure value in father element
601  RefineableQPVDElementWithPressure<2>* cast_father_element_pt=
602  dynamic_cast<RefineableQPVDElementWithPressure<2>*>(father_el_pt);
603  double press=cast_father_element_pt->interpolated_solid_p(s_father);
604 
605  // Pressure value gets copied straight into internal dof:
606  set_solid_p(0,press);
607 
608  // The slopes get copied from father and halved
609  set_solid_p(1,0.5*cast_father_element_pt->solid_p(1));
610  set_solid_p(2,0.5*cast_father_element_pt->solid_p(2));
611  }
612 
613 
614 //===============================================================
615 /// 3D rebuild from sons: reconstruct the pressure
616 //===============================================================
617 template<>
620 {
621  using namespace OcTreeNames;
622 
623  //Storage for the central solid pressure of each son
624  double centre_solid_p[8]= {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0};
625 
626  //Loop over the sons and assign the central solid pressure
627  for(unsigned ison=0;ison<8;ison++)
628  {
629  centre_solid_p[ison] =
630  octree_pt()->son_pt(ison)->object_pt()->
631  internal_data_pt(this->P_solid_internal_index)->value(0);
632  }
633 
634  //Central pressure value:
635  //-----------------------
636 
637  // Use average of the sons central pressure values
638  // Other options: Take average of the four (discontinuous)
639  // pressure values at the father's midpoint]
640  double av_press=0.0;
641 
642  // Loop over the sons and sum the centre pressures
643  for (unsigned ison=0;ison<8;ison++)
644  {
645  av_press += centre_solid_p[ison];
646  }
647 
648  // Use the average
649  internal_data_pt(this->P_solid_internal_index)->
650  set_value(0,0.125*av_press);
651 
652 
653  //Slope in s_0 direction
654  //----------------------
655 
656  // Use average of the 4 FD approximations based on the
657  // elements central pressure values
658  // [Other options: Take average of the four
659  // pressure derivatives]
660 
661  double slope1 = centre_solid_p[RDF] - centre_solid_p[LDF];
662  double slope2 = centre_solid_p[RUF] - centre_solid_p[LUF];
663  double slope3 = centre_solid_p[RDB] - centre_solid_p[LDB];
664  double slope4 = centre_solid_p[RUB] - centre_solid_p[LUB];
665 
666  // Use the average
667  internal_data_pt(this->P_solid_internal_index)->
668  set_value(1,0.25*(slope1+slope2+slope3+slope4));
669 
670 
671  //Slope in s_1 direction
672  //----------------------
673 
674  // Use average of the 4 FD approximations based on the
675  // elements central pressure values
676  // [Other options: Take average of the four
677  // pressure derivatives]
678  slope1 = centre_solid_p[LUB] - centre_solid_p[LDB];
679  slope2 = centre_solid_p[RUB] - centre_solid_p[RDB];
680  slope3 = centre_solid_p[LUF] - centre_solid_p[LDF];
681  slope4 = centre_solid_p[RUF] - centre_solid_p[RDF];
682 
683  // Use the average
684  internal_data_pt(this->P_solid_internal_index)->
685  set_value(2,0.25*(slope1+slope2+slope3+slope4));
686 
687 
688  //Slope in s_2 direction
689  //----------------------
690 
691  // Use average of the 4 FD approximations based on the
692  // elements central pressure values
693  // [Other options: Take average of the four
694  // pressure derivatives]
695  slope1 = centre_solid_p[LUF] - centre_solid_p[LUB];
696  slope2 = centre_solid_p[RUF] - centre_solid_p[RUB];
697  slope3 = centre_solid_p[LDF] - centre_solid_p[LDB];
698  slope4 = centre_solid_p[RDF] - centre_solid_p[RDB];
699 
700  // Use the average
701  internal_data_pt(this->P_solid_internal_index)->
702  set_value(3,0.25*(slope1+slope2+slope3+slope4));
703 }
704 
705 
706 //================================================================
707 /// 3D Further build: Interpolate the solid pressure values
708 //=================================================================
709 template<>
711 
712  {
714 
715  using namespace OcTreeNames;
716 
717  // What type of son am I? Ask my quadtree representation...
718  int son_type=octree_pt()->son_type();
719 
720  // Pointer to my father (in element impersonation)
721  RefineableQElement<3>* father_el_pt=
722  dynamic_cast<RefineableQElement<3>*>(
723  octree_pt()->father_pt()->object_pt());
724 
725  Vector<double> s_father(3);
726 
727  // Son midpoint is located at the following coordinates in father element:
728  for(unsigned i=0;i<3;i++)
729  {
730  s_father[i]=0.5*OcTree::Direction_to_vector[son_type][i];
731  }
732 
733  // Pressure value in father element
734  RefineableQPVDElementWithPressure<3>* cast_father_element_pt=
735  dynamic_cast<RefineableQPVDElementWithPressure<3>*>(father_el_pt);
736 
737  double press=cast_father_element_pt->interpolated_solid_p(s_father);
738 
739 
740  // Pressure value gets copied straight into internal dof:
741  set_solid_p(0,press);
742 
743  // The slopes get copied from father and halved
744  for(unsigned i=1;i<4;i++)
745  {
746  set_solid_p(i,0.5*cast_father_element_pt->solid_p(i));
747  }
748  }
749 
750 //=========================================================================
751 ///FaceGeometry of the 3D RefineableQPVDElementWithPressure
752 //========================================================================
753 template<>
755  public virtual SolidQElement<2,3>
756 {
757  public:
758  //Make sure that we call the constructor of the SolidQElement
759  //Only the Intel compiler seems to need this!
761 };
762 
763 
764 //========================================================================
765 ///FaceGeometry of the FaceGeometry of the 3D RefineableQPVDElementWithPressure
766 //==========================================================================
767 template<>
769 public virtual SolidQElement<1,3>
770 {
771  public:
772  //Make sure that we call the constructor of the SolidQElement
773  //Only the Intel compiler seems to need this!
775 };
776 
777 
778 
779 //===========================================================================
780 /// Class for refineable solid mechanics elements in near-incompressible/
781 /// incompressible formulation, so a pressure is included! These elements
782 /// include a continuously interpolated pressure a la Taylor Hood/
783 //===========================================================================
784 template<unsigned DIM>
786 public virtual QPVDElementWithContinuousPressure<DIM>,
787 public virtual RefineablePVDEquationsWithPressure<DIM>,
788 public virtual RefineableSolidQElement<DIM>
789 {
790 
791  public:
792 
793  /// \short Constructor:
799  RefineableSolidQElement<DIM>() {}
800 
801 
802  /// \short Overload the number of additional solid dofs at each node, we shall
803  /// always assign 1, otherwise it's a real pain
804  unsigned required_nvalue(const unsigned &n) const {return 1;}
805 
806  /// Number of continuously interpolated values (1) solid pressure
807  unsigned ncont_interpolated_values() const {return 1;}
808 
809  /// Empty rebuild from sons, empty
810  void rebuild_from_sons(Mesh* &mesh_pt) {}
811 
812  /// OK, interpolate the solid pressures
814  Vector<double>& values)
815  {
816  //There is only one solid pressure, initialise to zero
817  values.resize(1);
818 
819  //Get the interpolated value
820  values[0] = this->interpolated_solid_p(s);
821  }
822 
823  /// OK get the time-dependent verion
824  void get_interpolated_values(const unsigned &t, const Vector<double> &s,
825  Vector<double>& values)
826  {
827  //There is only one solid pressure, initialise to zero
828  values.resize(1);
829  //The solid pressure does not depend on time!
830  values[0] = this->interpolated_solid_p(s);
831  }
832 
833 
834  /// Unpin all pressure dofs
836  {
837  //find the index at which the pressure is stored
838  int solid_p_index = this->solid_p_nodal_index();
839  unsigned n_node = this->nnode();
840  // loop over nodes
841  for(unsigned n=0;n<n_node;n++)
842  {this->node_pt(n)->unpin(solid_p_index);}
843  }
844 
845 
846  /// Pin the redundant solid pressure
848  {
849  //Find the index of the solid pressure
850  int solid_p_index = this->solid_p_nodal_index();
851  //Let's pin all pressure nodes
852  unsigned n_node = this->nnode();
853  for(unsigned l=0;l<n_node;l++)
854  {
855  //Pin the solid pressure
856  this->node_pt(l)->pin(solid_p_index);
857  }
858 
859  //Now loop over the pressure nodes and unpin the solid pressures
860  unsigned n_solid_pres = this->npres_solid();
861  //Loop over these nodes and unpin the solid pressures
862  for(unsigned l=0;l<n_solid_pres;l++)
863  {
864  Node* nod_pt = this->solid_pressure_node_pt(l);
865  if(!nod_pt->is_hanging(solid_p_index))
866  {
867  nod_pt->unpin(solid_p_index);
868  }
869  }
870  }
871 
872 
873 /// \short Number of vertex nodes in the element
874  unsigned nvertex_node() const
876 
877  /// \short Pointer to the j-th vertex node in the element
878  Node* vertex_node_pt(const unsigned& j) const
879  {
881  }
882 
883  /// \short Order of recovery shape functions for Z2 error estimation:
884  /// Same order as shape functions.
885  unsigned nrecovery_order() {return 2;} //NNODE_1D-1;}
886 
887 
888  /// \short The pressure "nodes" are a
889  /// subset of the nodes, so when value_id==0, the n-th pressure
890  /// node is returned.
891  Node* interpolating_node_pt(const unsigned &n,
892  const int &value_id)
893 
894  {
895 #ifdef PARANOID
897 #endif
898  //If we are at the value return the solid pressure node
899  if(value_id==0)
900  {return this->solid_pressure_node_pt(n);}
901  //Otherwise return the nodal values
902  else {return this->node_pt(n);}
903  }
904 
905  /// \short The pressure nodes are the corner nodes, so when value_id==0,
906  /// the fraction is the same as the 1d node number, 0 or 1.
907  double local_one_d_fraction_of_interpolating_node(const unsigned &n1d,
908  const unsigned &i,
909  const int &value_id)
910  {
911 #ifdef PARANOID
913 #endif
914  //If it's the only value, we have the pressure
915  if(value_id==0)
916  {
917  //The pressure nodes are just located on the boundaries at 0 or 1
918  return double(n1d);
919  }
920  //Otherwise we have the geometric nodes
921  else
922  {
923  return this->local_one_d_fraction_of_node(n1d,i);
924  }
925  }
926 
927  /// \short The velocity nodes are the same as the geometric nodes. The
928  /// pressure nodes must be calculated by using the same methods as
929  /// the geometric nodes, but by recalling that there are only two pressure
930  /// nodes per edge.
932  const int &value_id)
933  {
934 #ifdef PARANOID
936 #endif
937 
938  //If we are calculating solid pressure nodes
939  if(value_id==0)
940  {
941  //Storage for the index of the pressure node
942  unsigned total_index=0;
943  //The number of nodes along each 1d edge is 2.
944  unsigned NNODE_1D = 2;
945  //Storage for the index along each boundary
946  Vector<int> index(DIM);
947  //Loop over the coordinates
948  for(unsigned i=0;i<DIM;i++)
949  {
950  //If we are at the lower limit, the index is zero
951  if(s[i]==-1.0)
952  {
953  index[i] = 0;
954  }
955  //If we are at the upper limit, the index is the number of nodes minus 1
956  else if(s[i] == 1.0)
957  {
958  index[i] = NNODE_1D-1;
959  }
960  //Otherwise, we have to calculate the index in general
961  else
962  {
963  //For uniformly spaced nodes the 0th node number would be
964  double float_index = 0.5*(1.0 + s[i])*(NNODE_1D-1);
965  index[i] = int(float_index);
966  //What is the excess. This should be safe because the
967  //taking the integer part rounds down
968  double excess = float_index - index[i];
969  //If the excess is bigger than our tolerance there is no node,
970  //return null
971  if((excess > FiniteElement::Node_location_tolerance) && (
972  (1.0 - excess) > FiniteElement::Node_location_tolerance))
973  {
974  return 0;
975  }
976  }
977  ///Construct the general pressure index from the components.
978  total_index += index[i]*
979  static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
980  static_cast<int>(i)));
981  }
982  //If we've got here we have a node, so let's return a pointer to it
983  return this->solid_pressure_node_pt(total_index);
984  }
985  //Otherwise velocity nodes are the same as pressure nodes
986  else
987  {
988  return this->get_node_at_local_coordinate(s);
989  }
990  }
991 
992 
993  /// \short The number of 1d pressure nodes is 2, otherwise we have
994  /// the positional nodes
995  unsigned ninterpolating_node_1d(const int &value_id)
996  {
997 #ifdef PARANOID
999 #endif
1000 
1001  if(value_id==0) {return 2;}
1002  else {return this->nnode_1d();}
1003  }
1004 
1005  /// \short The number of pressure nodes is 2^DIM. The number of
1006  /// velocity nodes is the same as the number of geometric nodes.
1007  unsigned ninterpolating_node(const int &value_id)
1008  {
1009 #ifdef PARANOID
1011 #endif
1012 
1013  if(value_id==0)
1014  {return static_cast<unsigned>(pow(2.0,static_cast<int>(DIM)));}
1015  else {return this->nnode();}
1016  }
1017 
1018  /// \short The basis interpolating the pressure is given by pshape().
1019  //// The basis interpolating the velocity is shape().
1021  Shape &psi,
1022  const int &value_id) const
1023  {
1024 #ifdef PARANOID
1026 #endif
1027 
1028  if(value_id==0)
1029  {return this->solid_pshape(s,psi);}
1030  else {return this->shape(s,psi);}
1031  }
1032 
1033 
1034  /// \short Perform additional hanging node procedures for variables
1035  /// that are not interpolated by all nodes.
1037  {
1038  this->setup_hang_for_value(this->solid_p_nodal_index());
1039  }
1040 
1041  //Return a pointer to the solid node at which pressure dof l2 is stored
1042  Node *solid_pressure_node_pt(const unsigned &l)
1043  {return this->node_pt(this->Pconv[l]);}
1044 
1045 };
1046 
1047 
1048 //=========================================================================
1049 ///FaceGeometry of the 2D RefineableQPVDElementWithContinuousPressure elements
1050 //=========================================================================
1051 template<>
1053  public virtual SolidQElement<1,3>
1054 {
1055  public:
1056  //Make sure that we call the constructor of the SolidQElement
1057  //Only the Intel compiler seems to need this!
1059 };
1060 
1061 //=========================================================================
1062 ///FaceGeometry of the FaceGeometry of the 2D
1063 ///RefineableQPVDElementWithContinuousPressure
1064 //=========================================================================
1065 template<>
1068  public virtual PointElement
1069 {
1070  public:
1071  //Make sure that we call the constructor of the SolidQElement
1072  //Only the Intel compiler seems to need this!
1074 };
1075 
1076 //=========================================================================
1077 ///FaceGeometry of the 3D RefineableQPVDElementWithContinuousPressure
1078 //=========================================================================
1079 template<>
1081  public virtual SolidQElement<2,3>
1082 {
1083  public:
1084  //Make sure that we call the constructor of the SolidQElement
1085  //Only the Intel compiler seems to need this!
1087 };
1088 
1089 //=========================================================================
1090 ///FaceGeometry of the FaceGeometry of the 3D
1091 ///RefineableQPVDElementWithContinuousPressue
1092 //=========================================================================
1093 template<>
1096 public virtual SolidQElement<1,3>
1097 {
1098  public:
1099  //Make sure that we call the constructor of the SolidQElement
1100  //Only the Intel compiler seems to need this!
1102 };
1103 
1104 }
1105 
1106 #endif
void fill_in_generic_residual_contribution_pvd_with_pressure(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
unsigned nvertex_node() const
Number of vertex nodes in the element.
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:386
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain tensor.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
No values are interpolated in this element (pure solid)
void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
The basis interpolating the pressure is given by pshape(). / The basis interpolating the velocity is ...
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values (1) solid pressure.
bool is_jacobian_evaluated_by_fd() const
Return the flag indicating whether the jacobian is evaluated by fd.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
OK get the time-dependent verion.
void unpin_elemental_solid_pressure_dofs()
Unpin all solid pressure dofs.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
cstr elem_len * i
Definition: cfortran.h:607
void pin_elemental_redundant_nodal_solid_pressures()
Pin the redundant solid pressure.
int solid_p_nodal_index() const
Set the value at which the solid pressure is stored in the nodes.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
No values are interpolated in this element (pure solid)
unsigned npres_solid() const
Return number of pressure values.
static void check_value_id(const int &n_continuously_interpolated_values, const int &value_id)
Static helper function that is used to check that the value_id is in range.
Node * interpolating_node_pt(const unsigned &n, const int &value_id)
The pressure "nodes" are a subset of the nodes, so when value_id==0, the n-th pressure node is return...
void further_setup_hanging_nodes()
No additional hanging node procedures are required for the solid elements.
ConstitutiveLaw * Constitutive_law_pt
Pointer to the constitutive law.
virtual Node * solid_pressure_node_pt(const unsigned &l)
char t
Definition: cfortran.h:572
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void rebuild_from_sons(Mesh *&mesh_pt)
Empty rebuild from sons, no need to reconstruct anything here.
void solid_pshape(const Vector< double > &s, Shape &psi) const
Return the pressure shape functions.
double *& lambda_sq_pt()
Access function for pointer to timescale ratio (nondim density)
void rebuild_from_sons(Mesh *&mesh_pt)
Empty rebuild from sons, empty.
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2139
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
unsigned npres_solid() const
Return number of pressure values.
bool Unsteady
Flag that switches inertia on/off.
unsigned nvertex_node() const
Number of vertex nodes in the element.
Class for Refineable PVD equations.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void further_setup_hanging_nodes()
No additional hanging node procedures are required for discontinuous solid pressures.
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 0 (pure solid problem)
void get_mass_matrix_diagonal(Vector< double > &mass_diag)
Compute the diagonal of the displacement mass matrix for LSC preconditioner.
void get_strain(const Vector< double > &s, DenseMatrix< double > &strain) const
Return the strain tensor.
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:342
static char t char * s
Definition: cfortran.h:572
unsigned required_nvalue(const unsigned &n) const
Overload the number of additional solid dofs at each node, we shall always assign 1...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
IsotropicGrowthFctPt & isotropic_growth_fct_pt()
Access function: Pointer to isotropic growth function.
bool is_inertia_enabled() const
Access function to flag that switches inertia on/off (const version)
virtual void further_build()
Further build: Pass the father's Use_undeformed_macro_element_for_new_lagrangian_coords flag down...
void rebuild_from_sons(Mesh *&mesh_pt)
Reconstruct the pressure from the sons Must be specialized for each dimension.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
in strain tensor.
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3804
Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
The velocity nodes are the same as the geometric nodes. The pressure nodes must be calculated by usin...
unsigned ninterpolating_node(const int &value_id)
The number of pressure nodes is 2^DIM. The number of velocity nodes is the same as the number of geom...
double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &value_id)
The pressure nodes are the corner nodes, so when value_id==0, the fraction is the same as the 1d node...
bool Incompressible
Boolean to determine whether the solid is incompressible or not.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1334
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
bool Evaluate_jacobian_by_fd
Use FD to evaluate Jacobian.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
bool is_incompressible() const
Return whether the material is incompressible.
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements...
Definition: elements.h:2371
IsotropicGrowthFctPt Isotropic_growth_fct_pt
Pointer to isotropic growth function.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
OK, interpolate the solid pressures.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
No values are interpolated in this element (pure solid)
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 0 (pure solid problem)
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
Class for refineable QPVDElement elements.
void further_build()
Pass the generic stuff down to the sons.
void unpin_elemental_solid_pressure_dofs()
Unpin all pressure dofs.
virtual unsigned nvertex_node() const
Definition: elements.h:2362
ConstitutiveLaw *& constitutive_law_pt()
Return the constitutive law pointer.
BodyForceFctPt Body_force_fct_pt
Pointer to body force function.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
No values are interpolated in this element (pure solid)
unsigned P_solid_internal_index
Internal index that indicates at which internal data value the solid presure is stored.
virtual Node * solid_pressure_node_pt(const unsigned &l)
unsigned nvertex_node() const
Number of vertex nodes in the element.
void further_build()
Further build: Interpolate the solid pressure values Again this must be specialised for each dimensio...
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Get the local fraction of any node in the n-th position in a one dimensional expansion along the i-th...
Definition: elements.h:1810
double interpolated_solid_p(const Vector< double > &s)
Return the interpolated_solid_pressure.
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 0 (pure solid problem)
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
A general mesh class.
Definition: mesh.h:74
void fill_in_generic_contribution_to_residuals_pvd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Call the residuals including hanging node cases.
void further_build()
Further build function, pass the pointers down to the sons.
unsigned ninterpolating_node_1d(const int &value_id)
The number of 1d pressure nodes is 2, otherwise we have the positional nodes.