mechanic_integrators.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#ifndef _MECH_INTEGRATORS_H
#define _MECH_INTEGRATORS_H
 
#include "SF_base.h"
#include "basics.h"
#include "fem_types.h"
#include "sf_interface.h"
#include "constitutive_model_library.h"
#include "constitutive_model.h"
#include "mechanics.h"
 
namespace opencarp {
 
#define EP_MAX_IPOINTS 128
#define EP_MAX_LPOINTS 8
 
class mech_stiffness_integrator : public SF::matrix_integrator<mesh_int_t,mesh_real_t>
{
  MaterialType & material;   
 
  SF::dmat<double> cond;          
  SF::dmat<double> shape;         
  SF::dmat<double> gshape;        
  double J[9];                    
  double F[9];                    
  double cauchy_stress[6];        
  SF::Point lpts[EP_MAX_LPOINTS]; 
 
  int p_order = 1;  
 
  SF::Point           ipts[EP_MAX_IPOINTS];    
  double              w   [EP_MAX_IPOINTS];    
 
  public:
  mech_stiffness_integrator(MaterialType & inp_mat) : material(inp_mat)
  {}
  void operator() (const SF::element_view<mesh_int_t, mesh_real_t> & elem, SF::dmat<double> & buff);
  void dpn(mesh_int_t & row_dpn, mesh_int_t & col_dpn);
};
 
class mech_internal_forces_integrator : public SF::vector_integrator<mesh_int_t, mesh_real_t>
{
  MaterialType & material;   
    
  SF::dmat<double> shape;         
  SF::dmat<double> gshape;        
  double J[9];                    
  double F[9];                    
  double iF[9];                   
  double S_PK1[9];                
  double S_PK1_xyz[9];            
  SF::Point lpts[EP_MAX_LPOINTS]; 
    
  int p_order = 1; 
    
  SF::Point ipts[EP_MAX_IPOINTS]; 
  double w[EP_MAX_IPOINTS]; 
    
  public:
  mech_internal_forces_integrator(MaterialType & inp_mat) : material(inp_mat)
  {}
  void operator() (const SF::element_view<mesh_int_t, mesh_real_t> & elem, double * buff);
  void dpn(mesh_int_t & dpn);
};
 
 
class mech_mass_integrator : public SF::matrix_integrator<mesh_int_t,mesh_real_t>
{
  SF::dmat<double> shape;         
  double J[9];                    
  SF::Point lpts[EP_MAX_LPOINTS]; 
 
  int p_order = 1;  
 
  SF::Point           ipts[EP_MAX_IPOINTS];    
  double              w   [EP_MAX_IPOINTS];    
 
  public:  void operator() (const SF::element_view<mesh_int_t, mesh_real_t> & elem, SF::dmat<double> & buff);
  void dpn(mesh_int_t & row_dpn, mesh_int_t & col_dpn);
};
 
void calc_deformation_gradient(const SF::element_view<mesh_int_t, mesh_real_t> & elem, int ipt, double * F, double * iF, double & detF, int index = -1, double eps = 0.0);
 
 
inline void calc_left_cauchy_green(double * F, double * b)
{
  b[0] = F[0]*F[0] + F[3]*F[3] + F[6]*F[6];
  b[1] = F[0]*F[1] + F[3]*F[4] + F[6]*F[7];
  b[2] = F[0]*F[2] + F[3]*F[5] + F[6]*F[8];
  b[3] = b[1];  
  b[4] = F[1]*F[1] + F[4]*F[4] + F[7]*F[7];
  b[5] = F[1]*F[2] + F[4]*F[5] + F[7]*F[8];
  b[6] = b[2];  
  b[7] = b[5];  
  b[8] = F[2]*F[2] + F[5]*F[5] + F[8]*F[8];
};
 
inline void calc_right_cauchy_green(double * F, double * C)
{
  C[0] = F[0]*F[0] + F[1]*F[1] + F[2]*F[2];
  C[1] = F[0]*F[3] + F[1]*F[4] + F[2]*F[5];
  C[2] = F[0]*F[6] + F[1]*F[7] + F[2]*F[8];
  C[3] = C[1];  
  C[4] = F[3]*F[3] + F[4]*F[4] + F[5]*F[5];
  C[5] = F[3]*F[6] + F[4]*F[7] + F[5]*F[8];
  C[6] = C[2];  
  C[7] = C[5];  
  C[8] = F[6]*F[6] + F[7]*F[7] + F[8]*F[8];
};
 
inline void calc_green_strain(double * F, double * E_Green) // Right now still with 9 entries, just for testing purposes
{
  E_Green[0] = 0.5 * (F[0]*F[0] + F[1]*F[1] + F[2]*F[2] - 1);
  E_Green[1] = 0.5 * (F[0]*F[3] + F[1]*F[4] + F[2]*F[5]);
  E_Green[2] = 0.5 * (F[0]*F[6] + F[1]*F[7] + F[2]*F[8]);
  E_Green[3] = E_Green[1];  
  E_Green[4] = 0.5 * (F[3]*F[3] + F[4]*F[4] + F[5]*F[5] - 1);
  E_Green[5] = 0.5 * (F[3]*F[6] + F[4]*F[7] + F[5]*F[8]);
  E_Green[6] = E_Green[2];  
  E_Green[7] = E_Green[5];  
  E_Green[8] = 0.5 * (F[6]*F[6] + F[7]*F[7] + F[8]*F[8] - 1);
};
 
inline void calc_green_strain_voigt(double * F, double * E_Green) // 6 entries: Voigt notation E = {E_11, E_22, E_33, 2*E12, 2*E13, 2*E23};
{
  E_Green[0] = 0.5 * (F[0]*F[0] + F[1]*F[1] + F[2]*F[2] - 1);
  E_Green[1] = 0.5 * (F[3]*F[3] + F[4]*F[4] + F[5]*F[5] - 1);
  E_Green[2] = 0.5 * (F[6]*F[6] + F[7]*F[7] + F[8]*F[8] - 1);
  E_Green[3] = F[0]*F[3] + F[1]*F[4] + F[2]*F[5];
  E_Green[4] = F[0]*F[6] + F[1]*F[7] + F[2]*F[8];
  E_Green[5] = F[3]*F[6] + F[4]*F[7] + F[5]*F[8];
};
 
inline Basic_constitutive_model * constitutive_model_setup(mechMat_t type, double * parameter_list)
{
  Basic_constitutive_model * constitutive_model;
  switch (type) {
    case NEOHOOKE:
      constitutive_model = new Neo_Hooke();
      constitutive_model->set_parameters(parameter_list);
      break;
    case GUCCIONE:
      constitutive_model = new Guccione();
      constitutive_model->set_parameters(parameter_list);
      break;
    default:
      break;
  }
  return constitutive_model;
}
}  // namespace opencarp
 
#endif