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-2022 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 generated by autoconfig
#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
{
  //==========================================================================
  //=========================================================================
  class VolumeConstraintElement : public GeneralisedElement
  {
  private:
    double* Prescribed_volume_pt;
 
    unsigned External_or_internal_data_index_of_traded_pressure;
 
    bool Traded_pressure_stored_as_internal_data;
 
    unsigned Index_of_traded_pressure_value;
 
    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);
      }
    }
 
 
    void fill_in_generic_contribution_to_residuals_volume_constraint(
      Vector<double>& residuals);
 
  public:
    VolumeConstraintElement(double* prescribed_volume_pt);
 
    VolumeConstraintElement(double* prescribed_volume_pt,
                            Data* p_traded_data_pt,
                            const unsigned& index_of_traded_pressure);
 
    ~VolumeConstraintElement() {}
 
    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);
      }
    }
 
    inline double p_traded()
    {
      return p_traded_data_pt()->value(Index_of_traded_pressure_value);
    }
 
    inline unsigned index_of_traded_pressure()
    {
      return Index_of_traded_pressure_value;
    }
 
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      this->fill_in_generic_contribution_to_residuals_volume_constraint(
        residuals);
    }
 
    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);
    }
 
    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);
    }
  };
 
 
 
  //=======================================================================
  //=======================================================================
  class VolumeConstraintBoundingElement : public virtual FaceElement
  {
  protected:
    unsigned Data_number_of_traded_pressure;
 
    unsigned Index_of_traded_pressure_value;
 
    bool Traded_pressure_stored_at_node;
 
    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);
      }
    }
 
    virtual void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals) = 0;
 
  public:
    VolumeConstraintBoundingElement() : Traded_pressure_stored_at_node(false) {}
 
    ~VolumeConstraintBoundingElement() {}
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      // Call the generic routine
      this->fill_in_generic_residual_contribution_volume_constraint(residuals);
    }
 
    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();
    }
  };
 
 
 
  //=======================================================================
  //=======================================================================
  class LineVolumeConstraintBoundingElement
    : public VolumeConstraintBoundingElement
  {
  protected:
    void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals);
 
  public:
    LineVolumeConstraintBoundingElement() : VolumeConstraintBoundingElement() {}
 
    ~LineVolumeConstraintBoundingElement() {}
 
    double contribution_to_enclosed_volume();
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class ElasticLineVolumeConstraintBoundingElement
    : public LineVolumeConstraintBoundingElement,
      public virtual FaceGeometry<ELEMENT>
  {
  public:
    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);
    }
 
    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);
    }
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class SpineLineVolumeConstraintBoundingElement
    : public LineVolumeConstraintBoundingElement,
      public virtual SpineElement<FaceGeometry<ELEMENT>>
  {
  public:
    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);
    }
 
 
    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);
    }
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
 
 
  //=======================================================================
  //=======================================================================
  class AxisymmetricVolumeConstraintBoundingElement
    : public VolumeConstraintBoundingElement
  {
  protected:
    void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals);
 
  public:
    AxisymmetricVolumeConstraintBoundingElement()
      : VolumeConstraintBoundingElement()
    {
    }
 
    ~AxisymmetricVolumeConstraintBoundingElement() {}
 
    double contribution_to_enclosed_volume();
 
    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;
    }
  };
 
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class ElasticAxisymmetricVolumeConstraintBoundingElement
    : public AxisymmetricVolumeConstraintBoundingElement,
      public virtual FaceGeometry<ELEMENT>
  {
  public:
    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);
    }
 
    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);
    }
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class SpineAxisymmetricVolumeConstraintBoundingElement
    : public AxisymmetricVolumeConstraintBoundingElement,
      public virtual SpineElement<FaceGeometry<ELEMENT>>
  {
  public:
    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);
    }
 
 
    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);
    }
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  class SurfaceVolumeConstraintBoundingElement
    : public VolumeConstraintBoundingElement
  {
  protected:
    void fill_in_generic_residual_contribution_volume_constraint(
      Vector<double>& residuals);
 
  public:
    SurfaceVolumeConstraintBoundingElement() : VolumeConstraintBoundingElement()
    {
    }
 
    ~SurfaceVolumeConstraintBoundingElement() {}
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class ElasticSurfaceVolumeConstraintBoundingElement
    : public SurfaceVolumeConstraintBoundingElement,
      public virtual FaceGeometry<ELEMENT>
  {
  public:
    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);
    }
 
 
    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);
    }
 
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<class ELEMENT>
  class SpineSurfaceVolumeConstraintBoundingElement
    : public SurfaceVolumeConstraintBoundingElement,
      public virtual SpineElement<FaceGeometry<ELEMENT>>
  {
  public:
    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);
    }
 
 
    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);
    }
 
 
    double zeta_nodal(const unsigned& n,
                      const unsigned& k,
                      const unsigned& i) const
    {
      return FaceElement::zeta_nodal(n, k, i);
    }
  };
 
 
} // namespace oomph
#endif