homo_lin_elasticity_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 general solid mechanics elements
 
//Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_HOMOGENISED_LINEAR_ELASTICITY_ELEMENTS_HEADER
#define OOMPH_HOMOGENISED_LINEAR_ELASTICITY_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/mesh.h"
#include "generic/block_preconditioner.h"
#include "linear_elasticity/elasticity_tensor.h"
 
 
namespace oomph
{
//=======================================================================
//=======================================================================
  class HomogenisedLinearElasticityEquationsBase : 
   public virtual FiniteElement
  {
    public:
 
   typedef void (*HomLinElasticityTensorFctPt)(const Vector<double> &x,
                                               ElasticityTensor* &E_pt);
   
   virtual inline unsigned u_index_linear_elasticity(const unsigned i) const
    {return i;}
   
   void interpolated_u_linear_elasticity(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<3;i++)
      {
       //Index at which the nodal value is stored
       unsigned u_nodal_index = u_index_linear_elasticity(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_linear_elasticity(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_linear_elasticity(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);
    }
   
   HomogenisedLinearElasticityEquationsBase() : Elasticity_tensor_fct_pt(0),
    Lambda_sq_pt(&Default_lambda_sq_value), P_pt(0), M_pt(0), Unsteady(false)
    {}
 
   unsigned* &p_pt() {return P_pt;}
   
   unsigned get_p()
   {
    if(P_pt==0) 
     {
      throw OomphLibError("Integer P pointer not set\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      return 0;
     }
    else {return *P_pt;}
   }
 
   unsigned* &m_pt() {return M_pt;}
 
   unsigned get_m()
   {
    if(M_pt==0) 
     {
      throw OomphLibError("Integer M pointer not set\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
      return 0;
     }
    else {return *M_pt;}
   }
 
   HomLinElasticityTensorFctPt& elasticity_tensor_fct_pt() 
    {return  Elasticity_tensor_fct_pt;}
   
   HomLinElasticityTensorFctPt elasticity_tensor_fct_pt() const
    {return  Elasticity_tensor_fct_pt;}
 
   inline void get_E_pt(const Vector<double> &x, ElasticityTensor*  &E_pt)
   {
    if(this->Elasticity_tensor_fct_pt==0)
     {
      throw OomphLibError("Elasticity tensor function pointer not set\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
    else
     {
      (*this->Elasticity_tensor_fct_pt)(x,E_pt);
     }
   }
 
   const double& lambda_sq() const {return *Lambda_sq_pt;}
   
   double* &lambda_sq_pt() {return Lambda_sq_pt;}
   
   bool& unsteady() {return Unsteady;}
   
   bool unsteady() const {return Unsteady;}
 
   virtual FaceElement* make_face_element(const unsigned &s_fixed_index,
                                          const int&s_limit)
    {
     throw OomphLibError(
      "make_face_element not implemented",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
  protected:
   
   HomLinElasticityTensorFctPt Elasticity_tensor_fct_pt;
 
   double* Lambda_sq_pt;
   
   unsigned *P_pt;
   
   unsigned *M_pt;
   
   bool Unsteady;
 
   static double Default_lambda_sq_value;
   
  };
 
 
//=======================================================================
//=======================================================================
template<unsigned DIM>
 class HomogenisedLinearElasticityEquations : 
   public HomogenisedLinearElasticityEquationsBase
   {
    public:
    
   HomogenisedLinearElasticityEquations() {}
 
   unsigned required_nvalue(const unsigned &n) const {return 3;}
   
   void fill_in_contribution_to_residuals(Vector<double> &residuals)
    {
     fill_in_generic_contribution_to_residuals_linear_elasticity(
      residuals,GeneralisedElement::Dummy_matrix,0);
    }
 
   virtual void fill_in_generic_contribution_to_residuals_linear_elasticity(
    Vector<double> &residuals, DenseMatrix<double> &jacobian, unsigned flag);
   
   void fill_in_contribution_to_jacobian(Vector<double> &residuals,
                                         DenseMatrix<double> &jacobian)
    {
     //Add the contribution to the residuals
     this->fill_in_generic_contribution_to_residuals_linear_elasticity(
      residuals,jacobian,1);
    }
 
   void calculate_effective_modulus(DenseMatrix<double> &H);
   
   void output(std::ostream &outfile) 
    {
     unsigned n_plot=5;
     output(outfile,n_plot);
    }
   
   void output(std::ostream &outfile, const unsigned &n_plot);
   
   
   void output(FILE* file_pt) 
    {
     unsigned n_plot=5;
     output(file_pt,n_plot);
    }
   
   void output(FILE* file_pt, const unsigned &n_plot);
   
  }; 
 
 
}
 
#endif