Structural_dynamic_mesh_model.h

/****************************************************************************
* Copyright (c) 2025, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#ifndef Structural_dynamic_mesh_model_included
#define Structural_dynamic_mesh_model_included
 
#include <Interprete_geometrique_base.h>
#include <stdexcept>
#include <vector>
#include <map>
 
#ifdef __NVCOMPILER
#pragma diag_suppress 68
#endif
#include <MGIS/Behaviour/Behaviour.hxx>
#include <MGIS/Behaviour/BehaviourDataView.hxx>
#include <MGIS/Behaviour/Integrate.hxx>
#ifdef __NVCOMPILER
#pragma diag_warning 68
#endif
 
// this enum defines the stress measure that the behavior is expected to return when calling 'integrate'
//   cauchy     -> cauchy stress tensor
//   pk1        -> first Piola-Kirchhoff tensor (aka Boussinesq tensor).
//                 This one is non-symmetric. It is such that P=pk1:dFdt, and Div(pk1^t)=f
//   pk2        -> second Piola-Kirchhoff tensor
//   conjugate  -> stress tensor conjugate with the input strain tensor, i.e. such that P=s:dedt
//                 where 'P' is the mechanical power, 's' is the stress tensor and 'dedt' the rate of the input strain tensor
enum class StressMeasureKind { cauchy, pk1, pk2, conjugate };
 
// this enum defines the tangent operator that the behavior is expected to return when calling 'integrate'
//   none      -> the tangent operator is not computed
//   dCauchydF -> derivative of the Cauchy stress tensor with respect to the gradient of the material deformation gradient
//   dPk2dE    -> derivative of the 2nd Piola-Kirchhoff stress tensor with respect to the Green-Lagrange strain tensor
//   dPk1dF    -> derivative of the 1st Piola-Kirchhoff stress tensor with respect to the gradient of the material deformation gradient
enum class TangentOperatorKind { none, dCauchydF, dPk2dE, dPk1dF };
 
class Domaine_ALE;
 
 
class Structural_dynamic_mesh_model : public Interprete_geometrique_base
{
  Declare_instanciable_sans_constructeur(Structural_dynamic_mesh_model);
 
 
public:
  Structural_dynamic_mesh_model();
  Entree& interpreter_(Entree&) override;
 
  // Mfront behaviour
  // ----------------
 
  inline void setMfrontLibraryPath(const std::string& l) ;
  inline void setMfrontModelName(const std::string& f) ;
  inline void setMfrontHypothesis(const std::string& h) ;
 
  // End Mfront behaviour
  // --------------------
 
  inline void setDensity(const double& rho) ;
  inline void setInertialDamping(const double& D) ;
  inline void setDtSafetyCoefficient(const double& c) ;
  inline void setGridDtMin(const double& dt) ;
  inline void setConfigurationResetDt(const double& dt) ;
  inline void setMaxAddedMassRatio(const double& r) ;
  inline void addMaterialProperty(const std::string name, const std::vector<double> val) ;
  inline int getNbNodesPerElem() ;
  inline double getDensity() ;
  inline double getDampingCoefficient() ;
  inline double getGridDtMin() ;
  inline void setMassElem(const double& elmass) ;
  inline const DoubleVect& getMeshPbPressure() const;
  inline const DoubleVect& getMeshPbVonMises() const;
  inline const DoubleTab& getMeshPbForceFace() const;
  inline const DoubleTab& getMeshDisplacement() const {return u; }
  inline const DoubleTab& getMeshVelocity() const {return v; }
  inline const DoubleTab& getMeshAcceleration() const {return a; }
  inline const DoubleTab& getMeshPosition() const {return x; }
  inline const DoubleTab& getMeshReferenceConfiguration() const {return B0_; }
  inline const DoubleTab& getMeshTransformationGradient() const {return Ft_; }
  inline const DoubleTab& getMeshStress() const {return Stress_; }
  inline const int& getMeshReferenceConfigurationNbComp() const {return dimB0_;}
  inline const int& getMeshTransformationGradientNbComp() const {return nSymSize_;}
  inline const int& getMeshStressNbComp() const {return symSize_;}
 
  inline void applyDtCoefficient() ;
 
  void initMfrontBehaviour() ;
  void initDynamicMeshProblem(const double temps, const int nsom, const int nelem, const int nface, const MD_Vector& md, const MD_Vector& mde, const MD_Vector& mdf) ;
 
  void setLocalFields(const int elnodes[4], const int elem) ;
  void computeInternalForces(double& volume, double& xlong, double& cSound, double& Pressure, double& VonMises) ;
  void setGlobalFields(const int elnodes[4], const double Pressure, const double VonMises) ;
  double computeCriticalDt(const double volume, const double xlong, const double cSound, const double dtMin, double& scaleMass) ;
  void computeForceFaces(const int nb_faces, const int nb_som_face, const IntTab& face_sommets) ;
  void checkElemOrientation(int elnodes[4], const int elem) ;
  void resumptionMesh(double, DoubleTab&, DoubleTab&, DoubleTab&, DoubleTab&, DoubleTab&, DoubleTab&, DoubleTab&);
  inline void saveBeginning_of_stepVariables() ;
  inline void resetVariablestoBeginning_of_step() ;
 
  void solveDynamicMeshProblem(const double, const DoubleTab&, const IntVect&,
                               DoubleTab&, const int , const int , const int ,
                               const IntTab& , const int , const int , const IntTab& , const Domaine_ALE&) ;
 
  // Global vectors and arrays for dynamic time integration
  // ------------------------------------------------------
 
  DoubleTab x ;
  DoubleTab u ;
  DoubleTab v ;
  DoubleTab vp ;
  DoubleTab a ;
  DoubleTab ff ;
 
  DoubleVect mass ;
  DoubleVect nodalScaleMass ;
  double totalMass = 0.;
 
  // Public variables
  // ----------------
 
  int gridNStep = 0.;
  double gridTime = 0.;
  double gridTime_n=0. ;
  double gridDt = 0.;
  bool isMassBuilt = false;
 
  double configurationResetDt = 0. ;
  double gridResetTime = 0.;
 
  double maxAddedMassRatio =0.;
  bool doConfigurationReset=false ;
 
  int resumption =0; //1 if resumption of calculation else 0
 
  // ICoCo + FSI: acceleration solution within implicit iterations
  // ----------------------------------------------------------------
 
  int iCoCoImplicitIteration = -1 ; // Default, no ICoCo implicit coupling
  bool acceleratedSolutionEnabled = true;
  int acceleratedSolution = 0 ;
  DoubleTab xLast ;
  DoubleTab vpLast ;
  double dtLast ;
  bool skipStressComputation = false ;
 
protected:
 
  // Mfront behaviour
  // ----------------
 
  mgis::behaviour::Behaviour mgisBehaviour_;
  mgis::behaviour::BehaviourDataView mgisBehaviourData_;
  std::vector<double> s0MaterialProperties_;
  std::vector<double> s1MaterialProperties_;
 
  static constexpr int stiffnessMatrixMinSize_ = 1 + mgis::behaviour::Behaviour::nopts;
 
  void load_behaviour_(StressMeasureKind smk, TangentOperatorKind tok);
 
 
  std::string l_;    // behavior library path
  std::string f_;    // name of the behavior in the library
  std::string h_;    // modelling hypothesis (tridimentional, plane strain, plane stress, ...)
  /*
  Valid values for h_ are:
   - `AxisymmetricalGeneralisedPlaneStrain`
   - `AxisymmetricalGeneralisedPlaneStress`
   - `Axisymmetrical`
   - `PlaneStress`
   - `PlaneStrain`
   - `GeneralisedPlaneStrain`
   - `Tridimensional`
  */
 
  std::map<std::string,std::vector<double>> matp_ ; // MFront material properties
 
  double rdt_ = 1.;    // this is used by mgis
 
  // numerotation, mfront convention:
  // for symetric tensors:
  //  - 2D: [t00, t11, t22, sqrt(2)*t01]
  //  - 3D: [t00, t11, t22, sqrt(2)*t01, sqrt(2)*t02, sqrt(2)*t12]
  // for non-symetric tensors:
  //  - 2D: [t00, t11, t22, t01, t10]
  //  - 3D: [t00, t11, t22, t01, t10, t02, t20, t12, t21]
  void integrate_behaviour_(double *const stress,                     // stress tensor
                            double *const ivar,                       // state variables
                            const double *const evar,                 // external variables
                            double *const stiff,                      // stiffness matrix
                            const double *const gradientOrStrain0,    // deformation gradient or strain at the beginning of the time step
                            const double *const gradientOrStrain1,    // deformation gradient or strain at the end of the time step
                            double *const vsound                      // speed of sound variables
                           );
 
  // End Mfront behaviour
  // --------------------
 
  // Parameters and global variables
  // -------------------------------
 
  double inertialDamping_=0. ;
  double density_ =0. ;
  double dtSafetyCoefficient_ = 0.5 ;
  double gridDtMin_ = 0 ;
 
  int nSymSize_ =0; // length of non-symetric tensor in Voigt notation
  int symSize_ =0; // length of symetric tensor in Voigt notation
  int dimB0_ =0; // lenght of reference configuration array
  DoubleTab Eta_ ; // Shape function derivatives
 
  // Local variables for internal forces computation
  // -----------------------------------------------
 
  int iel_ =0; // id of current element
  int nbn_ =0; // number of vertex nodes
  DoubleTab xl_ ; // local current coordinates
  DoubleTab ul_ ; // local nodal displacements
  DoubleTab fl_ ; // local nodal forces
 
  DoubleTab B0l_ ; // local B matrix on initial configuration
 
  IntVect invertNum_ ; // need for inverting elem numerotation or not
 
  // Global vectors and arrays for dynamic time integration
  // ------------------------------------------------------
 
  DoubleTab B0_ ;
  DoubleTab Ft_ ;
  DoubleTab Stress_ ;
 
  DoubleVect massElem_ ;
 
  // Fields for post-processing
  DoubleVect meshPbPressure_ ;
  DoubleVect meshPbVonMises_ ;
  DoubleTab  meshPbForceFace_ ;
 
  // ICoCo + FSI: save variables for implicit + subcycling iterations
  // ----------------------------------------------------------------
 
  double gridDt_n ;
  int gridNStep_n ;
  DoubleTab x_n ;
  DoubleTab u_n ;
  DoubleTab v_n ;
  DoubleTab a_n ;
  DoubleTab B0_n ;
  DoubleTab Ft_n ;
  DoubleTab Stress_n ;
 
  // ICoCo + FSI: acceleration solution within implicit iterations
  // ----------------------------------------------------------------
 
  DoubleVect cSoundElem_ ;
 
  // Functions
  // ---------
 
  void triangleSurfLength_(double& aire, double& xlong, const double Det);
  void tetrahedronVolLength_(double& vol, double& xlong, const double Det);
  void setB0_(const DoubleTab& B) ;
 
private:
 
  // Mfront behaviour
  // ----------------
 
  int nbIvars_=0;
  int nbEvars_ =0;
  int sizeEvars_ =0;
  int matpSize_ =0;
  int KSize_ =0;
  mgis::behaviour::Hypothesis hypothesis_;
 
  bool loaded_ = false;
 
  DoubleTab mfrontEvars_ ;
 
  // End Mfront behaviour
  // --------------------
 
};
 
inline void Structural_dynamic_mesh_model::setMfrontLibraryPath(const std::string& l)
{
  l_ = l ;
}
 
inline void Structural_dynamic_mesh_model::setMfrontModelName(const std::string& f)
{
  f_ = f ;
}
 
inline void Structural_dynamic_mesh_model::setMfrontHypothesis(const std::string& h)
{
  h_ = h ;
}
 
inline void Structural_dynamic_mesh_model::setDensity(const double& rho)
{
  density_ = rho ;
}
 
inline void Structural_dynamic_mesh_model::setInertialDamping(const double& D)
{
  inertialDamping_ = D ;
}
 
inline void Structural_dynamic_mesh_model::setDtSafetyCoefficient(const double& c)
{
  dtSafetyCoefficient_ = c ;
}
 
inline void Structural_dynamic_mesh_model::setGridDtMin(const double& dt)
{
  gridDtMin_ = dt ;
}
 
inline void Structural_dynamic_mesh_model::setConfigurationResetDt(const double& dt)
{
  configurationResetDt = dt ;
  gridResetTime = configurationResetDt ; // first reset time
}
 
inline void Structural_dynamic_mesh_model::setMaxAddedMassRatio(const double& r)
{
  maxAddedMassRatio = r ;
}
 
inline void Structural_dynamic_mesh_model::addMaterialProperty(const std::string name, const std::vector<double> val)
{
  matp_.insert(std::pair<std::string, std::vector<double>>(name, val)) ;
}
 
inline int Structural_dynamic_mesh_model::getNbNodesPerElem()
{
  return nbn_ ;
}
 
inline double Structural_dynamic_mesh_model::getDensity()
{
  return density_ ;
}
 
inline double Structural_dynamic_mesh_model::getDampingCoefficient()
{
  return inertialDamping_ ;
}
 
inline void Structural_dynamic_mesh_model::setMassElem(const double& elmass)
{
  massElem_[iel_] = elmass ;
}
 
inline const DoubleVect& Structural_dynamic_mesh_model::getMeshPbPressure() const
{
  return meshPbPressure_ ;
}
 
inline const DoubleVect& Structural_dynamic_mesh_model::getMeshPbVonMises() const
{
  return meshPbVonMises_ ;
}
 
inline const DoubleTab& Structural_dynamic_mesh_model::getMeshPbForceFace() const
{
  return meshPbForceFace_ ;
}
 
inline void Structural_dynamic_mesh_model::applyDtCoefficient()
{
  gridDt *= dtSafetyCoefficient_ ;
}
 
inline double Structural_dynamic_mesh_model::getGridDtMin()
{
  return gridDtMin_ ;
}
 
inline void Structural_dynamic_mesh_model::saveBeginning_of_stepVariables()
{
  Cerr << "implicit iterations with structural grid problem: save beginning-of-step variables" << finl ;
  gridTime_n = gridTime ;
  gridDt_n = gridDt ;
  gridNStep_n = gridNStep ;
  x_n = x ;
  u_n = u ;
  v_n = v ;
  a_n = a ;
  B0_n = B0_ ;
  Ft_n = Ft_ ;
  Stress_n = Stress_ ;
}
 
inline void Structural_dynamic_mesh_model::resetVariablestoBeginning_of_step()
{
  Cerr << "implicit iterations with structural grid problem: reset variables to beginning-of-step" << finl ;
  gridTime = gridTime_n ;
  gridDt = gridDt_n ;
  gridNStep = gridNStep_n ;
  x = x_n ;
  u = u_n ;
  v = v_n ;
  a = a_n ;
  B0_ = B0_n ;
  Ft_ = Ft_n ;
  Stress_ = Stress_n ;
}
 
#endif