poroelasticity_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//====================================================================
#ifndef OOMPH_POROELASTICITY_ELEMENTS_HEADER
#define OOMPH_POROELASTICITY_ELEMENTS_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "../generic/elements.h"
#include "../generic/shape.h"
 
#include "elasticity_tensor.h"
 
namespace oomph
{
  template<unsigned DIM>
  class PoroelasticityEquations : public virtual FiniteElement
  {
  public:
    typedef void (*SourceFctPt)(const double& time,
                                const Vector<double>& x,
                                Vector<double>& f);
 
    typedef void (*MassSourceFctPt)(const double& time,
                                    const Vector<double>& x,
                                    double& f);
 
    PoroelasticityEquations()
      : Elasticity_tensor_pt(0),
        Force_solid_fct_pt(0),
        Force_fluid_fct_pt(0),
        Mass_source_fct_pt(0),
        Lambda_sq_pt(&Default_lambda_sq_value),
        Density_ratio_pt(&Default_density_ratio_value),
        K_inv_pt(&Default_k_inv_value),
        Alpha_pt(&Default_alpha_value),
        Porosity_pt(&Default_porosity_value)
    {
    }
 
    const double& lambda_sq() const
    {
      return *Lambda_sq_pt;
    }
 
    double*& lambda_sq_pt()
    {
      return Lambda_sq_pt;
    }
 
    const double& density_ratio() const
    {
      return *Density_ratio_pt;
    }
 
    double*& density_ratio_pt()
    {
      return Density_ratio_pt;
    }
 
    const double& k_inv() const
    {
      return *K_inv_pt;
    }
 
    double*& k_inv_pt()
    {
      return K_inv_pt;
    }
 
    const double& alpha() const
    {
      return *Alpha_pt;
    }
 
    double*& alpha_pt()
    {
      return Alpha_pt;
    }
 
    const double& porosity() const
    {
      return *Porosity_pt;
    }
 
    double*& porosity_pt()
    {
      return Porosity_pt;
    }
 
    SourceFctPt& force_solid_fct_pt()
    {
      return Force_solid_fct_pt;
    }
 
    SourceFctPt force_solid_fct_pt() const
    {
      return Force_solid_fct_pt;
    }
 
    SourceFctPt& force_fluid_fct_pt()
    {
      return Force_fluid_fct_pt;
    }
 
    SourceFctPt force_fluid_fct_pt() const
    {
      return Force_fluid_fct_pt;
    }
 
    MassSourceFctPt& mass_source_fct_pt()
    {
      return Mass_source_fct_pt;
    }
 
    MassSourceFctPt mass_source_fct_pt() const
    {
      return Mass_source_fct_pt;
    }
 
    void force_solid(const double& time,
                     const Vector<double>& x,
                     Vector<double>& b) const
    {
      // If no function has been set, return zero vector
      if (Force_solid_fct_pt == 0)
      {
        // Get spatial dimension of element
        unsigned n = dim();
        for (unsigned i = 0; i < n; i++)
        {
          b[i] = 0.0;
        }
      }
      else
      {
        (*Force_solid_fct_pt)(time, x, b);
      }
    }
 
    void force_fluid(const double& time,
                     const Vector<double>& x,
                     Vector<double>& b) const
    {
      // If no function has been set, return zero vector
      if (Force_fluid_fct_pt == 0)
      {
        // Get spatial dimension of element
        unsigned n = dim();
        for (unsigned i = 0; i < n; i++)
        {
          b[i] = 0.0;
        }
      }
      else
      {
        (*Force_fluid_fct_pt)(time, x, b);
      }
    }
 
    void mass_source(const double& time,
                     const Vector<double>& x,
                     double& b) const
    {
      // If no function has been set, return zero vector
      if (Mass_source_fct_pt == 0)
      {
        b = 0.0;
      }
      else
      {
        (*Mass_source_fct_pt)(time, x, b);
      }
    }
 
    ElasticityTensor*& elasticity_tensor_pt()
    {
      return Elasticity_tensor_pt;
    }
 
    const double E(const unsigned& i,
                   const unsigned& j,
                   const unsigned& k,
                   const unsigned& l) const
    {
      return (*Elasticity_tensor_pt)(i, j, k, l);
    }
 
    void get_stress(const Vector<double>& s, DenseMatrix<double>& sigma) const;
 
    void get_strain(const Vector<double>& s, DenseMatrix<double>& strain) const;
 
    virtual unsigned required_nvalue(const unsigned& n) const = 0;
 
    virtual unsigned u_index(const unsigned& n) const = 0;
 
    virtual int q_edge_local_eqn(const unsigned& n) const = 0;
 
    virtual int q_internal_local_eqn(const unsigned& n) const = 0;
 
    virtual unsigned q_edge_index(const unsigned& n) const = 0;
 
    virtual unsigned q_internal_index() const = 0;
 
    virtual unsigned q_edge_node_number(const unsigned& n) const = 0;
 
    virtual double q_edge(const unsigned& n) const = 0;
 
    virtual double q_edge(const unsigned& t, const unsigned& n) const = 0;
 
    virtual double q_internal(const unsigned& n) const = 0;
 
    virtual double q_internal(const unsigned& t, const unsigned& n) const = 0;
 
    virtual unsigned nq_basis() const = 0;
 
    virtual unsigned nq_basis_edge() const = 0;
 
    virtual void get_q_basis_local(const Vector<double>& s,
                                   Shape& q_basis) const = 0;
 
    virtual void get_div_q_basis_local(const Vector<double>& s,
                                       Shape& div_q_basis_ds) const = 0;
 
    void get_q_basis(const Vector<double>& s, Shape& q_basis) const
    {
      const unsigned n_node = this->nnode();
      Shape psi(n_node, DIM);
      const unsigned n_q_basis = this->nq_basis();
      Shape q_basis_local(n_q_basis, DIM);
      this->get_q_basis_local(s, q_basis_local);
      (void)this->transform_basis(s, q_basis_local, psi, q_basis);
    }
 
    virtual unsigned nedge_gauss_point() const = 0;
 
    virtual double edge_gauss_point(const unsigned& edge,
                                    const unsigned& n) const = 0;
 
    virtual void edge_gauss_point_global(const unsigned& edge,
                                         const unsigned& n,
                                         Vector<double>& x) const = 0;
 
    virtual void pin_q_internal_value(const unsigned& n) = 0;
 
    virtual int p_local_eqn(const unsigned& n) const = 0;
 
    virtual double p_value(unsigned& n) const = 0;
 
    virtual unsigned np_basis() const = 0;
 
    virtual void get_p_basis(const Vector<double>& s, Shape& p_basis) const = 0;
 
    virtual void pin_p_value(const unsigned& n, const double& p) = 0;
 
    virtual void scale_basis(Shape& basis) const = 0;
 
    double transform_basis(const Vector<double>& s,
                           const Shape& q_basis_local,
                           Shape& psi,
                           DShape& dpsi,
                           Shape& q_basis) const;
 
    double transform_basis(const Vector<double>& s,
                           const Shape& q_basis_local,
                           Shape& psi,
                           Shape& q_basis) const
    {
      const unsigned n_node = this->nnode();
      DShape dpsi(n_node, DIM);
      return transform_basis(s, q_basis_local, psi, dpsi, q_basis);
    }
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      this->fill_in_generic_residual_contribution(
        residuals, GeneralisedElement::Dummy_matrix, 0);
    }
 
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      this->fill_in_generic_residual_contribution(residuals, jacobian, 1);
    }
 
    void interpolated_u(const Vector<double>& s, Vector<double>& disp) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      for (unsigned i = 0; i < DIM; i++)
      {
        // Index at which the nodal value is stored
        unsigned u_nodal_index = u_index(i);
 
        // Initialise value of u
        disp[i] = 0.0;
 
        // Loop over the local nodes and sum
        for (unsigned l = 0; l < n_node; l++)
        {
          disp[i] += nodal_value(l, u_nodal_index) * psi[l];
        }
      }
    }
 
    double interpolated_u(const Vector<double>& s, const unsigned& i) const
    {
      // Find number of nodes
      unsigned n_node = nnode();
 
      // Local shape function
      Shape psi(n_node);
 
      // Find values of shape function
      shape(s, psi);
 
      // Get nodal index at which i-th velocity is stored
      unsigned u_nodal_index = u_index(i);
 
      // Initialise value of u
      double interpolated_u = 0.0;
 
      // Loop over the local nodes and sum
      for (unsigned l = 0; l < n_node; l++)
      {
        interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
      }
 
      return (interpolated_u);
    }
 
    void interpolated_q(const Vector<double>& s, Vector<double>& u) const
    {
      unsigned n_q_basis = nq_basis();
      unsigned n_q_basis_edge = nq_basis_edge();
 
      Shape q_basis(n_q_basis, DIM);
 
      get_q_basis(s, q_basis);
      for (unsigned i = 0; i < DIM; i++)
      {
        u[i] = 0.0;
        for (unsigned l = 0; l < n_q_basis_edge; l++)
        {
          u[i] += q_edge(l) * q_basis(l, i);
        }
        for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
        {
          u[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
        }
      }
    }
 
    double interpolated_q(const Vector<double>& s, const unsigned i) const
    {
      unsigned n_q_basis = nq_basis();
      unsigned n_q_basis_edge = nq_basis_edge();
 
      Shape q_basis(n_q_basis, DIM);
 
      get_q_basis(s, q_basis);
      double q_i = 0.0;
      for (unsigned l = 0; l < n_q_basis_edge; l++)
      {
        q_i += q_edge(l) * q_basis(l, i);
      }
      for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
      {
        q_i += q_internal(l - n_q_basis_edge) * q_basis(l, i);
      }
 
      return q_i;
    }
 
    void interpolated_div_q(const Vector<double>& s, double& div_q) const
    {
      // Zero the divergence
      div_q = 0;
 
      // Get the number of nodes, q basis function, and q edge basis functions
      unsigned n_node = nnode();
      const unsigned n_q_basis = nq_basis();
      const unsigned n_q_basis_edge = nq_basis_edge();
 
      // Storage for the divergence basis
      Shape div_q_basis_ds(n_q_basis);
 
      // Storage for the geometric basis and it's derivatives
      Shape psi(n_node);
      DShape dpsi(n_node, DIM);
 
      // Call the geometric shape functions and their derivatives
      this->dshape_local(s, psi, dpsi);
 
      // Storage for the inverse of the geometric jacobian (just so we can call
      // the local to eulerian mapping)
      DenseMatrix<double> inverse_jacobian(DIM);
 
      // Get the determinant of the geometric mapping
      double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
 
      // Get the divergence basis (wrt local coords) at local coords s
      get_div_q_basis_local(s, div_q_basis_ds);
 
      // Add the contribution to the divergence from the edge basis functions
      for (unsigned l = 0; l < n_q_basis_edge; l++)
      {
        div_q += 1.0 / det * div_q_basis_ds(l) * q_edge(l);
      }
 
      // Add the contribution to the divergence from the internal basis
      // functions
      for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
      {
        div_q += 1.0 / det * div_q_basis_ds(l) * q_internal(l - n_q_basis_edge);
      }
    }
 
    double interpolated_div_q(const Vector<double>& s)
    {
      // Temporary storage for div u
      double div_q = 0;
 
      // Get the intepolated divergence
      interpolated_div_q(s, div_q);
 
      // Return it
      return div_q;
    }
 
    void interpolated_p(const Vector<double>& s, double& p) const
    {
      // Get the number of p basis functions
      unsigned n_p_basis = np_basis();
 
      // Storage for the p basis
      Shape p_basis(n_p_basis);
 
      // Call the p basis
      get_p_basis(s, p_basis);
 
      // Zero the pressure
      p = 0;
 
      // Add the contribution to the pressure from each basis function
      for (unsigned l = 0; l < n_p_basis; l++)
      {
        p += p_value(l) * p_basis(l);
      }
    }
 
    double interpolated_p(const Vector<double>& s) const
    {
      // Temporary storage for p
      double p = 0;
 
      // Get the interpolated pressure
      interpolated_p(s, p);
 
      // Return it
      return p;
    }
 
    double du_dt(const unsigned& n, const unsigned& i) const
    {
      // Get the timestepper
      TimeStepper* time_stepper_pt = node_pt(n)->time_stepper_pt();
 
      // Storage for the derivative - initialise to 0
      double du_dt = 0.0;
 
      // If we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // Get the nodal index
        const unsigned u_nodal_index = u_index(i);
 
        // Get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // Loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // Add the contribution to the derivative
          du_dt +=
            time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
        }
      }
 
      return du_dt;
    }
 
    double d2u_dt2(const unsigned& n, const unsigned& i) const
    {
      // Get the timestepper
      TimeStepper* time_stepper_pt = node_pt(n)->time_stepper_pt();
 
      // Storage for the derivative - initialise to 0
      double d2u_dt2 = 0.0;
 
      // If we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // Get the nodal index
        const unsigned u_nodal_index = u_index(i);
 
        // Get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // Loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // Add the contribution to the derivative
          d2u_dt2 +=
            time_stepper_pt->weight(2, t) * nodal_value(t, n, u_nodal_index);
        }
      }
 
      return d2u_dt2;
    }
 
    double dq_edge_dt(const unsigned& n) const
    {
      unsigned node_num = q_edge_node_number(n);
 
      // get the timestepper
      TimeStepper* time_stepper_pt = node_pt(node_num)->time_stepper_pt();
 
      // storage for the derivative - initialise to 0
      double dq_dt = 0.0;
 
      // if we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // add the contribution to the derivative
          dq_dt += time_stepper_pt->weight(1, t) * q_edge(t, n);
        }
      }
 
      return dq_dt;
    }
 
    double dq_internal_dt(const unsigned& n) const
    {
      // get the internal data index for q
      unsigned internal_index = q_internal_index();
 
      // get the timestepper
      TimeStepper* time_stepper_pt =
        internal_data_pt(internal_index)->time_stepper_pt();
 
      // storage for the derivative - initialise to 0
      double dq_dt = 0.0;
 
      // if we are doing an unsteady solve then calculate the derivative
      if (!time_stepper_pt->is_steady())
      {
        // get the number of values required to represent history
        const unsigned n_time = time_stepper_pt->ntstorage();
 
        // loop over history values
        for (unsigned t = 0; t < n_time; t++)
        {
          // add the contribution to the derivative
          dq_dt += time_stepper_pt->weight(1, t) * q_internal(t, n);
        }
      }
 
      return dq_dt;
    }
 
    void set_q_internal_timestepper(TimeStepper* const time_stepper_pt)
    {
      unsigned q_index = q_internal_index();
 
      this->internal_data_pt(q_index)->set_time_stepper(time_stepper_pt, false);
    }
 
    unsigned self_test()
    {
      return 0;
    }
 
    void output(std::ostream& outfile)
    {
      unsigned nplot = 5;
      output(outfile, nplot);
    }
 
    void output(std::ostream& outfile, const unsigned& nplot);
 
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    FiniteElement::SteadyExactSolutionFctPt exact_soln_pt);
 
    void output_fct(std::ostream& outfile,
                    const unsigned& nplot,
                    const double& time,
                    FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt);
 
    void compute_error(std::ostream& outfile,
                       FiniteElement::SteadyExactSolutionFctPt exact_soln_pt,
                       Vector<double>& error,
                       Vector<double>& norm);
 
    void compute_error(std::ostream& outfile,
                       FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt,
                       const double& time,
                       Vector<double>& error,
                       Vector<double>& norm);
 
  protected:
    virtual double shape_basis_test_local(const Vector<double>& s,
                                          Shape& psi,
                                          DShape& dpsi,
                                          Shape& u_basis,
                                          Shape& u_test,
                                          DShape& du_basis_dx,
                                          DShape& du_test_dx,
                                          Shape& q_basis,
                                          Shape& q_test,
                                          Shape& p_basis,
                                          Shape& p_test,
                                          Shape& div_q_basis_ds,
                                          Shape& div_q_test_ds) const = 0;
 
    virtual double shape_basis_test_local_at_knot(
      const unsigned& ipt,
      Shape& psi,
      DShape& dpsi,
      Shape& u_basis,
      Shape& u_test,
      DShape& du_basis_dx,
      DShape& du_test_dx,
      Shape& q_basis,
      Shape& q_test,
      Shape& p_basis,
      Shape& p_test,
      Shape& div_q_basis_ds,
      Shape& div_q_test_ds) const = 0;
 
    // fill in residuals and, if flag==true, jacobian
    virtual void fill_in_generic_residual_contribution(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, bool flag);
 
    ElasticityTensor* Elasticity_tensor_pt;
 
  private:
    SourceFctPt Force_solid_fct_pt;
 
    SourceFctPt Force_fluid_fct_pt;
 
    MassSourceFctPt Mass_source_fct_pt;
 
    double* Lambda_sq_pt;
 
    double* Density_ratio_pt;
 
    double* K_inv_pt;
 
    double* Alpha_pt;
 
    double* Porosity_pt;
 
    static double Default_lambda_sq_value;
 
    static double Default_density_ratio_value;
 
    static double Default_k_inv_value;
 
    static double Default_alpha_value;
 
    static double Default_porosity_value;
  };
 
} // namespace oomph
 
#endif