constrained_volume_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 elements that allow the imposition of a "constant volume"
31 // constraint in free surface problems.
32 
33 //Include guards, to prevent multiple includes
34 #ifndef CONSTRAINED_FLUID_VOLUME_ELEMENTS_HEADER
35 #define CONSTRAINED_FLUID_VOLUME_ELEMENTS_HEADER
36 
37 // Config header generated by autoconfig
38 #ifdef HAVE_CONFIG_H
39  #include <oomph-lib-config.h>
40 #endif
41 
42 //OOMPH-LIB headers
43 #include "../generic/Qelements.h"
44 #include "../generic/spines.h"
45 #include "../axisym_navier_stokes/axisym_navier_stokes_elements.h"
46 
47 //-------------------------------------------
48 // NOTE: This is still work
49 // in progress. Need to create versions
50 // for axisymmetric and cartesian 3D
51 // bulk equations, as well as spine
52 // version for all of these (resulting
53 // in six elements in total). These
54 // will gradually replace the very hacky
55 // fix_*.h files that currently live in
56 // various demo driver codes.
57 //-------------------------------------------
58 
59 namespace oomph
60 {
61 
62 //==========================================================================
63 /// A class that is used to implement the constraint that the fluid volume
64 /// in a region bounded by associated FaceElements (attached, e.g., to the
65 /// mesh boundaries that enclose a bubble) must take a specific value.
66 /// This GeneralisedElement is used only to store the desired volume and
67 /// a pointer to the (usually pressure) freedom that must be traded
68 /// for the volume constraint.
69 //=========================================================================
71  {
72  private:
73 
74  /// Pointer to the desired value of the volume
76 
77  /// \short The Data that contains the traded pressure is stored
78  /// as external or internal Data for the element. What is its index
79  /// in these containers?
81 
82  /// \short The Data that contains the traded pressure is stored
83  /// as external or internal Data for the element. Which one?
85 
86  /// Index of the value in traded pressure data that corresponds to the
87  /// traded pressure
89 
90  /// \short The local eqn number for the traded pressure
91  inline int ptraded_local_eqn()
92  {
94  {
95  return
96  this->internal_local_eqn(
99  }
100  else
101  {
102  return
103  this->external_local_eqn(
106  }
107  }
108 
109 
110  /// \short Fill in the residuals for the volume constraint
112  Vector<double> &residuals);
113 
114  public:
115 
116  /// \short Constructor: Pass pointer to target volume. "Pressure" value that
117  /// "traded" for the volume contraint is created internally (as a Data
118  /// item with a single pressure value)
119  VolumeConstraintElement(double* prescribed_volume_pt);
120 
121  /// \short Constructor: Pass pointer to target volume, pointer to Data
122  /// item whose value specified by index_of_traded_pressure represents
123  /// the "Pressure" value that "traded" for the volume contraint.
124  /// The Data is stored as external Data for this element.
125  VolumeConstraintElement(double* prescribed_volume_pt,
127  const unsigned& index_of_traded_pressure);
128 
129  /// \short Empty destructor
131 
132  /// Access to Data that contains the traded pressure
134  {
136  {
137  return internal_data_pt(
139  }
140  else
141  {
142  return external_data_pt(
144  }
145  }
146 
147  /// Return the traded pressure value
148  inline double p_traded()
150 
151  /// Return the index of Data object at which the traded pressure is stored
152  inline unsigned index_of_traded_pressure()
154 
155 
156  /// \short Fill in the residuals for the volume constraint
158  {
160  residuals);
161  }
162 
163  /// \short Fill in the residuals and jacobian for the volume constraint
165  DenseMatrix<double> &jacobian)
166  {
167  //No contribution to jacobian; see comment in that function
169  residuals);
170  }
171 
172  /// \short Fill in the residuals, jacobian and mass matrix for the volume
173  /// constraint.
175  Vector<double> &residuals,
176  DenseMatrix<double> &jacobian,
177  DenseMatrix<double> &mass_matrix)
178  {
179  //No contribution to jacobian or mass matrix; see comment in that function
181  residuals);
182  }
183 
184 };
185 
186 /////////////////////////////////////////////////////////////////////////
187 /////////////////////////////////////////////////////////////////////////
188 /////////////////////////////////////////////////////////////////////////
189 
190 
191 //=======================================================================
192 /// Base class for interface elements that allow the application of
193 /// a volume constraint on the region bounded by these elements. The
194 /// elements must be used together with the associated
195 /// VolumeConstraintElement which stores the value of the
196 /// target volume. Common functionality is provided in this base
197 /// for storing the external "pressure" value
198 /// that is traded for the volume constraint.
199 //=======================================================================
201  {
202 
203  protected:
204 
205  /// \short The Data that contains the traded pressure is usually stored
206  /// as external Data for this element, but may also be nodal Data
207  /// Which Data item is it?
209 
210  /// Index of the value in traded pressure data that corresponds to the
211  /// traded pressure
213 
214  /// \short Boolean to indicate whether the traded pressure is
215  /// stored externally or at a node (this can happen in Taylor-Hood
216  /// elements)
218 
219  /// \short The local eqn number for the traded pressure
220  inline int ptraded_local_eqn()
221  {
222  //Return the appropriate nodal value if required
224  {
227  }
228  else
229  {
232  }
233  }
234 
235  /// \short Helper function to fill in contributions to residuals
236  /// (remember that part of the residual is added by the the
237  /// associated VolumeConstraintElement). This is dimension/geometry
238  /// specific and must be implemented in derived classes for
239  /// 1D line, 2D surface and axisymmetric fluid boundaries
241  Vector<double> &residuals)=0;
242 
243  public:
244 
245  /// \short Constructor initialise the boolean flag
246  /// We expect the traded pressure data to be stored externally
248  {}
249 
250  /// \short Empty Destructor
252 
253  /// Fill in contribution to residuals and Jacobian
255  {
256  //Call the generic routine
258  }
259 
260  /// \short Set the "master" volume constraint element
261  /// The setup here is a bit more complicated than one might expect because
262  /// if an internal pressure on a boundary
263  /// is hijacked in TaylorHood elements then
264  /// that the "traded" value may already be stored as nodal data
265  /// in this element. This causes no problems apart from when running
266  /// with PARANOID in which case the resulting
267  /// repeated local equation numbers are spotted and an error is thrown.
268  /// The check is a finite amount of work and so can be avoided if
269  /// the boolean flag check_nodal_data is set to false.
271  &vol_constraint_el_pt,
272  const bool &check_nodal_data=true)
273  {
274 
275  //In order to buffer the case of nodal data, we (tediously) check that the
276  //traded pressure is not already nodal data of this element
277  if(check_nodal_data)
278  {
279  //Get memory address of the equation indexed by the
280  //traded pressure datum
281  long* global_eqn_number = vol_constraint_el_pt->p_traded_data_pt()
282  ->eqn_number_pt(vol_constraint_el_pt->index_of_traded_pressure());
283 
284  //Put in a check if the datum already corresponds to any other
285  //(nodal) data in the element by checking the memory addresses of
286  //their global equation numbers
287  const unsigned n_node = this->nnode();
288  for(unsigned n=0;n<n_node;n++)
289  {
290  //Cache the node pointer
291  Node* const nod_pt = this->node_pt(n);
292  //Find all nodal data values
293  unsigned n_value = nod_pt->nvalue();
294  //If we already have the data, set the
295  //lookup schemes accordingly and return from the function
296  for(unsigned i=0;i<n_value;i++)
297  {
298  if(nod_pt->eqn_number_pt(i) == global_eqn_number)
299  {
303  //Finished so exit the function
304  return;
305  }
306  }
307  }
308  }
309 
310  //Should only get here if the data is not nodal
311 
312  // Add "traded" pressure data as external data to this element
314  this->add_external_data(vol_constraint_el_pt->p_traded_data_pt());
315 
316  // Which value corresponds to the traded pressure
317  Index_of_traded_pressure_value=vol_constraint_el_pt->
318  index_of_traded_pressure();
319  }
320 
321 };
322 
323 /////////////////////////////////////////////////////////////////////////
324 /////////////////////////////////////////////////////////////////////////
325 /////////////////////////////////////////////////////////////////////////
326 
327 
328 //=======================================================================
329 /// One-dimensional interface elements that allow the application of
330 /// a volume constraint on the region bounded by these elements.
331 /// The volume is computed by integrating x.n around the boundary of
332 /// the domain and then dividing by two.
333 /// The sign is chosen so that the volume will be positive
334 /// when the elements surround a fluid domain.
335 ///
336 /// These elements must be used together with the associated
337 /// VolumeConstraintElement, which stores the value of the
338 /// target volume.
339 //=======================================================================
342  {
343  protected:
344 
345  /// \short Helper function to fill in contributions to residuals
346  /// (remember that part of the residual is added by the
347  /// the associated VolumeConstraintElement). This is specific for
348  /// 1D line elements that bound 2D cartesian fluid elements.
350  Vector<double> &residuals);
351 
352  public:
353 
354  /// \short Empty Contructor
356 
357  /// \short Empty Destructor
359 
360  /// Return this element's contribution to the total volume enclosed
362 
363  };
364 
365 
366 //////////////////////////////////////////////////////////////////////
367 //////////////////////////////////////////////////////////////////////
368 //////////////////////////////////////////////////////////////////////
369 
370 
371 //=======================================================================
372 /// The one-dimensional interface elements that allow imposition of a
373 /// volume constraint specialised for the case when the nodal positions of
374 /// the bulk elements are treated as solid degrees of freedom.
375 /// To enforce that a fluid volume has a
376 /// certain volume, attach these elements to all faces of the
377 /// (2D cartesian) bulk fluid elements (of type ELEMENT) that bound that region
378 /// and then specify the "pressure" value that is traded for the constraint.
379 //=======================================================================
380  template<class ELEMENT>
383  public virtual FaceGeometry<ELEMENT>
384  {
385 
386  public:
387 
388  /// \short Contructor: Specify bulk element and index of face to which
389  /// this face element is to be attached
391  const int &face_index) :
392  FaceGeometry<ELEMENT>(),
394  {
395  //Attach the geometrical information to the element, by
396  //making the face element from the bulk element
397  element_pt->build_face_element(face_index,this);
398  }
399 
400  /// Fill in contribution to residuals and Jacobian. This is specific
401  /// to solid-based elements in which derivatives w.r.t. to nodal
402  /// positions are evaluated by finite differencing
404  DenseMatrix<double> &jacobian)
405  {
406  //Call the generic routine
408 
409  // Shape derivatives
410  //Call the generic finite difference routine for the solid variables
411  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
412  }
413 
414  /// \short The "global" intrinsic coordinate of the element when
415  /// viewed as part of a geometric object should be given by
416  /// the FaceElement representation, by default
417  double zeta_nodal(const unsigned &n, const unsigned &k,
418  const unsigned &i) const
419  {return FaceElement::zeta_nodal(n,k,i);}
420 
421 };
422 
423 //=======================================================================
424 /// The one-dimensional interface elements that allow imposition of a
425 /// volume constraint specialised for the case when the nodal positions of
426 /// the bulk elements are adjusted using Spines.
427 /// To enforce that a fluid volume has a
428 /// certain volume, attach these elements to all faces of the
429 /// (2D cartesian) bulk fluid elements (of type ELEMENT) that bound that region
430 /// and then specify the "pressure" value that is traded for the constraint.
431 //=======================================================================
432  template<class ELEMENT>
435  public virtual SpineElement<FaceGeometry<ELEMENT> >
436  {
437 
438  public:
439 
440  /// \short Contructor: Specify bulk element and index of face to which
441  /// this face element is to be attached.
443  const int &face_index) :
444  SpineElement<FaceGeometry<ELEMENT> >(),
446  {
447  //Attach the geometrical information to the element, by
448  //making the face element from the bulk element
449  element_pt->build_face_element(face_index,this);
450  }
451 
452 
453  /// Fill in contribution to residuals and Jacobian. This is specific
454  /// to spine based elements in which the shape derivatives are evaluated
455  /// using geometric data
457  DenseMatrix<double> &jacobian)
458  {
459  //Call the generic routine
461 
462  //Call the generic routine to evaluate shape derivatives
463  this->fill_in_jacobian_from_geometric_data(jacobian);
464  }
465 
466  /// \short The "global" intrinsic coordinate of the element when
467  /// viewed as part of a geometric object should be given by
468  /// the FaceElement representation, by default
469  double zeta_nodal(const unsigned &n, const unsigned &k,
470  const unsigned &i) const
471  {return FaceElement::zeta_nodal(n,k,i);}
472 
473 };
474 
475 /////////////////////////////////////////////////////////////////////////
476 /////////////////////////////////////////////////////////////////////////
477 /////////////////////////////////////////////////////////////////////////
478 
479 
480 //=======================================================================
481 /// Axisymmetric (one-dimensional) interface elements that
482 /// allow the application of
483 /// a volume constraint on the region bounded by these elements.
484 /// The volume is computed by integrating x.n around the boundary of
485 /// the domain and then dividing by three.
486 /// The sign is chosen so that the volume will be positive
487 /// when the elements surround a fluid domain.
488 ///
489 /// These elements must be used together with the associated
490 /// VolumeConstraintElement, which stores the value of the
491 /// target volume.
492 //=======================================================================
495  {
496  protected:
497 
498  /// \short Helper function to fill in contributions to residuals
499  /// (remember that part of the residual is added by the
500  /// the associated VolumeConstraintElement). This is specific for
501  /// 1D line elements that bound 2D cartesian fluid elements.
503  Vector<double> &residuals);
504 
505  public:
506 
507  /// \short Empty Contructor
510 
511  /// \short Empty Destructor
513 
514  /// Return this element's contribution to the total volume enclosed
516 
517  /// \short Return this element's contribution to the volume flux over
518  /// the boundary
520  {
521  // Initialise
522  double vol=0.0;
523 
524  //Find out how many nodes there are
525  const unsigned n_node = this->nnode();
526 
527  //Set up memeory for the shape functions
528  Shape psif(n_node);
529  DShape dpsifds(n_node,1);
530 
531  //Set the value of n_intpt
532  const unsigned n_intpt = this->integral_pt()->nweight();
533 
534  //Storage for the local cooridinate
535  Vector<double> s(1);
536 
537  //Loop over the integration points
538  for(unsigned ipt=0;ipt<n_intpt;ipt++)
539  {
540  //Get the local coordinate at the integration point
541  s[0] = this->integral_pt()->knot(ipt,0);
542 
543  //Get the integral weight
544  double W = this->integral_pt()->weight(ipt);
545 
546  //Call the derivatives of the shape function at the knot point
547  this->dshape_local_at_knot(ipt,psif,dpsifds);
548 
549  // Get position, tangent vector and velocity vector (r and z
550  // components only)
551  Vector<double> interpolated_u(2,0.0);
552  Vector<double> interpolated_t1(2,0.0);
554  for(unsigned l=0;l<n_node;l++)
555  {
556  //Loop over directional components
557  for(unsigned i=0;i<2;i++)
558  {
559  interpolated_x[i] += this->nodal_position(l,i)*psif(l);
560  interpolated_u[i] += this->node_pt(l)->value(
561  dynamic_cast<AxisymmetricNavierStokesEquations*>(bulk_element_pt())
562  ->u_index_axi_nst(i))*psif(l);
563  interpolated_t1[i] += this->nodal_position(l,i)*dpsifds(l,0);
564  }
565  }
566 
567  //Calculate the length of the tangent Vector
568  double tlength = interpolated_t1[0]*interpolated_t1[0] +
569  interpolated_t1[1]*interpolated_t1[1];
570 
571  //Set the Jacobian of the line element
572  double J = sqrt(tlength)*interpolated_x[0];
573 
574  //Now calculate the normal Vector
575  Vector<double> interpolated_n(2);
576  this->outer_unit_normal(ipt,interpolated_n);
577 
578  // Assemble dot product
579  double dot = 0.0;
580  for(unsigned k=0;k<2;k++)
581  {
582  dot += interpolated_u[k]*interpolated_n[k];
583  }
584 
585  // Add to volume with sign chosen so that the volume is
586  // positive when the elements bound the fluid
587  vol += dot*W*J;
588  }
589 
590  return vol;
591  }
592 
593  };
594 
595 //////////////////////////////////////////////////////////////////////
596 //////////////////////////////////////////////////////////////////////
597 //////////////////////////////////////////////////////////////////////
598 
599 
600 //=======================================================================
601 /// The axisymmetric (one-dimensional) interface elements that allow
602 /// imposition of a
603 /// volume constraint specialised for the case when the nodal positions of
604 /// the bulk elements are treated as solid degrees of freedom.
605 /// To enforce that a fluid volume has a
606 /// certain volume, attach these elements to all faces of the
607 /// (2D axisymmetric) bulk fluid elements (of type ELEMENT)
608 /// that bound that region
609 /// and then specify the "pressure" value that is traded for the constraint.
610 //=======================================================================
611  template<class ELEMENT>
614  public virtual FaceGeometry<ELEMENT>
615  {
616 
617  public:
618 
619  /// \short Contructor: Specify bulk element and index of face to which
620  /// this face element is to be attached
622  FiniteElement* const &element_pt,
623  const int &face_index) :
624  FaceGeometry<ELEMENT>(),
626  {
627  //Attach the geometrical information to the element, by
628  //making the face element from the bulk element
629  element_pt->build_face_element(face_index,this);
630  }
631 
632  /// Fill in contribution to residuals and Jacobian. This is specific
633  /// to solid-based elements in which derivatives w.r.t. to nodal
634  /// positions are evaluated by finite differencing
636  DenseMatrix<double> &jacobian)
637  {
638  //Call the generic routine
640 
641  // Shape derivatives
642  //Call the generic finite difference routine for the solid variables
643  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
644  }
645 
646  /// \short The "global" intrinsic coordinate of the element when
647  /// viewed as part of a geometric object should be given by
648  /// the FaceElement representation, by default
649  double zeta_nodal(const unsigned &n, const unsigned &k,
650  const unsigned &i) const
651  {return FaceElement::zeta_nodal(n,k,i);}
652 
653 };
654 
655 //=======================================================================
656 /// The axisymmetric (one-dimensional) interface elements that
657 /// allow imposition of a
658 /// volume constraint specialised for the case when the nodal positions of
659 /// the bulk elements are adjusted using Spines.
660 /// To enforce that a fluid volume has a
661 /// certain volume, attach these elements to all faces of the
662 /// (2D axisymmetric) bulk fluid elements (of type ELEMENT) that bound
663 /// that region
664 /// and then specify the "pressure" value that is traded for the constraint.
665 //=======================================================================
666  template<class ELEMENT>
669  public virtual SpineElement<FaceGeometry<ELEMENT> >
670  {
671 
672  public:
673 
674  /// \short Contructor: Specify bulk element and index of face to which
675  /// this face element is to be attached.
677  FiniteElement* const &element_pt, const int &face_index) :
678  SpineElement<FaceGeometry<ELEMENT> >(),
680  {
681  //Attach the geometrical information to the element, by
682  //making the face element from the bulk element
683  element_pt->build_face_element(face_index,this);
684  }
685 
686 
687  /// Fill in contribution to residuals and Jacobian. This is specific
688  /// to spine based elements in which the shape derivatives are evaluated
689  /// using geometric data
691  DenseMatrix<double> &jacobian)
692  {
693  //Call the generic routine
695 
696  //Call the generic routine to evaluate shape derivatives
697  this->fill_in_jacobian_from_geometric_data(jacobian);
698  }
699 
700  /// \short The "global" intrinsic coordinate of the element when
701  /// viewed as part of a geometric object should be given by
702  /// the FaceElement representation, by default
703  double zeta_nodal(const unsigned &n, const unsigned &k,
704  const unsigned &i) const
705  {return FaceElement::zeta_nodal(n,k,i);}
706 
707 };
708 
709 
710 
711 /////////////////////////////////////////////////////////////////////////
712 /////////////////////////////////////////////////////////////////////////
713 /////////////////////////////////////////////////////////////////////////
714 
715 
716 //=======================================================================
717 /// Two-dimensional interface elements that allow the application of
718 /// a volume constraint on the region bounded by these elements.
719 /// The volume is computed by integrating x.n around the boundary of
720 /// the domain and then dividing by three.
721 /// The sign is chosen so that the volume will be positive
722 /// when the elements surround a fluid domain.
723 ///
724 /// These elements must be used together with the associated
725 /// VolumeConstraintElement, which stores the value of the
726 /// target volume.
727 //=======================================================================
730  {
731 
732  protected:
733 
734  /// \short Helper function to fill in contributions to residuals
735  /// (remember that part of the residual is added by the
736  /// the associated VolumeConstraintElement). This is specific for
737  /// 2D surface elements that bound 3D cartesian fluid elements.
739  Vector<double> &residuals);
740 
741  public:
742 
743  /// \short Empty Contructor
745  {}
746 
747  /// \short Empty Desctructor
749 
750 };
751 
752 
753 //////////////////////////////////////////////////////////////////////
754 //////////////////////////////////////////////////////////////////////
755 //////////////////////////////////////////////////////////////////////
756 
757 
758 //=======================================================================
759 /// The Two-dimensional interface elements that allow the application of
760 /// a volume constraint specialised for the case when the nodal positions
761 /// of the bulk elements are treated as solid degrees of freedom.
762 /// To enforce that a fluid volume has a
763 /// certain volume, attach these elements to all faces of the
764 /// (3D Cartesian) bulk fluid elements (of type ELEMENT) that bound that region
765 /// and then specify the "pressure" value that is traded for the constraint.
766 //=======================================================================
767  template<class ELEMENT>
770  public virtual FaceGeometry<ELEMENT>
771  {
772 
773  public:
774 
775  /// \short Contructor: Specify bulk element and index of face to which
776  /// this face element is to be attached.
778  FiniteElement* const &element_pt, const int &face_index) :
779  FaceGeometry<ELEMENT>(),
781  {
782  //Attach the geometrical information to the element, by
783  //making the face element from the bulk element
784  element_pt->build_face_element(face_index,this);
785  }
786 
787 
788 
789  /// Fill in contribution to residuals and Jacobian. This is specific
790  /// to solid-based elements in which derivatives w.r.t. to nodal
791  /// positions are evaluated by finite differencing
793  DenseMatrix<double> &jacobian)
794  {
795  //Call the generic routine
797 
798  // Shape derivatives
799  //Call the generic finite difference routine for the solid variables
800  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
801  }
802 
803 
804  /// \short The "global" intrinsic coordinate of the element when
805  /// viewed as part of a geometric object should be given by
806  /// the FaceElement representation, by default
807  double zeta_nodal(const unsigned &n, const unsigned &k,
808  const unsigned &i) const
809  {return FaceElement::zeta_nodal(n,k,i);}
810 
811  };
812 
813 
814 /////////////////////////////////////////////////////////////////////////
815 /////////////////////////////////////////////////////////////////////////
816 /////////////////////////////////////////////////////////////////////////
817 
818 
819 //=======================================================================
820 /// The Two-dimensional interface elements that allow the application of
821 /// a volume constraint specialised for the case when the nodal positions
822 /// of the bulk elements are adjusted using spines.
823 /// To enforce that a fluid volume has a
824 /// certain volume, attach these elements to all faces of the
825 /// (3D Cartesian) bulk fluid elements (of type ELEMENT) that bound that region
826 /// and then specify the "pressure" value that is traded for the constraint.
827 //=======================================================================
828  template<class ELEMENT>
831  public virtual SpineElement<FaceGeometry<ELEMENT> >
832 {
833 
834  public:
835 
836  /// \short Contructor: Specify bulk element and index of face to which
837  /// this face element is to be attached.
839  const int &face_index) :
840  SpineElement<FaceGeometry<ELEMENT> >(),
842  {
843  //Attach the geometrical information to the element, by
844  //making the face element from the bulk element
845  element_pt->build_face_element(face_index,this);
846  }
847 
848 
849 
850  /// Fill in contribution to residuals and Jacobian. This is specific
851  /// to spine based elements in which the shape derivatives are evaluated
852  /// using geometric data
854  DenseMatrix<double> &jacobian)
855  {
856  //Call the generic routine
858 
859  //Call the generic routine to evaluate shape derivatives
860  this->fill_in_jacobian_from_geometric_data(jacobian);
861  }
862 
863 
864  /// \short The "global" intrinsic coordinate of the element when
865  /// viewed as part of a geometric object should be given by
866  /// the FaceElement representation, by default
867  double zeta_nodal(const unsigned &n, const unsigned &k,
868  const unsigned &i) const
869  {return FaceElement::zeta_nodal(n,k,i);}
870 
871 };
872 
873 
874 }
875 #endif
876 
877 
878 
879 
880 
881 
VolumeConstraintBoundingElement()
Constructor initialise the boolean flag We expect the traded pressure data to be stored externally...
A Generalised Element class.
Definition: elements.h:76
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
ElasticAxisymmetricVolumeConstraintBoundingElement(FiniteElement *const &element_pt, const int &face_index)
Contructor: Specify bulk element and index of face to which this face element is to be attached...
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
double contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
long * eqn_number_pt(const unsigned &i)
Return the pointer to the equation number of the i-th stored variable.
Definition: nodes.h:356
SpineLineVolumeConstraintBoundingElement(FiniteElement *const &element_pt, const int &face_index)
Contructor: Specify bulk element and index of face to which this face element is to be attached...
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void fill_in_generic_residual_contribution_volume_constraint(Vector< double > &residuals)
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
void fill_in_generic_contribution_to_residuals_volume_constraint(Vector< double > &residuals)
Fill in the residuals for the volume constraint.
A general Finite Element class.
Definition: elements.h:1271
ElasticSurfaceVolumeConstraintBoundingElement(FiniteElement *const &element_pt, const int &face_index)
Contructor: Specify bulk element and index of face to which this face element is to be attached...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:5695
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the residuals and jacobian for the volume constraint.
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
bool Traded_pressure_stored_at_node
Boolean to indicate whether the traded pressure is stored externally or at a node (this can happen in...
bool Traded_pressure_stored_as_internal_data
The Data that contains the traded pressure is stored as external or internal Data for the element...
void fill_in_generic_residual_contribution_volume_constraint(Vector< double > &residuals)
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
SpineAxisymmetricVolumeConstraintBoundingElement(FiniteElement *const &element_pt, const int &face_index)
Contructor: Specify bulk element and index of face to which this face element is to be attached...
VolumeConstraintElement(double *prescribed_volume_pt)
Constructor: Pass pointer to target volume. "Pressure" value that "traded" for the volume contraint i...
double contribution_to_volume_flux()
Return this element's contribution to the volume flux over the boundary.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Fill in the residuals, jacobian and mass matrix for the volume constraint.
static char t char * s
Definition: cfortran.h:572
virtual void fill_in_generic_residual_contribution_volume_constraint(Vector< double > &residuals)=0
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
void set_volume_constraint_element(VolumeConstraintElement *const &vol_constraint_el_pt, const bool &check_nodal_data=true)
Set the "master" volume constraint element The setup here is a bit more complicated than one might ex...
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
double p_traded()
Return the traded pressure value.
double contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
SpineSurfaceVolumeConstraintBoundingElement(FiniteElement *const &element_pt, const int &face_index)
Contructor: Specify bulk element and index of face to which this face element is to be attached...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
int ptraded_local_eqn()
The local eqn number for the traded pressure.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary, when viewed as part of a compound geometric object is specified using the boundary coordinate defined by the mesh. Note: Boundary coordinates will have been set up when creating the underlying mesh, and their values will have been stored at the nodes.
Definition: elements.h:4251
double * Prescribed_volume_pt
Pointer to the desired value of the volume.
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
Definition: elements.cc:3174
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
unsigned Data_number_of_traded_pressure
The Data that contains the traded pressure is usually stored as external Data for this element...
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
Definition: elements.h:662
void fill_in_generic_residual_contribution_volume_constraint(Vector< double > &residuals)
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
double dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
Definition: Vector.h:289
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4470
The SpineElement<ELEMENT> class takes an existing element as a template parameter and adds the necess...
Definition: spines.h:423
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
int ptraded_local_eqn()
The local eqn number for the traded pressure.
unsigned External_or_internal_data_index_of_traded_pressure
The Data that contains the traded pressure is stored as external or internal Data for the element...
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
Data * p_traded_data_pt()
Access to Data that contains the traded pressure.
void fill_in_jacobian_from_geometric_data(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the Jacobian matrix from the geometric data. This version assumes that...
unsigned index_of_traded_pressure()
Return the index of Data object at which the traded pressure is stored.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in contribution to residuals and Jacobian.
int external_local_eqn(const unsigned &i, const unsigned &j)
Return the local equation number corresponding to the j-th value stored at the i-th external data...
Definition: elements.h:313
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals for the volume constraint.
ElasticLineVolumeConstraintBoundingElement(FiniteElement *const &element_pt, const int &face_index)
Contructor: Specify bulk element and index of face to which this face element is to be attached...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...