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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for elements that allow the imposition of a "constant volume"
27 // constraint in free surface problems.
28 
29 // Include guards, to prevent multiple includes
30 #ifndef CONSTRAINED_FLUID_VOLUME_ELEMENTS_HEADER
31 #define CONSTRAINED_FLUID_VOLUME_ELEMENTS_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35 #include <oomph-lib-config.h>
36 #endif
37 
38 // OOMPH-LIB headers
39 #include "../generic/Qelements.h"
40 #include "../generic/spines.h"
41 #include "../axisym_navier_stokes/axisym_navier_stokes_elements.h"
42 
43 //-------------------------------------------
44 // NOTE: This is still work
45 // in progress. Need to create versions
46 // for axisymmetric and cartesian 3D
47 // bulk equations, as well as spine
48 // version for all of these (resulting
49 // in six elements in total). These
50 // will gradually replace the very hacky
51 // fix_*.h files that currently live in
52 // various demo driver codes.
53 //-------------------------------------------
54 
55 namespace oomph
56 {
57  //==========================================================================
58  /// A class that is used to implement the constraint that the fluid volume
59  /// in a region bounded by associated FaceElements (attached, e.g., to the
60  /// mesh boundaries that enclose a bubble) must take a specific value.
61  /// This GeneralisedElement is used only to store the desired volume and
62  /// a pointer to the (usually pressure) freedom that must be traded
63  /// for the volume constraint.
64  //=========================================================================
65  class VolumeConstraintElement : public GeneralisedElement
66  {
67  private:
68  /// Pointer to the desired value of the volume
70 
71  /// The Data that contains the traded pressure is stored
72  /// as external or internal Data for the element. What is its index
73  /// in these containers?
75 
76  /// The Data that contains the traded pressure is stored
77  /// as external or internal Data for the element. Which one?
79 
80  /// Index of the value in traded pressure data that corresponds to the
81  /// traded pressure
83 
84  /// The local eqn number for the traded pressure
85  inline int ptraded_local_eqn()
86  {
88  {
89  return this->internal_local_eqn(
92  }
93  else
94  {
95  return this->external_local_eqn(
98  }
99  }
100 
101 
102  /// Fill in the residuals for the volume constraint
104  Vector<double>& residuals);
105 
106  public:
107  /// Constructor: Pass pointer to target volume. "Pressure" value that
108  /// "traded" for the volume contraint is created internally (as a Data
109  /// item with a single pressure value)
110  VolumeConstraintElement(double* prescribed_volume_pt);
111 
112  /// Constructor: Pass pointer to target volume, pointer to Data
113  /// item whose value specified by index_of_traded_pressure represents
114  /// the "Pressure" value that "traded" for the volume contraint.
115  /// The Data is stored as external Data for this element.
116  VolumeConstraintElement(double* prescribed_volume_pt,
117  Data* p_traded_data_pt,
118  const unsigned& index_of_traded_pressure);
119 
120  /// Empty destructor
122 
123  /// Access to Data that contains the traded pressure
124  inline Data* p_traded_data_pt()
125  {
127  {
128  return internal_data_pt(
130  }
131  else
132  {
133  return external_data_pt(
135  }
136  }
137 
138  /// Return the traded pressure value
139  inline double p_traded()
140  {
142  }
143 
144  /// Return the index of Data object at which the traded pressure is stored
145  inline unsigned index_of_traded_pressure()
146  {
148  }
149 
150 
151  /// Fill in the residuals for the volume constraint
152  void fill_in_contribution_to_residuals(Vector<double>& residuals)
153  {
155  residuals);
156  }
157 
158  /// Fill in the residuals and jacobian for the volume constraint
159  void fill_in_contribution_to_jacobian(Vector<double>& residuals,
160  DenseMatrix<double>& jacobian)
161  {
162  // No contribution to jacobian; see comment in that function
164  residuals);
165  }
166 
167  /// Fill in the residuals, jacobian and mass matrix for the volume
168  /// constraint.
170  Vector<double>& residuals,
171  DenseMatrix<double>& jacobian,
172  DenseMatrix<double>& mass_matrix)
173  {
174  // No contribution to jacobian or mass matrix; see comment in that
175  // function
177  residuals);
178  }
179  };
180 
181  /// //////////////////////////////////////////////////////////////////////
182  /// //////////////////////////////////////////////////////////////////////
183  /// //////////////////////////////////////////////////////////////////////
184 
185 
186  //=======================================================================
187  /// Base class for interface elements that allow the application of
188  /// a volume constraint on the region bounded by these elements. The
189  /// elements must be used together with the associated
190  /// VolumeConstraintElement which stores the value of the
191  /// target volume. Common functionality is provided in this base
192  /// for storing the external "pressure" value
193  /// that is traded for the volume constraint.
194  //=======================================================================
195  class VolumeConstraintBoundingElement : public virtual FaceElement
196  {
197  protected:
198  /// The Data that contains the traded pressure is usually stored
199  /// as external Data for this element, but may also be nodal Data
200  /// Which Data item is it?
202 
203  /// Index of the value in traded pressure data that corresponds to the
204  /// traded pressure
206 
207  /// Boolean to indicate whether the traded pressure is
208  /// stored externally or at a node (this can happen in Taylor-Hood
209  /// elements)
211 
212  /// The local eqn number for the traded pressure
213  inline int ptraded_local_eqn()
214  {
215  // Return the appropriate nodal value if required
217  {
218  return this->nodal_local_eqn(Data_number_of_traded_pressure,
220  }
221  else
222  {
223  return this->external_local_eqn(Data_number_of_traded_pressure,
225  }
226  }
227 
228  /// Helper function to fill in contributions to residuals
229  /// (remember that part of the residual is added by the the
230  /// associated VolumeConstraintElement). This is dimension/geometry
231  /// specific and must be implemented in derived classes for
232  /// 1D line, 2D surface and axisymmetric fluid boundaries
234  Vector<double>& residuals) = 0;
235 
236  public:
237  /// Constructor initialise the boolean flag
238  /// We expect the traded pressure data to be stored externally
240 
241  /// Empty Destructor
243 
244  /// Fill in contribution to residuals and Jacobian
245  void fill_in_contribution_to_residuals(Vector<double>& residuals)
246  {
247  // Call the generic routine
249  }
250 
251  /// Set the "master" volume constraint element
252  /// The setup here is a bit more complicated than one might expect because
253  /// if an internal pressure on a boundary
254  /// is hijacked in TaylorHood elements then
255  /// that the "traded" value may already be stored as nodal data
256  /// in this element. This causes no problems apart from when running
257  /// with PARANOID in which case the resulting
258  /// repeated local equation numbers are spotted and an error is thrown.
259  /// The check is a finite amount of work and so can be avoided if
260  /// the boolean flag check_nodal_data is set to false.
262  VolumeConstraintElement* const& vol_constraint_el_pt,
263  const bool& check_nodal_data = true)
264  {
265  // In order to buffer the case of nodal data, we (tediously) check that
266  // the traded pressure is not already nodal data of this element
267  if (check_nodal_data)
268  {
269  // Get memory address of the equation indexed by the
270  // traded pressure datum
271  long* global_eqn_number =
272  vol_constraint_el_pt->p_traded_data_pt()->eqn_number_pt(
273  vol_constraint_el_pt->index_of_traded_pressure());
274 
275  // Put in a check if the datum already corresponds to any other
276  //(nodal) data in the element by checking the memory addresses of
277  // their global equation numbers
278  const unsigned n_node = this->nnode();
279  for (unsigned n = 0; n < n_node; n++)
280  {
281  // Cache the node pointer
282  Node* const nod_pt = this->node_pt(n);
283  // Find all nodal data values
284  unsigned n_value = nod_pt->nvalue();
285  // If we already have the data, set the
286  // lookup schemes accordingly and return from the function
287  for (unsigned i = 0; i < n_value; i++)
288  {
289  if (nod_pt->eqn_number_pt(i) == global_eqn_number)
290  {
294  // Finished so exit the function
295  return;
296  }
297  }
298  }
299  }
300 
301  // Should only get here if the data is not nodal
302 
303  // Add "traded" pressure data as external data to this element
305  this->add_external_data(vol_constraint_el_pt->p_traded_data_pt());
306 
307  // Which value corresponds to the traded pressure
309  vol_constraint_el_pt->index_of_traded_pressure();
310  }
311  };
312 
313  /// //////////////////////////////////////////////////////////////////////
314  /// //////////////////////////////////////////////////////////////////////
315  /// //////////////////////////////////////////////////////////////////////
316 
317 
318  //=======================================================================
319  /// One-dimensional interface elements that allow the application of
320  /// a volume constraint on the region bounded by these elements.
321  /// The volume is computed by integrating x.n around the boundary of
322  /// the domain and then dividing by two.
323  /// The sign is chosen so that the volume will be positive
324  /// when the elements surround a fluid domain.
325  ///
326  /// These elements must be used together with the associated
327  /// VolumeConstraintElement, which stores the value of the
328  /// target volume.
329  //=======================================================================
332  {
333  protected:
334  /// Helper function to fill in contributions to residuals
335  /// (remember that part of the residual is added by the
336  /// the associated VolumeConstraintElement). This is specific for
337  /// 1D line elements that bound 2D cartesian fluid elements.
339  Vector<double>& residuals);
340 
341  public:
342  /// Empty Contructor
344 
345  /// Empty Destructor
347 
348  /// Return this element's contribution to the total volume enclosed
350  };
351 
352 
353  /// ///////////////////////////////////////////////////////////////////
354  /// ///////////////////////////////////////////////////////////////////
355  /// ///////////////////////////////////////////////////////////////////
356 
357 
358  //=======================================================================
359  /// The one-dimensional interface elements that allow imposition of a
360  /// volume constraint specialised for the case when the nodal positions of
361  /// the bulk elements are treated as solid degrees of freedom.
362  /// To enforce that a fluid volume has a
363  /// certain volume, attach these elements to all faces of the
364  /// (2D cartesian) bulk fluid elements (of type ELEMENT) that bound that
365  /// region and then specify the "pressure" value that is traded for the
366  /// constraint.
367  //=======================================================================
368  template<class ELEMENT>
371  public virtual FaceGeometry<ELEMENT>
372  {
373  public:
374  /// Contructor: Specify bulk element and index of face to which
375  /// this face element is to be attached
376  ElasticLineVolumeConstraintBoundingElement(FiniteElement* const& element_pt,
377  const int& face_index)
378  : FaceGeometry<ELEMENT>(), LineVolumeConstraintBoundingElement()
379  {
380  // Attach the geometrical information to the element, by
381  // making the face element from the bulk element
382  element_pt->build_face_element(face_index, this);
383  }
384 
385  /// Fill in contribution to residuals and Jacobian. This is specific
386  /// to solid-based elements in which derivatives w.r.t. to nodal
387  /// positions are evaluated by finite differencing
388  void fill_in_contribution_to_jacobian(Vector<double>& residuals,
389  DenseMatrix<double>& jacobian)
390  {
391  // Call the generic routine
393 
394  // Shape derivatives
395  // Call the generic finite difference routine for the solid variables
396  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
397  }
398 
399  /// The "global" intrinsic coordinate of the element when
400  /// viewed as part of a geometric object should be given by
401  /// the FaceElement representation, by default
402  double zeta_nodal(const unsigned& n,
403  const unsigned& k,
404  const unsigned& i) const
405  {
406  return FaceElement::zeta_nodal(n, k, i);
407  }
408  };
409 
410  //=======================================================================
411  /// The one-dimensional interface elements that allow imposition of a
412  /// volume constraint specialised for the case when the nodal positions of
413  /// the bulk elements are adjusted using Spines.
414  /// To enforce that a fluid volume has a
415  /// certain volume, attach these elements to all faces of the
416  /// (2D cartesian) bulk fluid elements (of type ELEMENT) that bound that
417  /// region and then specify the "pressure" value that is traded for the
418  /// constraint.
419  //=======================================================================
420  template<class ELEMENT>
423  public virtual SpineElement<FaceGeometry<ELEMENT>>
424  {
425  public:
426  /// Contructor: Specify bulk element and index of face to which
427  /// this face element is to be attached.
428  SpineLineVolumeConstraintBoundingElement(FiniteElement* const& element_pt,
429  const int& face_index)
430  : SpineElement<FaceGeometry<ELEMENT>>(),
432  {
433  // Attach the geometrical information to the element, by
434  // making the face element from the bulk element
435  element_pt->build_face_element(face_index, this);
436  }
437 
438 
439  /// Fill in contribution to residuals and Jacobian. This is specific
440  /// to spine based elements in which the shape derivatives are evaluated
441  /// using geometric data
442  void fill_in_contribution_to_jacobian(Vector<double>& residuals,
443  DenseMatrix<double>& jacobian)
444  {
445  // Call the generic routine
447 
448  // Call the generic routine to evaluate shape derivatives
449  this->fill_in_jacobian_from_geometric_data(jacobian);
450  }
451 
452  /// The "global" intrinsic coordinate of the element when
453  /// viewed as part of a geometric object should be given by
454  /// the FaceElement representation, by default
455  double zeta_nodal(const unsigned& n,
456  const unsigned& k,
457  const unsigned& i) const
458  {
459  return FaceElement::zeta_nodal(n, k, i);
460  }
461  };
462 
463  /// //////////////////////////////////////////////////////////////////////
464  /// //////////////////////////////////////////////////////////////////////
465  /// //////////////////////////////////////////////////////////////////////
466 
467 
468  //=======================================================================
469  /// Axisymmetric (one-dimensional) interface elements that
470  /// allow the application of
471  /// a volume constraint on the region bounded by these elements.
472  /// The volume is computed by integrating x.n around the boundary of
473  /// the domain and then dividing by three.
474  /// The sign is chosen so that the volume will be positive
475  /// when the elements surround a fluid domain.
476  ///
477  /// These elements must be used together with the associated
478  /// VolumeConstraintElement, which stores the value of the
479  /// target volume.
480  //=======================================================================
483  {
484  protected:
485  /// Helper function to fill in contributions to residuals
486  /// (remember that part of the residual is added by the
487  /// the associated VolumeConstraintElement). This is specific for
488  /// 1D line elements that bound 2D cartesian fluid elements.
490  Vector<double>& residuals);
491 
492  public:
493  /// Empty Contructor
496  {
497  }
498 
499  /// Empty Destructor
501 
502  /// Return this element's contribution to the total volume enclosed
504 
505  /// Return this element's contribution to the volume flux over
506  /// the boundary
508  {
509  // Initialise
510  double vol = 0.0;
511 
512  // Find out how many nodes there are
513  const unsigned n_node = this->nnode();
514 
515  // Set up memeory for the shape functions
516  Shape psif(n_node);
517  DShape dpsifds(n_node, 1);
518 
519  // Set the value of n_intpt
520  const unsigned n_intpt = this->integral_pt()->nweight();
521 
522  // Storage for the local cooridinate
523  Vector<double> s(1);
524 
525  // Loop over the integration points
526  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
527  {
528  // Get the local coordinate at the integration point
529  s[0] = this->integral_pt()->knot(ipt, 0);
530 
531  // Get the integral weight
532  double W = this->integral_pt()->weight(ipt);
533 
534  // Call the derivatives of the shape function at the knot point
535  this->dshape_local_at_knot(ipt, psif, dpsifds);
536 
537  // Get position, tangent vector and velocity vector (r and z
538  // components only)
539  Vector<double> interpolated_u(2, 0.0);
540  Vector<double> interpolated_t1(2, 0.0);
541  Vector<double> interpolated_x(2, 0.0);
542  for (unsigned l = 0; l < n_node; l++)
543  {
544  // Loop over directional components
545  for (unsigned i = 0; i < 2; i++)
546  {
547  interpolated_x[i] += this->nodal_position(l, i) * psif(l);
548  interpolated_u[i] +=
549  this->node_pt(l)->value(
550  dynamic_cast<AxisymmetricNavierStokesEquations*>(
551  bulk_element_pt())
552  ->u_index_axi_nst(i)) *
553  psif(l);
554  interpolated_t1[i] += this->nodal_position(l, i) * dpsifds(l, 0);
555  }
556  }
557 
558  // Calculate the length of the tangent Vector
559  double tlength = interpolated_t1[0] * interpolated_t1[0] +
560  interpolated_t1[1] * interpolated_t1[1];
561 
562  // Set the Jacobian of the line element
563  double J = sqrt(tlength) * interpolated_x[0];
564 
565  // Now calculate the normal Vector
566  Vector<double> interpolated_n(2);
567  this->outer_unit_normal(ipt, interpolated_n);
568 
569  // Assemble dot product
570  double dot = 0.0;
571  for (unsigned k = 0; k < 2; k++)
572  {
573  dot += interpolated_u[k] * interpolated_n[k];
574  }
575 
576  // Add to volume with sign chosen so that the volume is
577  // positive when the elements bound the fluid
578  vol += dot * W * J;
579  }
580 
581  return vol;
582  }
583  };
584 
585  /// ///////////////////////////////////////////////////////////////////
586  /// ///////////////////////////////////////////////////////////////////
587  /// ///////////////////////////////////////////////////////////////////
588 
589 
590  //=======================================================================
591  /// The axisymmetric (one-dimensional) interface elements that allow
592  /// imposition of a
593  /// volume constraint specialised for the case when the nodal positions of
594  /// the bulk elements are treated as solid degrees of freedom.
595  /// To enforce that a fluid volume has a
596  /// certain volume, attach these elements to all faces of the
597  /// (2D axisymmetric) bulk fluid elements (of type ELEMENT)
598  /// that bound that region
599  /// and then specify the "pressure" value that is traded for the constraint.
600  //=======================================================================
601  template<class ELEMENT>
604  public virtual FaceGeometry<ELEMENT>
605  {
606  public:
607  /// Contructor: Specify bulk element and index of face to which
608  /// this face element is to be attached
610  FiniteElement* const& element_pt, const int& face_index)
611  : FaceGeometry<ELEMENT>(), AxisymmetricVolumeConstraintBoundingElement()
612  {
613  // Attach the geometrical information to the element, by
614  // making the face element from the bulk element
615  element_pt->build_face_element(face_index, this);
616  }
617 
618  /// Fill in contribution to residuals and Jacobian. This is specific
619  /// to solid-based elements in which derivatives w.r.t. to nodal
620  /// positions are evaluated by finite differencing
621  void fill_in_contribution_to_jacobian(Vector<double>& residuals,
622  DenseMatrix<double>& jacobian)
623  {
624  // Call the generic routine
626 
627  // Shape derivatives
628  // Call the generic finite difference routine for the solid variables
629  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
630  }
631 
632  /// The "global" intrinsic coordinate of the element when
633  /// viewed as part of a geometric object should be given by
634  /// the FaceElement representation, by default
635  double zeta_nodal(const unsigned& n,
636  const unsigned& k,
637  const unsigned& i) const
638  {
639  return FaceElement::zeta_nodal(n, k, i);
640  }
641  };
642 
643  //=======================================================================
644  /// The axisymmetric (one-dimensional) interface elements that
645  /// allow imposition of a
646  /// volume constraint specialised for the case when the nodal positions of
647  /// the bulk elements are adjusted using Spines.
648  /// To enforce that a fluid volume has a
649  /// certain volume, attach these elements to all faces of the
650  /// (2D axisymmetric) bulk fluid elements (of type ELEMENT) that bound
651  /// that region
652  /// and then specify the "pressure" value that is traded for the constraint.
653  //=======================================================================
654  template<class ELEMENT>
657  public virtual SpineElement<FaceGeometry<ELEMENT>>
658  {
659  public:
660  /// Contructor: Specify bulk element and index of face to which
661  /// this face element is to be attached.
663  FiniteElement* const& element_pt, const int& face_index)
664  : SpineElement<FaceGeometry<ELEMENT>>(),
666  {
667  // Attach the geometrical information to the element, by
668  // making the face element from the bulk element
669  element_pt->build_face_element(face_index, this);
670  }
671 
672 
673  /// Fill in contribution to residuals and Jacobian. This is specific
674  /// to spine based elements in which the shape derivatives are evaluated
675  /// using geometric data
676  void fill_in_contribution_to_jacobian(Vector<double>& residuals,
677  DenseMatrix<double>& jacobian)
678  {
679  // Call the generic routine
681 
682  // Call the generic routine to evaluate shape derivatives
683  this->fill_in_jacobian_from_geometric_data(jacobian);
684  }
685 
686  /// The "global" intrinsic coordinate of the element when
687  /// viewed as part of a geometric object should be given by
688  /// the FaceElement representation, by default
689  double zeta_nodal(const unsigned& n,
690  const unsigned& k,
691  const unsigned& i) const
692  {
693  return FaceElement::zeta_nodal(n, k, i);
694  }
695  };
696 
697 
698  /// //////////////////////////////////////////////////////////////////////
699  /// //////////////////////////////////////////////////////////////////////
700  /// //////////////////////////////////////////////////////////////////////
701 
702 
703  //=======================================================================
704  /// Two-dimensional interface elements that allow the application of
705  /// a volume constraint on the region bounded by these elements.
706  /// The volume is computed by integrating x.n around the boundary of
707  /// the domain and then dividing by three.
708  /// The sign is chosen so that the volume will be positive
709  /// when the elements surround a fluid domain.
710  ///
711  /// These elements must be used together with the associated
712  /// VolumeConstraintElement, which stores the value of the
713  /// target volume.
714  //=======================================================================
717  {
718  protected:
719  /// Helper function to fill in contributions to residuals
720  /// (remember that part of the residual is added by the
721  /// the associated VolumeConstraintElement). This is specific for
722  /// 2D surface elements that bound 3D cartesian fluid elements.
724  Vector<double>& residuals);
725 
726  public:
727  /// Empty Contructor
729  {
730  }
731 
732  /// Empty Desctructor
734  };
735 
736 
737  /// ///////////////////////////////////////////////////////////////////
738  /// ///////////////////////////////////////////////////////////////////
739  /// ///////////////////////////////////////////////////////////////////
740 
741 
742  //=======================================================================
743  /// The Two-dimensional interface elements that allow the application of
744  /// a volume constraint specialised for the case when the nodal positions
745  /// of the bulk elements are treated as solid degrees of freedom.
746  /// To enforce that a fluid volume has a
747  /// certain volume, attach these elements to all faces of the
748  /// (3D Cartesian) bulk fluid elements (of type ELEMENT) that bound that
749  /// region and then specify the "pressure" value that is traded for the
750  /// constraint.
751  //=======================================================================
752  template<class ELEMENT>
755  public virtual FaceGeometry<ELEMENT>
756  {
757  public:
758  /// Contructor: Specify bulk element and index of face to which
759  /// this face element is to be attached.
761  FiniteElement* const& element_pt, const int& face_index)
762  : FaceGeometry<ELEMENT>(), SurfaceVolumeConstraintBoundingElement()
763  {
764  // Attach the geometrical information to the element, by
765  // making the face element from the bulk element
766  element_pt->build_face_element(face_index, this);
767  }
768 
769 
770  /// Fill in contribution to residuals and Jacobian. This is specific
771  /// to solid-based elements in which derivatives w.r.t. to nodal
772  /// positions are evaluated by finite differencing
773  void fill_in_contribution_to_jacobian(Vector<double>& residuals,
774  DenseMatrix<double>& jacobian)
775  {
776  // Call the generic routine
778 
779  // Shape derivatives
780  // Call the generic finite difference routine for the solid variables
781  this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
782  }
783 
784 
785  /// The "global" intrinsic coordinate of the element when
786  /// viewed as part of a geometric object should be given by
787  /// the FaceElement representation, by default
788  double zeta_nodal(const unsigned& n,
789  const unsigned& k,
790  const unsigned& i) const
791  {
792  return FaceElement::zeta_nodal(n, k, i);
793  }
794  };
795 
796 
797  /// //////////////////////////////////////////////////////////////////////
798  /// //////////////////////////////////////////////////////////////////////
799  /// //////////////////////////////////////////////////////////////////////
800 
801 
802  //=======================================================================
803  /// The Two-dimensional interface elements that allow the application of
804  /// a volume constraint specialised for the case when the nodal positions
805  /// of the bulk elements are adjusted using spines.
806  /// To enforce that a fluid volume has a
807  /// certain volume, attach these elements to all faces of the
808  /// (3D Cartesian) bulk fluid elements (of type ELEMENT) that bound that
809  /// region and then specify the "pressure" value that is traded for the
810  /// constraint.
811  //=======================================================================
812  template<class ELEMENT>
815  public virtual SpineElement<FaceGeometry<ELEMENT>>
816  {
817  public:
818  /// Contructor: Specify bulk element and index of face to which
819  /// this face element is to be attached.
821  FiniteElement* const& element_pt, const int& face_index)
822  : SpineElement<FaceGeometry<ELEMENT>>(),
824  {
825  // Attach the geometrical information to the element, by
826  // making the face element from the bulk element
827  element_pt->build_face_element(face_index, this);
828  }
829 
830 
831  /// Fill in contribution to residuals and Jacobian. This is specific
832  /// to spine based elements in which the shape derivatives are evaluated
833  /// using geometric data
834  void fill_in_contribution_to_jacobian(Vector<double>& residuals,
835  DenseMatrix<double>& jacobian)
836  {
837  // Call the generic routine
839 
840  // Call the generic routine to evaluate shape derivatives
841  this->fill_in_jacobian_from_geometric_data(jacobian);
842  }
843 
844 
845  /// The "global" intrinsic coordinate of the element when
846  /// viewed as part of a geometric object should be given by
847  /// the FaceElement representation, by default
848  double zeta_nodal(const unsigned& n,
849  const unsigned& k,
850  const unsigned& i) const
851  {
852  return FaceElement::zeta_nodal(n, k, i);
853  }
854  };
855 
856 
857 } // namespace oomph
858 #endif
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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 contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
double contribution_to_volume_flux()
Return this element's contribution to the volume flux over the boundary.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
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 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(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian. This is specific to solid-based elements in which der...
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
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.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian. This is specific to solid-based elements in which der...
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(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian. This is specific to solid-based elements in which der...
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.
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_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 contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
The axisymmetric (one-dimensional) interface elements that allow imposition of a volume constraint sp...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian. This is specific to spine based elements in which the...
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.
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 ...
The one-dimensional interface elements that allow imposition of a volume constraint specialised for t...
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.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian. This is specific to spine based elements in which the...
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 ...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution to residuals and Jacobian. This is specific to spine based elements in which the...
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_residual_contribution_volume_constraint(Vector< double > &residuals)
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in contribution to residuals and Jacobian.
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...
unsigned Index_of_traded_pressure_value
Index of the value in traded pressure data that corresponds to the traded pressure.
VolumeConstraintBoundingElement()
Constructor initialise the boolean flag We expect the traded pressure data to be stored externally.
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...
bool Traded_pressure_stored_at_node
Boolean to indicate whether the traded pressure is stored externally or at a node (this can happen in...
int ptraded_local_eqn()
The local eqn number for the traded pressure.
unsigned Data_number_of_traded_pressure
The Data that contains the traded pressure is usually stored as external Data for this element,...
A class that is used to implement the constraint that the fluid volume in a region bounded by associa...
void fill_in_generic_contribution_to_residuals_volume_constraint(Vector< double > &residuals)
Fill in the residuals for the volume constraint.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals for the volume constraint.
unsigned Index_of_traded_pressure_value
Index of the value in traded pressure data that corresponds to 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....
Data * p_traded_data_pt()
Access to Data that contains the traded pressure.
VolumeConstraintElement(double *prescribed_volume_pt)
Constructor: Pass pointer to target volume. "Pressure" value that "traded" for the volume contraint i...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the residuals and jacobian for the volume constraint.
double p_traded()
Return the traded pressure value.
bool Traded_pressure_stored_as_internal_data
The Data that contains the traded pressure is stored as external or internal Data for the element....
int ptraded_local_eqn()
The local eqn number for the traded pressure.
unsigned index_of_traded_pressure()
Return the index of Data object at which the traded pressure is stored.
double * Prescribed_volume_pt
Pointer to the desired value of the volume.
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.