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// 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 refineable solid mechanics elements
27
28// Include guard to prevent multiple inclusions of this header
29#ifndef OOMPH_REFINEABLE_ELASTICITY_ELEMENTS_HEADER
30#define OOMPH_REFINEABLE_ELASTICITY_ELEMENTS_HEADER
31
32// oomph-lib headers
33#include "solid_elements.h"
34#include "../generic/refineable_quad_element.h"
35#include "../generic/refineable_brick_element.h"
36#include "../generic/error_estimator.h"
37
38namespace oomph
39{
40 //========================================================================
41 /// Class for Refineable PVD equations
42 //========================================================================
43 template<unsigned DIM>
44 class RefineablePVDEquations : public virtual PVDEquations<DIM>,
45 public virtual RefineableSolidElement,
46 public virtual ElementWithZ2ErrorEstimator
47 {
48 public:
49 /// Constructor
51 : PVDEquations<DIM>(),
55 {
56 }
57
58 /// Call the residuals including hanging node cases
60 Vector<double>& residuals,
61 DenseMatrix<double>& jacobian,
62 const unsigned& flag);
63
64 /// No values are interpolated in this element (pure solid)
65 void get_interpolated_values(const unsigned& t,
66 const Vector<double>& s,
67 Vector<double>& values)
68 {
69 values.clear();
70 }
71
72 /// No values are interpolated in this element (pure solid)
74 Vector<double>& values)
75 {
76 values.clear();
77 }
78
79 /// Number of 'flux' terms for Z2 error estimation
81 {
82 // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
83 return DIM + DIM * (DIM - 1) / 2;
84 }
85
86 /// Get 'flux' for Z2 error recovery: Upper triangular entries
87 /// in strain tensor.
89 {
90#ifdef PARANOID
91 unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
92 if (flux.size() != num_entries)
93 {
94 std::ostringstream error_message;
95 error_message << "The flux vector has the wrong number of entries, "
96 << flux.size() << ", whereas it should be " << num_entries
97 << std::endl;
98 throw OomphLibError(error_message.str(),
99 OOMPH_CURRENT_FUNCTION,
100 OOMPH_EXCEPTION_LOCATION);
101 }
102#endif
103
104 // Get strain matrix
105 DenseMatrix<double> strain(DIM);
106 this->get_strain(s, strain);
107
108 // Pack into flux Vector
109 unsigned icount = 0;
110
111 // Start with diagonal terms
112 for (unsigned i = 0; i < DIM; i++)
113 {
114 flux[icount] = strain(i, i);
115 icount++;
116 }
117
118 // Off diagonals row by row
119 for (unsigned i = 0; i < DIM; i++)
120 {
121 for (unsigned j = i + 1; j < DIM; j++)
122 {
123 flux[icount] = strain(i, j);
124 icount++;
125 }
126 }
127 }
128
129 /// Number of continuously interpolated values: 0 (pure solid problem)
131 {
132 return 0;
133 }
134
135 // Return a pointer to the solid node at which pressure dof l2 is stored
136 // This is only required so that the generic templating in
137 // PseudoSolidNodeUpdateElements works OK
138 virtual Node* solid_pressure_node_pt(const unsigned& l)
139 {
140 return 0;
141 }
142
143 /// Further build function, pass the pointers down to the sons
145 {
146 RefineablePVDEquations<DIM>* cast_father_element_pt =
147 dynamic_cast<RefineablePVDEquations<DIM>*>(this->father_element_pt());
148
149 // Do whatever needs to be done in the base class
151
152 // Set pointer to isotropic growth function
154 cast_father_element_pt->isotropic_growth_fct_pt();
155
156 // Set pointer to body force function
157 this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
158
159 // Set pointer to the contitutive law
160 this->Constitutive_law_pt = cast_father_element_pt->constitutive_law_pt();
161
162 // Set the timescale ratio (non-dim. density)
163 this->Lambda_sq_pt = cast_father_element_pt->lambda_sq_pt();
164
165 // Set the flag that switches inertia on/off
166 this->Unsteady = cast_father_element_pt->is_inertia_enabled();
167
168 // Evaluation of Jacobian by same method as father
170 cast_father_element_pt->is_jacobian_evaluated_by_fd();
171 }
172 };
173
174 //========================================================================
175 /// Class for refineable QPVDElement elements
176 //========================================================================
177 template<unsigned DIM, unsigned NNODE_1D>
178 class RefineableQPVDElement : public virtual QPVDElement<DIM, NNODE_1D>,
179 public virtual RefineablePVDEquations<DIM>,
180 public virtual RefineableSolidQElement<DIM>
181 {
182 public:
183 /// Constructor:
185 : QPVDElement<DIM, NNODE_1D>(),
190 {
191 }
192
193 /// Empty rebuild from sons, no need to reconstruct anything here
194 void rebuild_from_sons(Mesh*& mesh_pt) {}
195
196 /// Number of vertex nodes in the element
197 unsigned nvertex_node() const
198 {
200 }
201
202 /// Pointer to the j-th vertex node in the element
203 Node* vertex_node_pt(const unsigned& j) const
204 {
206 }
207
208 /// Order of recovery shape functions for Z2 error estimation:
209 /// Same order as shape functions.
211 {
212 return NNODE_1D - 1;
213 }
214
215 /// No additional hanging node procedures are required
216 /// for the solid elements.
218 };
219
220 //==============================================================
221 /// FaceGeometry of the 2D RefineableQPVDElement elements
222 //==============================================================
223 template<unsigned NNODE_1D>
225 : public virtual SolidQElement<1, NNODE_1D>
226 {
227 public:
228 // Make sure that we call the constructor of the SolidQElement
229 // Only the Intel compiler seems to need this!
230 FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
231 };
232
233 //==============================================================
234 /// FaceGeometry of the FaceGeometry of the 2D RefineableQPVDElement
235 //==============================================================
236 template<unsigned NNODE_1D>
238 : public virtual PointElement
239 {
240 public:
241 // Make sure that we call the constructor of the SolidQElement
242 // Only the Intel compiler seems to need this!
244 };
245
246
247 //==============================================================
248 /// FaceGeometry of the 3D RefineableQPVDElement elements
249 //==============================================================
250 template<unsigned NNODE_1D>
252 : public virtual SolidQElement<2, 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<2, NNODE_1D>() {}
258 };
259
260 //==============================================================
261 /// FaceGeometry of the FaceGeometry of the 3D RefineableQPVDElement
262 //==============================================================
263 template<unsigned NNODE_1D>
265 : public virtual SolidQElement<1, NNODE_1D>
266 {
267 public:
268 // Make sure that we call the constructor of the SolidQElement
269 // Only the Intel compiler seems to need this!
270 FaceGeometry() : SolidQElement<1, NNODE_1D>() {}
271 };
272
273
274 //===========================================================================
275 /// Class for Refineable solid mechanics elements in near-incompressible/
276 /// incompressible formulation, so a pressure is included! In this case,
277 /// the pressure interpolation is discontinuous, a la Crouzeix Raviart
278 //===========================================================================
279 template<unsigned DIM>
281 : public virtual PVDEquationsWithPressure<DIM>,
282 public virtual RefineableSolidElement,
283 public virtual ElementWithZ2ErrorEstimator
284 {
285 public:
286 /// Constructor:
292 {
293 }
294
295 /// Add element's contribution to elemental residual vector and/or
296 /// Jacobian matrix
297 /// flag=1: compute both
298 /// flag=0: compute only residual vector
300 Vector<double>& residuals,
301 DenseMatrix<double>& jacobian,
302 DenseMatrix<double>& mass_matrix,
303 const unsigned& flag);
304
305 /// No values are interpolated in this element (pure solid)
306 void get_interpolated_values(const unsigned& t,
307 const Vector<double>& s,
308 Vector<double>& values)
309 {
310 values.clear();
311 }
312
313 /// No values are interpolated in this element (pure solid)
315 Vector<double>& values)
316 {
317 values.clear();
318 }
319
320 /// Number of 'flux' terms for Z2 error estimation
322 {
323 // DIM Diagonal strain rates and DIM*(DIM-1)/2 off diagonal terms
324 return DIM + DIM * (DIM - 1) / 2;
325 }
326
327 // Get 'flux' for Z2 error recovery: Upper triangular entries
328 /// in strain tensor.
329 //----------------------------------------------------------------
331 {
332 // Find the dimension of the problem
333#ifdef PARANOID
334 unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
335 if (flux.size() != num_entries)
336 {
337 std::ostringstream error_message;
338 error_message << "The flux vector has the wrong number of entries, "
339 << flux.size() << ", whereas it should be " << num_entries
340 << std::endl;
341 throw OomphLibError(error_message.str(),
342 OOMPH_CURRENT_FUNCTION,
343 OOMPH_EXCEPTION_LOCATION);
344 }
345#endif
346
347 // Get strain matrix
348 DenseMatrix<double> strain(DIM);
349 this->get_strain(s, strain);
350
351 // Pack into flux Vector
352 unsigned icount = 0;
353
354 // Start with diagonal terms
355 for (unsigned i = 0; i < DIM; i++)
356 {
357 flux[icount] = strain(i, i);
358 icount++;
359 }
360
361 // Off diagonals row by row
362 for (unsigned i = 0; i < DIM; i++)
363 {
364 for (unsigned j = i + 1; j < DIM; j++)
365 {
366 flux[icount] = strain(i, j);
367 icount++;
368 }
369 }
370 }
371
372 /// Number of continuously interpolated values: 0 (pure solid problem)
374 {
375 return 0;
376 }
377
378 // Return a pointer to the solid node at which pressure dof l2 is stored
379 virtual Node* solid_pressure_node_pt(const unsigned& l)
380 {
381 return 0;
382 }
383
384
385 /// Pass the generic stuff down to the sons
387 {
388 RefineablePVDEquationsWithPressure<DIM>* cast_father_element_pt =
390 this->father_element_pt());
391
392 // Do whatever needs to be done in the base class
394
395 // Set pointer to isotropic growth function
397 cast_father_element_pt->isotropic_growth_fct_pt();
398
399 // Set pointer to body force function
400 this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
401
402 // Set pointer to the contitutive law
403 this->Constitutive_law_pt = cast_father_element_pt->constitutive_law_pt();
404
405 // Set the timescale ratio (non-dim. density)
406 this->Lambda_sq_pt = cast_father_element_pt->lambda_sq_pt();
407
408 // Set the flag that switches inertia on/off
409 this->Unsteady = cast_father_element_pt->is_inertia_enabled();
410
411 // Set the incompressibility flag
412 this->Incompressible = cast_father_element_pt->is_incompressible();
413
414 // Evaluation of Jacobian by same method as father
416 cast_father_element_pt->is_jacobian_evaluated_by_fd();
417 }
418
419
420 /// Compute the diagonal of the displacement mass matrix for
421 /// LSC preconditioner
423 };
424
425 //===========================================================================
426 /// Class for refineable solid mechanics elements in near-incompressible/
427 /// incompressible formulation, so a pressure is included! In this case,
428 /// the pressure interpolation is discontinuous, a la Crouzeix Raviart,
429 /// and the displacement is always quadratic.
430 //===========================================================================
431 template<unsigned DIM>
433 : public virtual QPVDElementWithPressure<DIM>,
434 public virtual RefineablePVDEquationsWithPressure<DIM>,
435 public virtual RefineableSolidQElement<DIM>
436 {
437 private:
438 /// Unpin all solid pressure dofs
440 {
441 unsigned n_pres = this->npres_solid();
442 // loop over pressure dofs and unpin them
443 for (unsigned l = 0; l < n_pres; l++)
444 {
446 }
447 }
448
449 public:
450 /// Constructor:
457 {
458 }
459
460
461 /// Reconstruct the pressure from the sons
462 /// Must be specialized for each dimension
463 inline void rebuild_from_sons(Mesh*& mesh_pt);
464
465 /// Number of vertex nodes in the element
466 unsigned nvertex_node() const
467 {
469 }
470
471 /// Pointer to the j-th vertex node in the element
472 Node* vertex_node_pt(const unsigned& j) const
473 {
475 }
476
477 /// Order of recovery shape functions for Z2 error estimation:
478 /// Same order as shape functions.
480 {
481 return 2;
482 } // NNODE_1D-1;}
483
484 /// No additional hanging node procedures are required for
485 /// discontinuous solid pressures.
487
488 /// Further build: Interpolate the solid pressure values
489 /// Again this must be specialised for each dimension
490 inline void further_build();
491
492
493 /// Number of continuously interpolated values: 0 (pure solid problem)
495 {
496 return 0;
497 }
498 };
499
500 //======================================================================
501 /// FaceGeometry of the 2D RefineableQPVDElementWithPressure
502 //=======================================================================
503 template<>
505 : public virtual SolidQElement<1, 3>
506 {
507 public:
508 // Make sure that we call the constructor of the SolidQElement
509 // Only the Intel compiler seems to need this!
511 };
512
513
514 //===========================================================================
515 /// FaceGeometry of the FaceGeometry of the 2D
516 /// RefineableQPVDElementWithPressure
517 //============================================================================
518 template<>
520 : public virtual PointElement
521 {
522 public:
523 // Make sure that we call the constructor of the SolidQElement
524 // Only the Intel compiler seems to need this!
526 };
527
528 //====================================================================
529 /// 2D rebuild from sons: reconstruct the solid pressure
530 //===================================================================
531 template<>
533 Mesh*& mesh_pt)
534 {
535 using namespace QuadTreeNames;
536
537 // Storage for the solid pressure of each son
538 double centre_solid_p[4] = {0.0, 0.0, 0.0, 0.0};
539
540 // Loop over the sons and assign the central solid pressure
541 for (unsigned ison = 0; ison < 4; ison++)
542 {
543 // Add the sons midnode pressure
544 centre_solid_p[ison] =
546 quadtree_pt()->son_pt(ison)->object_pt())
547 ->solid_p(0);
548 }
549
550 // Use the average for the central solid pressure
551 double p_value = 0.25 * (centre_solid_p[0] + centre_solid_p[1] +
552 centre_solid_p[2] + centre_solid_p[3]);
553 // Actually set the pressure
554 set_solid_p(0, p_value);
555
556
557 // Slope in s_0 direction
558 //----------------------
559
560 // Use average of the 2 FD approximations based on the
561 // elements central pressure values
562 // [Other options: Take average of the four
563 // pressure derivatives]
564 double slope1 = centre_solid_p[SE] - centre_solid_p[SW];
565 double slope2 = centre_solid_p[NE] - centre_solid_p[NW];
566
567
568 // Use the average value
569 p_value = 0.5 * (slope1 + slope2);
570 set_solid_p(1, p_value);
571
572 // Slope in s_1 direction
573 //----------------------
574
575 // Use average of the 2 FD approximations based on the
576 // elements central pressure values
577 // [Other options: Take average of the four
578 // pressure derivatives]
579
580 slope1 = centre_solid_p[NE] - centre_solid_p[SE];
581 slope2 = centre_solid_p[NW] - centre_solid_p[SW];
582
583 // Use the average
584 p_value = 0.5 * (slope1 + slope2);
585 set_solid_p(2, p_value);
586 }
587
588
589 //==============================================================
590 /// 2D further build interpolates the internal solid pressure
591 /// from the father.
592 //=============================================================
593 template<>
595 {
597
598 using namespace QuadTreeNames;
599
600 // What type of son am I? Ask my quadtree representation...
601 int son_type = quadtree_pt()->son_type();
602
603 // Pointer to my father (in element impersonation)
604 RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
605
606 Vector<double> s_father(2);
607
608 // Son midpoint is located at the following coordinates in father element:
609
610 // South west son
611 if (son_type == SW)
612 {
613 s_father[0] = -0.5;
614 s_father[1] = -0.5;
615 }
616 // South east son
617 else if (son_type == SE)
618 {
619 s_father[0] = 0.5;
620 s_father[1] = -0.5;
621 }
622 // North east son
623 else if (son_type == NE)
624 {
625 s_father[0] = 0.5;
626 s_father[1] = 0.5;
627 }
628
629 // North west son
630 else if (son_type == NW)
631 {
632 s_father[0] = -0.5;
633 s_father[1] = 0.5;
634 }
635
636 // Pressure value in father element
637 RefineableQPVDElementWithPressure<2>* cast_father_element_pt =
638 dynamic_cast<RefineableQPVDElementWithPressure<2>*>(father_el_pt);
639 double press = cast_father_element_pt->interpolated_solid_p(s_father);
640
641 // Pressure value gets copied straight into internal dof:
642 set_solid_p(0, press);
643
644 // The slopes get copied from father and halved
645 set_solid_p(1, 0.5 * cast_father_element_pt->solid_p(1));
646 set_solid_p(2, 0.5 * cast_father_element_pt->solid_p(2));
647 }
648
649
650 //===============================================================
651 /// 3D rebuild from sons: reconstruct the pressure
652 //===============================================================
653 template<>
655 Mesh*& mesh_pt)
656 {
657 using namespace OcTreeNames;
658
659 // Storage for the central solid pressure of each son
660 double centre_solid_p[8] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
661
662 // Loop over the sons and assign the central solid pressure
663 for (unsigned ison = 0; ison < 8; ison++)
664 {
665 centre_solid_p[ison] = octree_pt()
666 ->son_pt(ison)
667 ->object_pt()
668 ->internal_data_pt(this->P_solid_internal_index)
669 ->value(0);
670 }
671
672 // Central pressure value:
673 //-----------------------
674
675 // Use average of the sons central pressure values
676 // Other options: Take average of the four (discontinuous)
677 // pressure values at the father's midpoint]
678 double av_press = 0.0;
679
680 // Loop over the sons and sum the centre pressures
681 for (unsigned ison = 0; ison < 8; ison++)
682 {
683 av_press += centre_solid_p[ison];
684 }
685
686 // Use the average
687 internal_data_pt(this->P_solid_internal_index)
688 ->set_value(0, 0.125 * av_press);
689
690
691 // Slope in s_0 direction
692 //----------------------
693
694 // Use average of the 4 FD approximations based on the
695 // elements central pressure values
696 // [Other options: Take average of the four
697 // pressure derivatives]
698
699 double slope1 = centre_solid_p[RDF] - centre_solid_p[LDF];
700 double slope2 = centre_solid_p[RUF] - centre_solid_p[LUF];
701 double slope3 = centre_solid_p[RDB] - centre_solid_p[LDB];
702 double slope4 = centre_solid_p[RUB] - centre_solid_p[LUB];
703
704 // Use the average
705 internal_data_pt(this->P_solid_internal_index)
706 ->set_value(1, 0.25 * (slope1 + slope2 + slope3 + slope4));
707
708
709 // Slope in s_1 direction
710 //----------------------
711
712 // Use average of the 4 FD approximations based on the
713 // elements central pressure values
714 // [Other options: Take average of the four
715 // pressure derivatives]
716 slope1 = centre_solid_p[LUB] - centre_solid_p[LDB];
717 slope2 = centre_solid_p[RUB] - centre_solid_p[RDB];
718 slope3 = centre_solid_p[LUF] - centre_solid_p[LDF];
719 slope4 = centre_solid_p[RUF] - centre_solid_p[RDF];
720
721 // Use the average
722 internal_data_pt(this->P_solid_internal_index)
723 ->set_value(2, 0.25 * (slope1 + slope2 + slope3 + slope4));
724
725
726 // Slope in s_2 direction
727 //----------------------
728
729 // Use average of the 4 FD approximations based on the
730 // elements central pressure values
731 // [Other options: Take average of the four
732 // pressure derivatives]
733 slope1 = centre_solid_p[LUF] - centre_solid_p[LUB];
734 slope2 = centre_solid_p[RUF] - centre_solid_p[RUB];
735 slope3 = centre_solid_p[LDF] - centre_solid_p[LDB];
736 slope4 = centre_solid_p[RDF] - centre_solid_p[RDB];
737
738 // Use the average
739 internal_data_pt(this->P_solid_internal_index)
740 ->set_value(3, 0.25 * (slope1 + slope2 + slope3 + slope4));
741 }
742
743
744 //================================================================
745 /// 3D Further build: Interpolate the solid pressure values
746 //=================================================================
747 template<>
749
750 {
752
753 using namespace OcTreeNames;
754
755 // What type of son am I? Ask my quadtree representation...
756 int son_type = octree_pt()->son_type();
757
758 // Pointer to my father (in element impersonation)
759 RefineableQElement<3>* father_el_pt = dynamic_cast<RefineableQElement<3>*>(
760 octree_pt()->father_pt()->object_pt());
761
762 Vector<double> s_father(3);
763
764 // Son midpoint is located at the following coordinates in father element:
765 for (unsigned i = 0; i < 3; i++)
766 {
767 s_father[i] = 0.5 * OcTree::Direction_to_vector[son_type][i];
768 }
769
770 // Pressure value in father element
771 RefineableQPVDElementWithPressure<3>* cast_father_element_pt =
772 dynamic_cast<RefineableQPVDElementWithPressure<3>*>(father_el_pt);
773
774 double press = cast_father_element_pt->interpolated_solid_p(s_father);
775
776
777 // Pressure value gets copied straight into internal dof:
778 set_solid_p(0, press);
779
780 // The slopes get copied from father and halved
781 for (unsigned i = 1; i < 4; i++)
782 {
783 set_solid_p(i, 0.5 * cast_father_element_pt->solid_p(i));
784 }
785 }
786
787 //=========================================================================
788 /// FaceGeometry of the 3D RefineableQPVDElementWithPressure
789 //========================================================================
790 template<>
792 : public virtual SolidQElement<2, 3>
793 {
794 public:
795 // Make sure that we call the constructor of the SolidQElement
796 // Only the Intel compiler seems to need this!
798 };
799
800
801 //========================================================================
802 /// FaceGeometry of the FaceGeometry of the 3D
803 /// RefineableQPVDElementWithPressure
804 //==========================================================================
805 template<>
807 : public virtual SolidQElement<1, 3>
808 {
809 public:
810 // Make sure that we call the constructor of the SolidQElement
811 // Only the Intel compiler seems to need this!
813 };
814
815
816 //===========================================================================
817 /// Class for refineable solid mechanics elements in near-incompressible/
818 /// incompressible formulation, so a pressure is included! These elements
819 /// include a continuously interpolated pressure a la Taylor Hood/
820 //===========================================================================
821 template<unsigned DIM>
823 : public virtual QPVDElementWithContinuousPressure<DIM>,
824 public virtual RefineablePVDEquationsWithPressure<DIM>,
825 public virtual RefineableSolidQElement<DIM>
826 {
827 public:
828 /// Constructor:
835 {
836 }
837
838
839 /// Overload the number of additional solid dofs at each node, we
840 /// shall always assign 1, otherwise it's a real pain
841 unsigned required_nvalue(const unsigned& n) const
842 {
843 return 1;
844 }
845
846 /// Number of continuously interpolated values (1) solid pressure
848 {
849 return 1;
850 }
851
852 /// Empty rebuild from sons, empty
853 void rebuild_from_sons(Mesh*& mesh_pt) {}
854
855 /// OK, interpolate the solid pressures
857 Vector<double>& values)
858 {
859 // There is only one solid pressure, initialise to zero
860 values.resize(1);
861
862 // Get the interpolated value
863 values[0] = this->interpolated_solid_p(s);
864 }
865
866 /// OK get the time-dependent verion
867 void get_interpolated_values(const unsigned& t,
868 const Vector<double>& s,
869 Vector<double>& values)
870 {
871 // There is only one solid pressure, initialise to zero
872 values.resize(1);
873 // The solid pressure does not depend on time!
874 values[0] = this->interpolated_solid_p(s);
875 }
876
877
878 /// Unpin all pressure dofs
880 {
881 // find the index at which the pressure is stored
882 int solid_p_index = this->solid_p_nodal_index();
883 unsigned n_node = this->nnode();
884 // loop over nodes
885 for (unsigned n = 0; n < n_node; n++)
886 {
887 this->node_pt(n)->unpin(solid_p_index);
888 }
889 }
890
891
892 /// Pin the redundant solid pressure
894 {
895 // Find the index of the solid pressure
896 int solid_p_index = this->solid_p_nodal_index();
897 // Let's pin all pressure nodes
898 unsigned n_node = this->nnode();
899 for (unsigned l = 0; l < n_node; l++)
900 {
901 // Pin the solid pressure
902 this->node_pt(l)->pin(solid_p_index);
903 }
904
905 // Now loop over the pressure nodes and unpin the solid pressures
906 unsigned n_solid_pres = this->npres_solid();
907 // Loop over these nodes and unpin the solid pressures
908 for (unsigned l = 0; l < n_solid_pres; l++)
909 {
910 Node* nod_pt = this->solid_pressure_node_pt(l);
911 if (!nod_pt->is_hanging(solid_p_index))
912 {
913 nod_pt->unpin(solid_p_index);
914 }
915 }
916 }
917
918
919 /// Number of vertex nodes in the element
920 unsigned nvertex_node() const
921 {
923 }
924
925 /// Pointer to the j-th vertex node in the element
926 Node* vertex_node_pt(const unsigned& j) const
927 {
929 }
930
931 /// Order of recovery shape functions for Z2 error estimation:
932 /// Same order as shape functions.
934 {
935 return 2;
936 } // NNODE_1D-1;}
937
938
939 /// The pressure "nodes" are a
940 /// subset of the nodes, so when value_id==0, the n-th pressure
941 /// node is returned.
942 Node* interpolating_node_pt(const unsigned& n, const int& value_id)
943
944 {
945#ifdef PARANOID
947#endif
948 // If we are at the value return the solid pressure node
949 if (value_id == 0)
950 {
951 return this->solid_pressure_node_pt(n);
952 }
953 // Otherwise return the nodal values
954 else
955 {
956 return this->node_pt(n);
957 }
958 }
959
960 /// The pressure nodes are the corner nodes, so when value_id==0,
961 /// the fraction is the same as the 1d node number, 0 or 1.
963 const unsigned& i,
964 const int& value_id)
965 {
966#ifdef PARANOID
968#endif
969 // If it's the only value, we have the pressure
970 if (value_id == 0)
971 {
972 // The pressure nodes are just located on the boundaries at 0 or 1
973 return double(n1d);
974 }
975 // Otherwise we have the geometric nodes
976 else
977 {
978 return this->local_one_d_fraction_of_node(n1d, i);
979 }
980 }
981
982 /// The velocity nodes are the same as the geometric nodes. The
983 /// pressure nodes must be calculated by using the same methods as
984 /// the geometric nodes, but by recalling that there are only two pressure
985 /// nodes per edge.
987 const int& value_id)
988 {
989#ifdef PARANOID
991#endif
992
993 // If we are calculating solid pressure nodes
994 if (value_id == 0)
995 {
996 // Storage for the index of the pressure node
997 unsigned total_index = 0;
998 // The number of nodes along each 1d edge is 2.
999 unsigned NNODE_1D = 2;
1000 // Storage for the index along each boundary
1001 Vector<int> index(DIM);
1002 // Loop over the coordinates
1003 for (unsigned i = 0; i < DIM; i++)
1004 {
1005 // If we are at the lower limit, the index is zero
1006 if (s[i] == -1.0)
1007 {
1008 index[i] = 0;
1009 }
1010 // If we are at the upper limit, the index is the number of nodes
1011 // minus 1
1012 else if (s[i] == 1.0)
1013 {
1014 index[i] = NNODE_1D - 1;
1015 }
1016 // Otherwise, we have to calculate the index in general
1017 else
1018 {
1019 // For uniformly spaced nodes the 0th node number would be
1020 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
1021 index[i] = int(float_index);
1022 // What is the excess. This should be safe because the
1023 // taking the integer part rounds down
1024 double excess = float_index - index[i];
1025 // If the excess is bigger than our tolerance there is no node,
1026 // return null
1028 ((1.0 - excess) > FiniteElement::Node_location_tolerance))
1029 {
1030 return 0;
1031 }
1032 }
1033 /// Construct the general pressure index from the components.
1034 total_index +=
1035 index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
1036 static_cast<int>(i)));
1037 }
1038 // If we've got here we have a node, so let's return a pointer to it
1039 return this->solid_pressure_node_pt(total_index);
1040 }
1041 // Otherwise velocity nodes are the same as pressure nodes
1042 else
1043 {
1044 return this->get_node_at_local_coordinate(s);
1045 }
1046 }
1047
1048
1049 /// The number of 1d pressure nodes is 2, otherwise we have
1050 /// the positional nodes
1051 unsigned ninterpolating_node_1d(const int& value_id)
1052 {
1053#ifdef PARANOID
1055#endif
1056
1057 if (value_id == 0)
1058 {
1059 return 2;
1060 }
1061 else
1062 {
1063 return this->nnode_1d();
1064 }
1065 }
1066
1067 /// The number of pressure nodes is 2^DIM. The number of
1068 /// velocity nodes is the same as the number of geometric nodes.
1069 unsigned ninterpolating_node(const int& value_id)
1070 {
1071#ifdef PARANOID
1073#endif
1074
1075 if (value_id == 0)
1076 {
1077 return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM)));
1078 }
1079 else
1080 {
1081 return this->nnode();
1082 }
1083 }
1084
1085 /// The basis interpolating the pressure is given by pshape().
1086 /// / The basis interpolating the velocity is shape().
1088 Shape& psi,
1089 const int& value_id) const
1090 {
1091#ifdef PARANOID
1093#endif
1094
1095 if (value_id == 0)
1096 {
1097 return this->solid_pshape(s, psi);
1098 }
1099 else
1100 {
1101 return this->shape(s, psi);
1102 }
1103 }
1104
1105
1106 /// Perform additional hanging node procedures for variables
1107 /// that are not interpolated by all nodes.
1109 {
1110 this->setup_hang_for_value(this->solid_p_nodal_index());
1111 }
1112
1113 // Return a pointer to the solid node at which pressure dof l2 is stored
1114 Node* solid_pressure_node_pt(const unsigned& l)
1115 {
1116 return this->node_pt(this->Pconv[l]);
1117 }
1118 };
1119
1120
1121 //=========================================================================
1122 /// FaceGeometry of the 2D RefineableQPVDElementWithContinuousPressure
1123 /// elements
1124 //=========================================================================
1125 template<>
1127 : public virtual SolidQElement<1, 3>
1128 {
1129 public:
1130 // Make sure that we call the constructor of the SolidQElement
1131 // Only the Intel compiler seems to need this!
1133 };
1134
1135 //=========================================================================
1136 /// FaceGeometry of the FaceGeometry of the 2D
1137 /// RefineableQPVDElementWithContinuousPressure
1138 //=========================================================================
1139 template<>
1142 : public virtual PointElement
1143 {
1144 public:
1145 // Make sure that we call the constructor of the SolidQElement
1146 // Only the Intel compiler seems to need this!
1148 };
1149
1150 //=========================================================================
1151 /// FaceGeometry of the 3D RefineableQPVDElementWithContinuousPressure
1152 //=========================================================================
1153 template<>
1155 : public virtual SolidQElement<2, 3>
1156 {
1157 public:
1158 // Make sure that we call the constructor of the SolidQElement
1159 // Only the Intel compiler seems to need this!
1161 };
1162
1163 //=========================================================================
1164 /// FaceGeometry of the FaceGeometry of the 3D
1165 /// RefineableQPVDElementWithContinuousPressue
1166 //=========================================================================
1167 template<>
1170 : public virtual SolidQElement<1, 3>
1171 {
1172 public:
1173 // Make sure that we call the constructor of the SolidQElement
1174 // Only the Intel compiler seems to need this!
1176 };
1177
1178} // namespace oomph
1179
1180#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
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
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
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:3882
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1374
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...
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
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:2218
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:2500
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:1858
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
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:353
An OomphLibError object which should be thrown when an run-time error is encountered....
bool Unsteady
Flag that switches inertia on/off.
ConstitutiveLaw *& constitutive_law_pt()
Return the constitutive law pointer.
BodyForceFctPt & body_force_fct_pt()
Access function: Pointer to body force function.
BodyForceFctPt Body_force_fct_pt
Pointer to body force function.
bool is_jacobian_evaluated_by_fd() const
Return the flag indicating whether the jacobian is evaluated by fd.
ConstitutiveLaw * Constitutive_law_pt
Pointer to the constitutive law.
IsotropicGrowthFctPt & isotropic_growth_fct_pt()
Access function: Pointer to isotropic growth function.
bool is_inertia_enabled() const
Access function to flag that switches inertia on/off (const version)
void 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)
bool Evaluate_jacobian_by_fd
Use FD to evaluate Jacobian.
double * Lambda_sq_pt
Timescale ratio (non-dim. density)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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.
bool is_incompressible() const
Return whether the material is incompressible.
A class for elements that solve the equations of solid mechanics, based on the principle of virtual d...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
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_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.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned npres_solid() const
Return number of pressure values.
double solid_p(const unsigned &l)
Return the lth pressure value.
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...
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
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.
Class for Refineable solid mechanics elements in near-incompressible/ incompressible formulation,...
void get_interpolated_values(const unsigned &t, 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 num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
in strain tensor.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
No values are interpolated in this element (pure solid)
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...
virtual Node * solid_pressure_node_pt(const unsigned &l)
void further_build()
Pass the generic stuff down to the sons.
void get_mass_matrix_diagonal(Vector< double > &mass_diag)
Compute the diagonal of the displacement mass matrix for LSC preconditioner.
Class for Refineable PVD equations.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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 Vector< double > &s, Vector< double > &values)
No values are interpolated in this element (pure solid)
virtual Node * solid_pressure_node_pt(const unsigned &l)
void further_build()
Further build function, pass the pointers down to the sons.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
No values are interpolated in this element (pure solid)
void fill_in_generic_contribution_to_residuals_pvd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Call the residuals including hanging node cases.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 0 (pure solid problem)
Refineable version of QElement<3,NNODE_1D>.
Class for refineable solid mechanics elements in near-incompressible/ incompressible formulation,...
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...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void rebuild_from_sons(Mesh *&mesh_pt)
Empty rebuild from sons, empty.
unsigned required_nvalue(const unsigned &n) const
Overload the number of additional solid dofs at each node, we shall always assign 1,...
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 ...
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...
unsigned nvertex_node() const
Number of vertex nodes in the element.
void pin_elemental_redundant_nodal_solid_pressures()
Pin the redundant solid pressure.
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...
unsigned ninterpolating_node_1d(const int &value_id)
The number of 1d pressure nodes is 2, otherwise we have the positional nodes.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
OK, interpolate the solid pressures.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values (1) solid pressure.
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...
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
OK get the time-dependent verion.
Class for refineable solid mechanics elements in near-incompressible/ incompressible formulation,...
void further_setup_hanging_nodes()
No additional hanging node procedures are required for discontinuous solid pressures.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nvertex_node() const
Number of vertex nodes in the element.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 0 (pure solid problem)
void rebuild_from_sons(Mesh *&mesh_pt)
Reconstruct the pressure from the sons Must be specialized for each dimension.
void unpin_elemental_solid_pressure_dofs()
Unpin all solid pressure dofs.
void further_build()
Further build: Interpolate the solid pressure values Again this must be specialised for each dimensio...
Class for refineable QPVDElement elements.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void further_setup_hanging_nodes()
No additional hanging node procedures are required for the solid elements.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void rebuild_from_sons(Mesh *&mesh_pt)
Empty rebuild from sons, no need to reconstruct anything here.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
virtual void further_build()
Further build: Pass the father's Use_undeformed_macro_element_for_new_lagrangian_coords flag down,...
A class that is used to template the solid refineable Q elements by dimension. It's really nothing mo...
Definition: Qelements.h:2286
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
SolidQElement elements are quadrilateral elements whose derivatives also include those based upon the...
Definition: Qelements.h:1742
//////////////////////////////////////////////////////////////////// ////////////////////////////////...