constrained_volume_elements.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Header file for elements that allow the imposition of a "constant volume"
// constraint in free surface problems.
 
// Include guards, to prevent multiple includes
#ifndef CONSTRAINED_FLUID_VOLUME_ELEMENTS_HEADER
#define CONSTRAINED_FLUID_VOLUME_ELEMENTS_HEADER
 
// Config header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// OOMPH-LIB headers
#include "generic/Qelements.h"
#include "generic/spines.h"
#include "../axisym_navier_stokes/axisym_navier_stokes_elements.h"
 
//-------------------------------------------
// NOTE: This is still work
// in progress. Need to create versions
// for axisymmetric and cartesian 3D
// bulk equations, as well as spine
// version for all of these (resulting
// in six elements in total). These
// will gradually replace the very hacky
// fix_*.h files that currently live in
// various demo driver codes.
//-------------------------------------------
 
namespace oomph
{
  //==========================================================================
  /// A class that is used to implement the constraint that the fluid volume
  /// in a region bounded by associated FaceElements (attached, e.g., to the
  /// mesh boundaries that enclose a bubble) must take a specific value.
  /// This GeneralisedElement is used only to store the desired volume and
  /// a pointer to the (usually pressure) freedom that must be traded
  /// for the volume constraint.
  //=========================================================================
  class VolumeConstraintElement : public GeneralisedElement
  {
  private:
    /// Pointer to the desired value of the volume
    double* Prescribed_volume_pt;
 
    /// The Data that contains the traded pressure is stored
    /// as external or internal Data for the element. What is its index
    /// in these containers?
    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. Which one?
    bool Traded_pressure_stored_as_internal_data;
 
    /// Index of the value in traded pressure data that corresponds to the
    /// traded pressure
    unsigned Index_of_traded_pressure_value;
 
    /// The local eqn number for the traded pressure
    inline int ptraded_local_eqn()
    {
      if (Traded_pressure_stored_as_internal_data)
      {
        return this->internal_local_eqn(
          External_or_internal_data_index_of_traded_pressure,
          Index_of_traded_pressure_value);
      }
      else
      {
        return this->external_local_eqn(
          External_or_internal_data_index_of_traded_pressure,
          Index_of_traded_pressure_value);
      }
    }
 
 
    /// Fill in the residuals for the volume constraint
    void fill_in_generic_contribution_to_residuals_volume_constraint(
      Vector<double>& residuals);
 
  public:
    /// Constructor: Pass pointer to target volume. "Pressure" value that
    /// "traded" for the volume contraint is created internally (as a Data
    /// item with a single pressure value)
    VolumeConstraintElement(double* prescribed_volume_pt);
 
    /// Constructor: Pass pointer to target volume, pointer to Data
    /// item whose value specified by index_of_traded_pressure represents
    /// the "Pressure" value that "traded" for the volume contraint.
    /// The Data is stored as external Data for this element.
    VolumeConstraintElement(double* prescribed_volume_pt,
                            Data* p_traded_data_pt,
                            const unsigned& index_of_traded_pressure);
 
    /// Empty destructor
    ~VolumeConstraintElement() {}
 
    /// Access to Data that contains the traded pressure
    inline Data* p_traded_data_pt()
    {
      if (Traded_pressure_stored_as_internal_data)
      {
        return internal_data_pt(
          External_or_internal_data_index_of_traded_pressure);
      }
      else
      {
        return external_data_pt(
          External_or_internal_data_index_of_traded_pressure);
      }
    }
 
    /// Return the traded pressure value
    inline double p_traded()
    {
      return p_traded_data_pt()->value(Index_of_traded_pressure_value);
    }
 
    /// Return the index of Data object  at which the traded pressure is stored
    inline unsigned index_of_traded_pressure()
    {
      return Index_of_traded_pressure_value;
    }
 
 
    /// Fill in the residuals for the volume constraint
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      this->fill_in_generic_contribution_to_residuals_volume_constraint(
        residuals);
    }
 
    /// Fill in the residuals and jacobian for the volume constraint
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // No contribution to jacobian; see comment in that function
      this->fill_in_generic_contribution_to_residuals_volume_constraint(
        residuals);
    }
 
    /// Fill in the residuals, jacobian and mass matrix for the volume
    /// constraint.
    void fill_in_contribution_to_jacobian_and_mass_matrix(
      Vector<double>& residuals,
      DenseMatrix<double>& jacobian,
      DenseMatrix<double>& mass_matrix)
    {
      // No contribution to jacobian or mass matrix; see comment in that
      // function
      this->fill_in_generic_contribution_to_residuals_volume_constraint(
        residuals);
    }
  };
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Base class for interface elements that allow the application of
  /// a volume constraint on the region bounded by these elements. The
  /// elements must be used together with the associated
  /// VolumeConstraintElement which stores the value of the
  /// target volume. Common functionality is provided in this base
  /// for storing the external "pressure" value
  /// that is traded for the volume constraint.
  //=======================================================================
  class VolumeConstraintBoundingElement : public virtual FaceElement
  {
  protected:
    /// The Data that contains the traded pressure is usually stored
    /// as external Data for this element, but may also be nodal Data
    /// Which Data item is it?
    unsigned Data_number_of_traded_pressure;
 
    /// Index of the value in traded pressure data that corresponds to the
    /// traded pressure
    unsigned Index_of_traded_pressure_value;
 
    /// Boolean to indicate whether the traded pressure is
    /// stored externally or at a node (this can happen in Taylor-Hood
    /// elements)
    bool Traded_pressure_stored_at_node;
 
    /// The local eqn number for the traded pressure
    inline int ptraded_local_eqn()
    {
      // Return the appropriate nodal value if required
      if (Traded_pressure_stored_at_node)
      {
        return this->nodal_local_eqn(Data_number_of_traded_pressure,
                                     Index_of_traded_pressure_value);
      }
      else
      {
        return this->external_local_eqn(Data_number_of_traded_pressure,
                                        Index_of_traded_pressure_value);
      }
    }
 
    /// Helper function to fill in contributions to residuals
    /// (remember that part of the residual is added by the the
    /// associated VolumeConstraintElement). This is dimension/geometry
    /// specific and must be implemented in derived classes for
    /// 1D line, 2D surface and axisymmetric fluid boundaries
    virtual void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals) = 0;
 
  public:
    /// Constructor initialise the boolean flag
    /// We expect the traded pressure data to be stored externally
    VolumeConstraintBoundingElement() : Traded_pressure_stored_at_node(false) {}
 
    /// Empty Destructor
    ~VolumeConstraintBoundingElement() {}
 
    /// Fill in contribution to residuals and Jacobian
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
    }
 
    /// Set the "master" volume constraint element
    /// The setup here is a bit more complicated than one might expect because
    /// if an internal pressure on a boundary
    /// is hijacked in TaylorHood elements then
    /// that the "traded" value may already be stored as nodal data
    /// in this element. This causes no problems apart from when running
    /// with PARANOID in which case the resulting
    /// repeated local equation numbers are spotted and an error is thrown.
    /// The check is a finite amount of work and so can be avoided if
    /// the boolean flag check_nodal_data is set to false.
    void set_volume_constraint_element(
      VolumeConstraintElement* const& vol_constraint_el_pt,
      const bool& check_nodal_data = true)
    {
      // In order to buffer the case of nodal data, we (tediously) check that
      // the traded pressure is not already nodal data of this element
      if (check_nodal_data)
      {
        // Get memory address of the equation indexed by the
        // traded pressure datum
        long* global_eqn_number =
          vol_constraint_el_pt->p_traded_data_pt()->eqn_number_pt(
            vol_constraint_el_pt->index_of_traded_pressure());
 
        // Put in a check if the datum already corresponds to any other
        //(nodal) data in the element by checking the memory addresses of
        // their global equation numbers
        const unsigned n_node = this->nnode();
        for (unsigned n = 0; n < n_node; n++)
        {
          // Cache the node pointer
          Node* const nod_pt = this->node_pt(n);
          // Find all nodal data values
          unsigned n_value = nod_pt->nvalue();
          // If we already have the data, set the
          // lookup schemes accordingly and return from the function
          for (unsigned i = 0; i < n_value; i++)
          {
            if (nod_pt->eqn_number_pt(i) == global_eqn_number)
            {
              Traded_pressure_stored_at_node = true;
              Data_number_of_traded_pressure = n;
              Index_of_traded_pressure_value = i;
              // Finished so exit the function
              return;
            }
          }
        }
      }
 
      // Should only get here if the data is not nodal
 
      // Add "traded" pressure data as external data to this element
      Data_number_of_traded_pressure =
        this->add_external_data(vol_constraint_el_pt->p_traded_data_pt());
 
      // Which value corresponds to the traded pressure
      Index_of_traded_pressure_value =
        vol_constraint_el_pt->index_of_traded_pressure();
    }
  };
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// One-dimensional interface elements that allow the application of
  /// a volume constraint on the region bounded by these elements.
  /// The volume is computed by integrating x.n around the boundary of
  /// the domain and then dividing by two.
  /// The sign is chosen so that the volume will be positive
  /// when the elements surround a fluid domain.
  ///
  /// These elements must be used together with the associated
  /// VolumeConstraintElement, which stores the value of the
  /// target volume.
  //=======================================================================
  class LineVolumeConstraintBoundingElement
    : public VolumeConstraintBoundingElement
  {
  protected:
    /// Helper function to fill in contributions to residuals
    /// (remember that part of the residual is added by the
    /// the associated VolumeConstraintElement). This is specific for
    /// 1D line elements that bound 2D cartesian fluid elements.
    void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals);
 
  public:
    /// Empty Contructor
    LineVolumeConstraintBoundingElement() : VolumeConstraintBoundingElement() {}
 
    /// Empty Destructor
    ~LineVolumeConstraintBoundingElement() {}
 
    /// Return this element's contribution to the total volume enclosed
    double contribution_to_enclosed_volume();
  };
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// The one-dimensional interface elements that allow imposition of a
  /// volume constraint specialised for the case when the nodal positions of
  /// the bulk elements are treated as solid degrees of freedom.
  /// To enforce that a fluid volume has a
  /// certain volume, attach these elements to all faces of the
  /// (2D cartesian) bulk fluid elements (of type ELEMENT) that bound that
  /// region and then specify the "pressure" value that is traded for the
  /// constraint.
  //=======================================================================
  template<class ELEMENT>
  class ElasticLineVolumeConstraintBoundingElement
    : public LineVolumeConstraintBoundingElement,
      public virtual FaceGeometry<ELEMENT>
  {
  public:
    /// Contructor: Specify bulk element and index of face to which
    /// this face element is to be attached
    ElasticLineVolumeConstraintBoundingElement(FiniteElement* const& element_pt,
                                               const int& face_index)
      : FaceGeometry<ELEMENT>(), LineVolumeConstraintBoundingElement()
    {
      // Attach the geometrical information to the element, by
      // making the face element from the bulk element
      element_pt->build_face_element(face_index, this);
    }
 
    /// Fill in contribution to residuals and Jacobian. This is specific
    /// to solid-based elements in which derivatives w.r.t. to nodal
    /// positions are evaluated by finite differencing
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
 
      // Shape derivatives
      // Call the generic finite difference routine for the solid variables
      this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
    }
 
    /// The "global" intrinsic coordinate of the element when
    /// viewed as part of a geometric object should be given by
    /// the FaceElement representation, by default
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
  //=======================================================================
  /// The one-dimensional interface elements that allow imposition of a
  /// volume constraint specialised for the case when the nodal positions of
  /// the bulk elements are adjusted using Spines.
  /// To enforce that a fluid volume has a
  /// certain volume, attach these elements to all faces of the
  /// (2D cartesian) bulk fluid elements (of type ELEMENT) that bound that
  /// region and then specify the "pressure" value that is traded for the
  /// constraint.
  //=======================================================================
  template<class ELEMENT>
  class SpineLineVolumeConstraintBoundingElement
    : public LineVolumeConstraintBoundingElement,
      public virtual SpineElement<FaceGeometry<ELEMENT>>
  {
  public:
    /// Contructor: Specify bulk element and index of face to which
    /// this face element is to be attached.
    SpineLineVolumeConstraintBoundingElement(FiniteElement* const& element_pt,
                                             const int& face_index)
      : SpineElement<FaceGeometry<ELEMENT>>(),
        LineVolumeConstraintBoundingElement()
    {
      // Attach the geometrical information to the element, by
      // making the face element from the bulk element
      element_pt->build_face_element(face_index, this);
    }
 
 
    /// Fill in contribution to residuals and Jacobian. This is specific
    /// to spine based elements in which the shape derivatives are evaluated
    /// using geometric data
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
 
      // Call the generic routine to evaluate shape derivatives
      this->fill_in_jacobian_from_geometric_data(jacobian);
    }
 
    /// The "global" intrinsic coordinate of the element when
    /// viewed as part of a geometric object should be given by
    /// the FaceElement representation, by default
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Axisymmetric (one-dimensional) interface elements that
  /// allow the application of
  /// a volume constraint on the region bounded by these elements.
  /// The volume is computed by integrating x.n around the boundary of
  /// the domain and then dividing by three.
  /// The sign is chosen so that the volume will be positive
  /// when the elements surround a fluid domain.
  ///
  /// These elements must be used together with the associated
  /// VolumeConstraintElement, which stores the value of the
  /// target volume.
  //=======================================================================
  class AxisymmetricVolumeConstraintBoundingElement
    : public VolumeConstraintBoundingElement
  {
  protected:
    /// Helper function to fill in contributions to residuals
    /// (remember that part of the residual is added by the
    /// the associated VolumeConstraintElement). This is specific for
    /// 1D line elements that bound 2D cartesian fluid elements.
    void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals);
 
  public:
    /// Empty Contructor
    AxisymmetricVolumeConstraintBoundingElement()
      : VolumeConstraintBoundingElement()
    {
    }
 
    /// Empty Destructor
    ~AxisymmetricVolumeConstraintBoundingElement() {}
 
    /// Return this element's contribution to the total volume enclosed
    double contribution_to_enclosed_volume();
 
    /// Return this element's contribution to the volume flux over
    /// the boundary
    double contribution_to_volume_flux()
    {
      // Initialise
      double vol = 0.0;
 
      // Find out how many nodes there are
      const unsigned n_node = this->nnode();
 
      // Set up memeory for the shape functions
      Shape psif(n_node);
      DShape dpsifds(n_node, 1);
 
      // Set the value of n_intpt
      const unsigned n_intpt = this->integral_pt()->nweight();
 
      // Storage for the local cooridinate
      Vector<double> s(1);
 
      // Loop over the integration points
      for (unsigned ipt = 0; ipt < n_intpt; ipt++)
      {
        // Get the local coordinate at the integration point
        s[0] = this->integral_pt()->knot(ipt, 0);
 
        // Get the integral weight
        double W = this->integral_pt()->weight(ipt);
 
        // Call the derivatives of the shape function at the knot point
        this->dshape_local_at_knot(ipt, psif, dpsifds);
 
        // Get position, tangent vector and velocity vector (r and z
        // components only)
        Vector<double> interpolated_u(2, 0.0);
        Vector<double> interpolated_t1(2, 0.0);
        Vector<double> interpolated_x(2, 0.0);
        for (unsigned l = 0; l < n_node; l++)
        {
          // Loop over directional components
          for (unsigned i = 0; i < 2; i++)
          {
            interpolated_x[i] += this->nodal_position(l, i) * psif(l);
            interpolated_u[i] +=
              this->node_pt(l)->value(
                dynamic_cast<AxisymmetricNavierStokesEquations*>(
                  bulk_element_pt())
                  ->u_index_axi_nst(i)) *
              psif(l);
            interpolated_t1[i] += this->nodal_position(l, i) * dpsifds(l, 0);
          }
        }
 
        // Calculate the length of the tangent Vector
        double tlength = interpolated_t1[0] * interpolated_t1[0] +
                         interpolated_t1[1] * interpolated_t1[1];
 
        // Set the Jacobian of the line element
        double J = sqrt(tlength) * interpolated_x[0];
 
        // Now calculate the normal Vector
        Vector<double> interpolated_n(2);
        this->outer_unit_normal(ipt, interpolated_n);
 
        // Assemble dot product
        double dot = 0.0;
        for (unsigned k = 0; k < 2; k++)
        {
          dot += interpolated_u[k] * interpolated_n[k];
        }
 
        // Add to volume with sign chosen so that the volume is
        // positive when the elements bound the fluid
        vol += dot * W * J;
      }
 
      return vol;
    }
  };
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// The axisymmetric (one-dimensional) interface elements that allow
  /// imposition of a
  /// volume constraint specialised for the case when the nodal positions of
  /// the bulk elements are treated as solid degrees of freedom.
  /// To enforce that a fluid volume has a
  /// certain volume, attach these elements to all faces of the
  /// (2D axisymmetric) bulk fluid elements (of type ELEMENT)
  /// that bound that region
  /// and then specify the "pressure" value that is traded for the constraint.
  //=======================================================================
  template<class ELEMENT>
  class ElasticAxisymmetricVolumeConstraintBoundingElement
    : public AxisymmetricVolumeConstraintBoundingElement,
      public virtual FaceGeometry<ELEMENT>
  {
  public:
    /// Contructor: Specify bulk element and index of face to which
    /// this face element is to be attached
    ElasticAxisymmetricVolumeConstraintBoundingElement(
      FiniteElement* const& element_pt, const int& face_index)
      : FaceGeometry<ELEMENT>(), AxisymmetricVolumeConstraintBoundingElement()
    {
      // Attach the geometrical information to the element, by
      // making the face element from the bulk element
      element_pt->build_face_element(face_index, this);
    }
 
    /// Fill in contribution to residuals and Jacobian. This is specific
    /// to solid-based elements in which derivatives w.r.t. to nodal
    /// positions are evaluated by finite differencing
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
 
      // Shape derivatives
      // Call the generic finite difference routine for the solid variables
      this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
    }
 
    /// The "global" intrinsic coordinate of the element when
    /// viewed as part of a geometric object should be given by
    /// the FaceElement representation, by default
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
  //=======================================================================
  /// The axisymmetric (one-dimensional) interface elements that
  /// allow imposition of a
  /// volume constraint specialised for the case when the nodal positions of
  /// the bulk elements are adjusted using Spines.
  /// To enforce that a fluid volume has a
  /// certain volume, attach these elements to all faces of the
  /// (2D axisymmetric) bulk fluid elements (of type ELEMENT) that bound
  /// that region
  /// and then specify the "pressure" value that is traded for the constraint.
  //=======================================================================
  template<class ELEMENT>
  class SpineAxisymmetricVolumeConstraintBoundingElement
    : public AxisymmetricVolumeConstraintBoundingElement,
      public virtual SpineElement<FaceGeometry<ELEMENT>>
  {
  public:
    /// Contructor: Specify bulk element and index of face to which
    /// this face element is to be attached.
    SpineAxisymmetricVolumeConstraintBoundingElement(
      FiniteElement* const& element_pt, const int& face_index)
      : SpineElement<FaceGeometry<ELEMENT>>(),
        AxisymmetricVolumeConstraintBoundingElement()
    {
      // Attach the geometrical information to the element, by
      // making the face element from the bulk element
      element_pt->build_face_element(face_index, this);
    }
 
 
    /// Fill in contribution to residuals and Jacobian. This is specific
    /// to spine based elements in which the shape derivatives are evaluated
    /// using geometric data
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
 
      // Call the generic routine to evaluate shape derivatives
      this->fill_in_jacobian_from_geometric_data(jacobian);
    }
 
    /// The "global" intrinsic coordinate of the element when
    /// viewed as part of a geometric object should be given by
    /// the FaceElement representation, by default
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// Two-dimensional interface elements that allow the application of
  /// a volume constraint on the region bounded by these elements.
  /// The volume is computed by integrating x.n around the boundary of
  /// the domain and then dividing by three.
  /// The sign is chosen so that the volume will be positive
  /// when the elements surround a fluid domain.
  ///
  /// These elements must be used together with the associated
  /// VolumeConstraintElement, which stores the value of the
  /// target volume.
  //=======================================================================
  class SurfaceVolumeConstraintBoundingElement
    : public VolumeConstraintBoundingElement
  {
  protected:
    /// Helper function to fill in contributions to residuals
    /// (remember that part of the residual is added by the
    /// the associated VolumeConstraintElement). This is specific for
    /// 2D surface elements that bound 3D cartesian fluid elements.
    void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals);
 
  public:
    /// Empty Contructor
    SurfaceVolumeConstraintBoundingElement() : VolumeConstraintBoundingElement()
    {
    }
 
    /// Empty Desctructor
    ~SurfaceVolumeConstraintBoundingElement() {}
  };
 
 
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// The Two-dimensional interface elements that allow the application of
  /// a volume constraint specialised for the case when the nodal positions
  /// of the bulk elements are treated as solid degrees of freedom.
  /// To enforce that a fluid volume has a
  /// certain volume, attach these elements to all faces of the
  /// (3D Cartesian) bulk fluid elements (of type ELEMENT) that bound that
  /// region and then specify the "pressure" value that is traded for the
  /// constraint.
  //=======================================================================
  template<class ELEMENT>
  class ElasticSurfaceVolumeConstraintBoundingElement
    : public SurfaceVolumeConstraintBoundingElement,
      public virtual FaceGeometry<ELEMENT>
  {
  public:
    /// Contructor: Specify bulk element and index of face to which
    /// this face element is to be attached.
    ElasticSurfaceVolumeConstraintBoundingElement(
      FiniteElement* const& element_pt, const int& face_index)
      : FaceGeometry<ELEMENT>(), SurfaceVolumeConstraintBoundingElement()
    {
      // Attach the geometrical information to the element, by
      // making the face element from the bulk element
      element_pt->build_face_element(face_index, this);
    }
 
 
    /// Fill in contribution to residuals and Jacobian. This is specific
    /// to solid-based elements in which derivatives w.r.t. to nodal
    /// positions are evaluated by finite differencing
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
 
      // Shape derivatives
      // Call the generic finite difference routine for the solid variables
      this->fill_in_jacobian_from_solid_position_by_fd(jacobian);
    }
 
 
    /// The "global" intrinsic coordinate of the element when
    /// viewed as part of a geometric object should be given by
    /// the FaceElement representation, by default
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
 
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////
 
 
  //=======================================================================
  /// The Two-dimensional interface elements that allow the application of
  /// a volume constraint specialised for the case when the nodal positions
  /// of the bulk elements are adjusted using spines.
  /// To enforce that a fluid volume has a
  /// certain volume, attach these elements to all faces of the
  /// (3D Cartesian) bulk fluid elements (of type ELEMENT) that bound that
  /// region and then specify the "pressure" value that is traded for the
  /// constraint.
  //=======================================================================
  template<class ELEMENT>
  class SpineSurfaceVolumeConstraintBoundingElement
    : public SurfaceVolumeConstraintBoundingElement,
      public virtual SpineElement<FaceGeometry<ELEMENT>>
  {
  public:
    /// Contructor: Specify bulk element and index of face to which
    /// this face element is to be attached.
    SpineSurfaceVolumeConstraintBoundingElement(
      FiniteElement* const& element_pt, const int& face_index)
      : SpineElement<FaceGeometry<ELEMENT>>(),
        SurfaceVolumeConstraintBoundingElement()
    {
      // Attach the geometrical information to the element, by
      // making the face element from the bulk element
      element_pt->build_face_element(face_index, this);
    }
 
 
    /// Fill in contribution to residuals and Jacobian. This is specific
    /// to spine based elements in which the shape derivatives are evaluated
    /// using geometric data
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
 
      // Call the generic routine to evaluate shape derivatives
      this->fill_in_jacobian_from_geometric_data(jacobian);
    }
 
 
    /// The "global" intrinsic coordinate of the element when
    /// viewed as part of a geometric object should be given by
    /// the FaceElement representation, by default
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
 
} // namespace oomph
#endif