Structural_dynamic_mesh_model.cpp

/****************************************************************************
* 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.
*
*****************************************************************************/
 
#include <Structural_dynamic_mesh_model.h>
 
#include <Domaine_ALE.h>
#include <stdexcept>
 
Implemente_instanciable_sans_constructeur(Structural_dynamic_mesh_model,"Structural_dynamic_mesh_model",Interprete_geometrique_base);
 
// Syntaxe:
//Structural_dynamic_mesh_model NOMDOMAINE
//{
//        Mfront_library path_to_libBehaviour.so
//        Mfront_model_name Ogden|SaintVenantKirchhoffElasticity
//        Mfront_material_property
//        [ Grid_dt_min ]
//        [ YoungModulus ]
//        [ Density ]
//        [ Inertial_Damping ]
//}
 
// XD Structural_dynamic_mesh_model interprete Structural_dynamic_mesh_model NO_BRACE Fictitious structural model for
// XD_CONT mesh motion. Link with MGIS library
// XD attr dom ref_domaine dom REQ domain name
// XD attr bloc bloc_lecture_Structural_dynamic_mesh_model bloc REQ not_set
 
// XD bloc_lecture_Structural_dynamic_mesh_model objet_lecture nul NO_BRACE bloc
// XD attr aco chaine(into=["{"]) aco REQ Opening curly bracket.
// XD attr Mfront_library chaine(into=["Mfront_library"]) Mfront_library OPT Keyword to specify the
// XD_CONT path_to_libBehaviour.so If the user does not define the MFront library path, we use the default one instead:
// XD_CONT path = \$project_directory/share/MFront_material_libraries/<MFront_library_keyword>/src/libBehaviour.so
// XD attr Mfront_model_name chaine(into=["Mfront_model_name"]) Mfront_model OPT keyword to specify the Mfront model.
// XD_CONT Choice between Ogden and SaintVenantKirchhoffElasticity.
// XD attr Mfront_material_property chaine(into=["Mfront_material_property"]) Mfront_material_property OPT keyword to
// XD_CONT specify the material property. Eg. Ogden_alpha_, Ogden_mu_, Ogden_K
// XD attr YoungModulus floattant young OPT Young Module
// XD attr Density floattant rho OPT fictitious structural density
// XD attr Inertial_Damping floattant inertial_Damping OPT fictitious structural inertial damping
// XD attr Grid_dt_min floattant grid_dt_min OPT fictitious structural time step
// XD attr acof chaine(into=["}"]) acof REQ Closing curly bracket.
 
Structural_dynamic_mesh_model::Structural_dynamic_mesh_model()
{
 
  l_ = "unset string" ;
  f_ = "unset string" ;
  h_ = "unset string" ;
 
  nSymSize_ = (dimension == 2) ? 5 : 9;
  symSize_  = (dimension == 2) ? 4 : 6;
  nbn_ = (dimension == 2) ? 3 : 4;
  dimB0_ =nbn_*dimension;
 
  x.reset();
  u.reset();
  v.reset();
  vp.reset();
  a.reset();
  ff.reset();
  mass.reset() ;
  nodalScaleMass.reset() ;
  Eta_.reset() ;
  xl_.reset() ;
  ul_.reset() ;
  fl_.reset() ;
  B0l_ .reset() ;
  invertNum_.reset() ;
  B0_.reset() ;
  Ft_.reset() ;
  Stress_.reset() ;
  massElem_ .reset();
  meshPbPressure_.reset() ;
  meshPbVonMises_.reset() ;
  meshPbForceFace_.reset() ;
  mfrontEvars_.reset() ;
}
 
Sortie& Structural_dynamic_mesh_model::printOn(Sortie& os) const
{
  return os;
}
 
 
Entree& Structural_dynamic_mesh_model::readOn(Entree& is)
{
  return is;
}
 
Entree& Structural_dynamic_mesh_model::interpreter_(Entree& is)
{
 
  associer_domaine(is);
  Domaine_ALE& dom=ref_cast(Domaine_ALE, domaine());
  dom.reading_structural_dynamic_mesh_model(is);
  return is;
}
 
void Structural_dynamic_mesh_model::initMfrontBehaviour()
{
 
  if ( l_ == "unset string")
    {
      std::string path;
      if(f_ != "unset string")
        {
          // If the user does not define the MFront library path,
          // use the default one instead:
          // path = $project_directory/share/MFront_material_libraries/<MFront_library_keyword>/src/libBehaviour.so
          path = std::string(getenv("TrioCFD_project_directory")) + "/share/Mfront_material_libraries/" + f_+ "/src/libBehaviour.so";
          l_=path;
        }
      else
        {
          Cerr << "MFront model name is undefined (keyword: Mfront_model_name)." << finl;
          Cerr << "Cannot set the default MFront library path: "
               << "$project_directory/share/MFront_material_libraries/<Mfront_library_keyword>/src/libBehaviour.so" << finl;
          Cerr << "Please define either the MFront library path (keyword: Mfront_library) "
               << "or the MFront model name (keyword: Mfront_model_name)." << finl;
          Process::exit() ;
        }
    }
 
  if ( f_ == "unset string")
    {
      Cerr << "Error: Mfront model name undefined (Mfront_model_name keyword)" << finl ;
      Process::exit() ;
    }
 
  if (h_ == "unset string")
    {
      switch (dimension)
        {
        case (2) :
          Cerr << "Use default Mfront hypothesis for 2D (PlaneStrain)" << finl ;
          h_ = "PlaneStrain" ;
          break ;
        case (3) :
          Cerr << "Use default Mfront hypothesis for 3D (Tridimensional)" << finl ;
          h_ = "Tridimensional" ;
          break ;
        }
 
    }
  hypothesis_= mgis::behaviour::fromString(h_) ;
 
  // dt/rdt
  mgisBehaviourData_.dt  = 0.;
  mgisBehaviourData_.rdt = &rdt_;
 
  load_behaviour_(StressMeasureKind::cauchy, TangentOperatorKind::none) ;
  if ( nbIvars_ > 0)
    {
      Cerr << "Mfront behaviour: internal state variables:" << finl ;
      for (auto isvs : mgisBehaviour_.isvs)
        {
          size_t vOffset = mgis::behaviour::getVariableOffset(mgisBehaviour_.isvs, isvs.name, hypothesis_);
          size_t vSize   = mgis::behaviour::getVariableSize(isvs, hypothesis_);
          Cerr << "    - " << isvs.name << ": " << vOffset << " -> " << vOffset + vSize << finl ;
        }
 
      Cerr << "Wrong Mfront behaviour for structural dynamic mesh motion" << finl ;
      Cerr << "Only pure elastic materials are allowed" << finl ;
      Process::exit();
    }
  if ( nbEvars_ > 0)
    {
      Cerr << "Mfront behaviour: external state variables:" << finl ;
      sizeEvars_ = 0 ;
      for (auto esvs : mgisBehaviour_.esvs)
        {
          size_t vOffset = mgis::behaviour::getVariableOffset(mgisBehaviour_.esvs, esvs.name, hypothesis_);
          size_t vSize   = mgis::behaviour::getVariableSize(esvs, hypothesis_);
          Cerr << "    - " << esvs.name << ": " << vOffset << " -> " << vOffset + vSize << finl ;
          sizeEvars_ += static_cast<int> (vSize) ;
        }
 
      Cerr << "Warning: external variables ignored for mesh motion" << finl ;
    }
 
  // Check and fill material properties for Mfront behaviour
  matpSize_                                 = static_cast<int> (mgis::behaviour::getArraySize(mgisBehaviour_.mps, hypothesis_));
  s0MaterialProperties_.resize(matpSize_);
  s1MaterialProperties_.resize(matpSize_);
  mgisBehaviourData_.s0.material_properties = s0MaterialProperties_.data();
  mgisBehaviourData_.s1.material_properties = s1MaterialProperties_.data();
  for (auto mps : mgisBehaviour_.mps)
    {
      size_t vOffset = mgis::behaviour::getVariableOffset(mgisBehaviour_.mps, mps.name, hypothesis_);
      size_t vSize   = mgis::behaviour::getVariableSize(mps, hypothesis_);
      auto it        = matp_.find(mps.name);
      if (it != matp_.end())
        {
          if (vSize != (it->second).size())
            {
              Cerr << "Error with input for Mfront behaviour " << f_ << finl ;
              Cerr << "The number of scalar values must be " << vSize << " for property " << mps.name << finl ;
              Process::exit();
            }
          for (size_t i = 0; i < vSize; ++i)
            {
              s0MaterialProperties_[vOffset + i] = it->second[i];
              s1MaterialProperties_[vOffset + i] = it->second[i];
            }
        }
      else // missing property: set its value(s) to zero with a warning
        {
          Cerr << "Warning with input for Mfront behaviour " << f_ << finl ;
          Cerr << "Missing data for for property " << mps.name << ", default value(s) set to zero" << finl ;
          for (size_t i = 0; i < vSize; ++i)
            {
              s0MaterialProperties_[vOffset + i] = 0;
              s1MaterialProperties_[vOffset + i] = 0;
            }
        }
    }
  mgisBehaviourData_.s0.mass_density = &density_;
  mgisBehaviourData_.s1.mass_density = &density_;
 
}
 
void Structural_dynamic_mesh_model::initDynamicMeshProblem(const double temps, const int nsom, const int nelem, const int nface, const MD_Vector& mds, const MD_Vector& mde, const MD_Vector& mdf)
{
  switch (dimension)
    {
    case (2) :
      nbn_ = 3 ;
      nSymSize_ = 5 ;
      symSize_ = 4 ;
 
      Eta_.resize(dimension,nbn_) ;
      Eta_(0,0) = -1. ;
      Eta_(0,1) =  1. ;
      Eta_(0,2) =  0. ;
      Eta_(1,0) = -1. ;
      Eta_(1,1) =  0. ;
      Eta_(1,2) =  1. ;
 
      break ;
    case (3) :
      nbn_ = 4 ;
      nSymSize_ = 9 ;
      symSize_ = 6 ;
 
      Eta_.resize(dimension,nbn_) ;
      Eta_(0,0) = -1. ;
      Eta_(0,1) =  1. ;
      Eta_(0,2) =  0. ;
      Eta_(0,3) =  0. ;
      Eta_(1,0) = -1. ;
      Eta_(1,1) =  0. ;
      Eta_(1,2) =  1. ;
      Eta_(1,3) =  0. ;
      Eta_(2,0) = -1. ;
      Eta_(2,1) =  0. ;
      Eta_(2,2) =  0. ;
      Eta_(2,3) =  1. ;
 
      break ;
    }
 
  if(!resumption)
    {
      u.resize(nsom,dimension) ;
      MD_Vector_tools::creer_tableau_distribue(mds, u);
      v.resize(nsom,dimension) ;
      MD_Vector_tools::creer_tableau_distribue(mds, v);
 
      a.resize(nsom,dimension) ;
      MD_Vector_tools::creer_tableau_distribue(mds, a);
 
      x.resize(nsom,dimension) ;
      MD_Vector_tools::creer_tableau_distribue(mds, x);
 
    }
 
  vp.resize(nsom,dimension) ;
  ff.resize(nsom,dimension) ;
 
  mass.resize(nsom) ;
  if (getGridDtMin() > 0.) nodalScaleMass.resize(nsom) ;
 
  dimB0_=dimension*nbn_;
  if(!resumption)
    {
      B0_.resize(nelem,dimB0_) ;
      MD_Vector_tools::creer_tableau_distribue(mde, B0_);
      Ft_.resize(nelem, nSymSize_) ;
      MD_Vector_tools::creer_tableau_distribue(mde, Ft_);
      Stress_.resize(nelem, symSize_) ;
      MD_Vector_tools::creer_tableau_distribue(mde, Stress_);
 
    }
 
  invertNum_.resize(nelem) ;
 
  massElem_.resize(nelem) ;
 
  cSoundElem_.resize(nelem) ;
 
  meshPbPressure_.resize(nelem) ;
  meshPbVonMises_.resize(nelem) ;
  meshPbForceFace_.resize(nface,dimension) ;
 
  mfrontEvars_.resize(nelem, sizeEvars_) ;
 
  // Only ff and mass need to be built as distributed arrays
  MD_Vector_tools::creer_tableau_distribue(mds, ff) ;
  MD_Vector_tools::creer_tableau_distribue(mds, mass) ;
  if (getGridDtMin() > 0.) MD_Vector_tools::creer_tableau_distribue(mds, nodalScaleMass) ;
 
  if(!resumption)
    {
      u=0 ;
      v=0 ;
      a=0 ;
    }
  vp=0 ;
  ff=0 ;
  mass=0. ;
  if (getGridDtMin() > 0.) nodalScaleMass=0. ;
 
  if(!resumption)
    {
      B0_ = 0. ;
      Stress_ = 0. ;
    }
  invertNum_ = 0 ; // 0 = no numerotation inversion by default
  massElem_ = 0. ;
  mfrontEvars_ = 0. ;
 
  // Arrays for postprocessing
  meshPbPressure_ = 0. ;
  meshPbVonMises_ = 0. ;
  meshPbForceFace_ = 0. ;
  MD_Vector_tools::creer_tableau_distribue(mde, meshPbPressure_) ;
  MD_Vector_tools::creer_tableau_distribue(mde, meshPbVonMises_) ;
  MD_Vector_tools::creer_tableau_distribue(mdf, meshPbForceFace_) ;
 
  // Initial transformation gradient is identity
  if(!resumption)
    {
      for(int i=0; i<nelem; i++)
        {
          switch (dimension)
            {
            case (2) :
 
              Ft_(i,0) = 1. ;
              Ft_(i,1) = 1. ;
              Ft_(i,2) = 1. ;
              Ft_(i,3) = 0. ;
              Ft_(i,4) = 0. ;
 
              break ;
 
            case (3) :
 
              Ft_(i,0) = 1. ;
              Ft_(i,1) = 1. ;
              Ft_(i,2) = 1. ;
              Ft_(i,3) = 0. ;
              Ft_(i,4) = 0. ;
              Ft_(i,5) = 0. ;
              Ft_(i,6) = 0. ;
              Ft_(i,7) = 0. ;
              Ft_(i,8) = 0. ;
 
              break ;
            }
        }
    }
 
  xl_.resize(nbn_,dimension) ;
  ul_.resize(nbn_,dimension) ;
  fl_.resize(nbn_,dimension) ;
 
  B0l_.resize(dimension,nbn_) ;
 
  xLast = x ;
  vpLast =vp ;
 
}
 
void Structural_dynamic_mesh_model::computeInternalForces(double& Vol, double& Xlong, double& cSound, double& Pressure, double& VonMises)
{
 
  // Integration with 1 single gauss point on triangular (dimension=2) or tetrahedral (dimension=3) finite element
 
  //  Compute jacobian and B matrices
 
  DoubleTab Jac(dimension, dimension) ;
 
  Jac=0 ;
  for (int k=0; k < nbn_; k++)
    {
      for (int j=0; j < dimension; j++)
        {
          for (int i=0; i < dimension ; i++)
            {
              Jac(j,i) += Eta_(j,k) * xl_(k,i) ;
            }
        }
    }
 
  double Det ;
  DoubleTab InvJac(dimension, dimension) ;
  double aux ;
  switch (dimension)
    {
    case (2) :
      Det = Jac(0,0) * Jac(1,1) - Jac(0,1) * Jac(1,0) ;
      aux = 1. / Det ;
      InvJac(0,0) =  Jac(1,1) * aux ;
      InvJac(0,1) = -Jac(0,1) * aux ;
      InvJac(1,0) = -Jac(1,0) * aux ;
      InvJac(1,1) =  Jac(0,0) * aux ;
 
      triangleSurfLength_(Vol, Xlong, Det) ;
 
      break ;
    case (3) :
      Det= Jac(0,0) * (Jac(1,1) * Jac(2,2) - Jac(1,2) * Jac(2,1))
           - Jac(0,1) * (Jac(1,0) * Jac(2,2) - Jac(1,2) * Jac(2,0))
           + Jac(0,2) * (Jac(1,0) * Jac(2,1) - Jac(1,1) * Jac(2,0)) ;
      aux = 1. / Det ;
      InvJac(0,0) =  (Jac(1,1) * Jac(2,2) - Jac(2,1) * Jac(1,2)) * aux ;
      InvJac(0,1) = -(Jac(0,1) * Jac(2,2) - Jac(2,1) * Jac(0,2)) * aux ;
      InvJac(0,2) =  (Jac(0,1) * Jac(1,2) - Jac(1,1) * Jac(0,2)) * aux ;
      InvJac(1,0) = -(Jac(1,0) * Jac(2,2) - Jac(2,0) * Jac(1,2)) * aux ;
      InvJac(1,1) =  (Jac(0,0) * Jac(2,2) - Jac(2,0) * Jac(0,2)) * aux ;
      InvJac(1,2) = -(Jac(0,0) * Jac(1,2) - Jac(1,0) * Jac(0,2)) * aux ;
      InvJac(2,0) =  (Jac(1,0) * Jac(2,1) - Jac(2,0) * Jac(1,1)) * aux ;
      InvJac(2,1) = -(Jac(0,0) * Jac(2,1) - Jac(2,0) * Jac(0,1)) * aux ;
      InvJac(2,2) =  (Jac(0,0) * Jac(1,1) - Jac(1,0) * Jac(0,1)) * aux ;
 
      tetrahedronVolLength_(Vol, Xlong, Det) ;
 
      break ;
    }
 
  DoubleTab B(dimension, nbn_) ;
  B=0 ;
  for (int k=0; k < dimension; k++)
    {
      for (int j=0; j < dimension; j++)
        {
          for (int i=0; i < nbn_ ; i++)
            {
              B(j,i) += InvJac(j,k) * Eta_(k,i) ;
            }
        }
    }
 
  if (gridNStep == 0)
    {
      setB0_(B) ; // save B matrix on initial configuration
      B0l_ = B ;
    }
 
  if (doConfigurationReset)  // take current configuration as the reference configuration
    // for further mesh deformation
    {
      setB0_(B) ;
      B0l_ = B ;
 
      for (int i=0; i<symSize_; i++) Stress_(iel_,i) = 0. ;
      switch (dimension)
        {
        case (2) :
 
          Ft_(iel_,0) = 1. ;
          Ft_(iel_,1) = 1. ;
          Ft_(iel_,2) = 1. ;
          Ft_(iel_,3) = 0. ;
          Ft_(iel_,4) = 0. ;
 
          break ;
 
        case (3) :
 
          Ft_(iel_,0) = 1. ;
          Ft_(iel_,1) = 1. ;
          Ft_(iel_,2) = 1. ;
          Ft_(iel_,3) = 0. ;
          Ft_(iel_,4) = 0. ;
          Ft_(iel_,5) = 0. ;
          Ft_(iel_,6) = 0. ;
          Ft_(iel_,7) = 0. ;
          Ft_(iel_,8) = 0. ;
 
          break ;
        }
    }
 
#define _xikbjk(i,j) xl_(k,i)*B0l_(j,k)
  std::vector <double> stress(symSize_) ;
 
  if (!skipStressComputation)
    {
 
      //  Next transformation gradient
 
      std::vector <double> ftpdt(nSymSize_) ;
      for (int i=0; i<nSymSize_; i++) ftpdt[i] = 0 ;
 
      for (int k=0; k<nbn_; k++)
        {
          ftpdt[0] += _xikbjk(0,0); // F11=x1k*dNk/dx01
          ftpdt[1] += _xikbjk(1,1); // F22=x2k*dNk/dx02
          ftpdt[3] += _xikbjk(0,1); // F12=x1k*dNk/dx02
          ftpdt[4] += _xikbjk(1,0); // F21=x2k*dNk/dx01
          if (dimension == 2)
            {
              ftpdt[2] = 1. ;
            }
          else if (dimension == 3)
            {
              ftpdt[2] += _xikbjk(2,2); // F33=x3k*dNk/dx03
              ftpdt[5] += _xikbjk(0,2); // F13=x1k*dNk/dx03
              ftpdt[6] += _xikbjk(2,0); // F31=x3k*dNk/dx01
              ftpdt[7] += _xikbjk(1,2); // F23=x2k*dNk/dx03
              ftpdt[8] += _xikbjk(2,1); // F32=x3k*dNk/dx02
            }
        }
 
      //  Integrate material behavior
 
      std::vector <double> ft(nSymSize_) ;
      for (int i=0; i<nSymSize_; i++) ft[i] = Ft_(iel_,i) ;
 
      for (int i=0; i<symSize_; i++) stress[i] = Stress_(iel_,i) ;
 
      std::vector <double> evars(sizeEvars_) ;
      for (int i=0; i<sizeEvars_; i++) evars[i] = mfrontEvars_(iel_,i) ;
 
      double voidInternalVars[1] ; // 0 size forbidden
      std::vector <double>  K(stiffnessMatrixMinSize_) ;
      K[0] = 100; // Stiffness matrix not computed, see mfront/mgis convention
      // Value >= 50: compute speed of sound
 
      double speedOfSound[1] ;
 
      integrate_behaviour_(&stress[0], voidInternalVars, &evars[0], &K[0], &ft[0], &ftpdt[0], speedOfSound) ;
 
      if (speedOfSound[0] <= 0.)
        {
          Cerr << "Error (mesh motion): speed of sound not computed in provided Mfront behaviour" << finl ;
          Process::exit();
        }
      else
        {
          cSound = speedOfSound[0];
          cSoundElem_(iel_) = cSound;
        }
 
      // Update stress and transformation gradient
      for (int i=0; i<symSize_; i++) Stress_(iel_,i) = stress[i] ;
      for (int i=0; i<nSymSize_; i++) Ft_(iel_,i) = ftpdt[i] ;
    }
  else
    {
      // Skipped stress computation: take stored values for stress and cSound
      for (int i=0; i<symSize_; i++) stress[i] = Stress_(iel_,i) ;
      cSound = cSoundElem_(iel_) ;
    }
 
  //  Compute local forces
  aux = 1./sqrt(2.) ;
  fl_ = 0 ;
  std::vector <double> cauchy(symSize_) ;
  switch (dimension)
    {
    case (2) :
      // Switch from Voigt notation (Mfront) to EPX notation
      cauchy[0] = stress[0] ;
      cauchy[1] = stress[1] ;
      cauchy[2] = stress[3] * aux ;
      cauchy[3] = stress[2] ;
 
      for (int k=0; k<nbn_; k++)
        {
          fl_(k,0) += Vol * (B(0,k) * cauchy[0] + B(1,k) * cauchy[2]) ;
          fl_(k,1) += Vol * (B(0,k) * cauchy[2] + B(1,k) * cauchy[1]) ;
        }
 
      Pressure = (cauchy[0] + cauchy[1] + cauchy[3]) / 3. ;
      // abs to suppress the risk of negative values due to truncation with very small values
      VonMises = sqrt(abs(cauchy[0] * (cauchy[0] - cauchy[1]) + cauchy[1] * (cauchy[1] - cauchy[3])
                          + cauchy[3] * (cauchy[3] - cauchy[0]) + 3. * cauchy[2] * cauchy[2])) ;
 
      break ;
    case (3) :
      // Switch from Voigt notation (Mfront) to EPX notation
      cauchy[0] = stress[0] ;
      cauchy[1] = stress[1] ;
      cauchy[2] = stress[2] ;
      cauchy[3] = stress[3] * aux ;
      cauchy[4] = stress[5] * aux ;
      cauchy[5] = stress[4] * aux ;
 
      for (int k=0; k<nbn_; k++)
        {
          fl_(k,0) += Vol * (B(0,k) * cauchy[0] + B(1,k) * cauchy[3] + B(2,k) * cauchy[5]) ;
          fl_(k,1) += Vol * (B(0,k) * cauchy[3] + B(1,k) * cauchy[1] + B(2,k) * cauchy[4]) ;
          fl_(k,2) += Vol * (B(0,k) * cauchy[5] + B(1,k) * cauchy[4] + B(2,k) * cauchy[2]) ;
        }
 
      Pressure = (cauchy[0] + cauchy[1] + cauchy[2]) / 3. ;
      // abs to suppress the risk of negative values due to truncation with very small values
      VonMises = sqrt(abs(cauchy[0] * (cauchy[0] - cauchy[1]) + cauchy[1] * (cauchy[1] - cauchy[2]) + cauchy[2] * (cauchy[2] - cauchy[0])
                          + 3. * (cauchy[3] * cauchy[3] + cauchy[4] * cauchy[4] + cauchy[5] * cauchy[5]))) ;
 
      break ;
    }
 
}
 
void Structural_dynamic_mesh_model::triangleSurfLength_(double& aire, double& xlong, const double Det)
{
 
  aire = Det / 2. ;
 
  double v12[2] ;
  double v13[2] ;
  double v23[2] ;
 
  for (int i=0; i<dimension; i++)
    {
      v12[i] = xl_(1,i) - xl_(0,i) ;
      v13[i] = xl_(2,i) - xl_(0,i) ;
      v23[i] = xl_(2,i) - xl_(1,i) ;
    }
 
  double w[3] ;
 
  w[0] = v12[0] * v12[0] + v12[1] * v12[1] ;
  w[1] = v13[0] * v13[0] + v13[1] * v13[1] ;
  w[2] = v23[0] * v23[0] + v23[1] * v23[1] ;
 
  double max_value = std::max({w[0], w[1], w[2]}) ;
  xlong = Det / sqrt(max_value) ;
 
}
 
void Structural_dynamic_mesh_model::tetrahedronVolLength_(double& vol, double& xlong, const double Det)
{
 
  vol = Det / 6. ;
 
  double v21[3] ;
  double v31[3] ;
  double v41[3] ;
  double v32[3] ;
  double v42[3] ;
 
  for (int i=0; i<dimension; i++)
    {
      v21[i] = xl_(0,i) - xl_(1,i) ;
      v31[i] = xl_(0,i) - xl_(2,i) ;
      v41[i] = xl_(0,i) - xl_(3,i) ;
      v32[i] = xl_(1,i) - xl_(2,i) ;
      v42[i] = xl_(1,i) - xl_(3,i) ;
    }
 
  double w[3] ;
  double s[4] ;
 
  // face 1 2 3 : cross product 21^31
  w[0] = v21[1]*v31[2] - v21[2]*v31[1] ;
  w[1] = v21[2]*v31[0] - v21[0]*v31[2] ;
  w[2] = v21[0]*v31[1] - v21[1]*v31[0] ;
  s[0] = w[0]*w[0] + w[1]*w[1] + w[2]*w[2] ;
  // face 1 2 4 : cross product 21^41
  w[0] = v21[1]*v41[2] - v21[2]*v41[1] ;
  w[1] = v21[2]*v41[0] - v21[0]*v41[2] ;
  w[2] = v21[0]*v41[1] - v21[1]*v41[0] ;
  s[1] = w[0]*w[0] + w[1]*w[1] + w[2]*w[2] ;
  // face 1 3 4 : cross product 31^41
  w[0] = v31[1]*v41[2] - v31[2]*v41[1] ;
  w[1] = v31[2]*v41[0] - v31[0]*v41[2] ;
  w[2] = v31[0]*v41[1] - v31[1]*v41[0] ;
  s[2] = w[0]*w[0] + w[1]*w[1] + w[2]*w[2] ;
  // face 2 3 4 : cross product 32^42
  w[0] = v32[1]*v42[2] - v32[2]*v42[1] ;
  w[1] = v32[2]*v42[0] - v32[0]*v42[2] ;
  w[2] = v32[0]*v42[1] - v32[1]*v42[0] ;
  s[3] = w[0]*w[0] + w[1]*w[1] + w[2]*w[2] ;
 
  double max_value = std::max({s[0], s[1], s[2], s[3]}) ;
  xlong = Det / sqrt(max_value) ;
 
}
 
void Structural_dynamic_mesh_model::setB0_(const DoubleTab& B)
{
 
  // Save B operator at initial step to compute transformation gradient with respect to initial configuration
 
  int pos = 0 ;
 
  for (int j=0; j < dimension; j++)
    {
      for (int i=0; i < nbn_ ; i++)
        {
          B0_(iel_,pos) = B(j,i) ;
          pos++ ;
        }
    }
 
}
 
void Structural_dynamic_mesh_model::load_behaviour_(StressMeasureKind smk, TangentOperatorKind tok)
{
  if (loaded_)
    return;
  else
    loaded_ = true;
 
  auto getMgisStressMeasureKind = [](StressMeasureKind k)
  {
    switch (k)
      {
      case StressMeasureKind::cauchy:
        return mgis::behaviour::FiniteStrainBehaviourOptions::CAUCHY;
      case StressMeasureKind::pk1:
        return mgis::behaviour::FiniteStrainBehaviourOptions::PK1;
      case StressMeasureKind::pk2:
        return mgis::behaviour::FiniteStrainBehaviourOptions::PK2;
      case StressMeasureKind::conjugate:
        return mgis::behaviour::FiniteStrainBehaviourOptions::CAUCHY;    // to avoid warning
      }
    return mgis::behaviour::FiniteStrainBehaviourOptions::CAUCHY;    // to avoid warning
  };
 
  auto getMgisTangentOperatorKind = [](TangentOperatorKind k)
  {
    switch (k)
      {
      case TangentOperatorKind::none:
        return mgis::behaviour::FiniteStrainBehaviourOptions::DSIG_DF;    // default mgis tangent operator kind
      case TangentOperatorKind::dCauchydF:
        return mgis::behaviour::FiniteStrainBehaviourOptions::DSIG_DF;
      case TangentOperatorKind::dPk2dE:
        return mgis::behaviour::FiniteStrainBehaviourOptions::DS_DEGL;
      case TangentOperatorKind::dPk1dF:
        return mgis::behaviour::FiniteStrainBehaviourOptions::DPK1_DF;
      }
    return mgis::behaviour::FiniteStrainBehaviourOptions::DSIG_DF;    // to avoid warning
  };
 
  mgisBehaviour_ = mgis::behaviour::load(
  {getMgisStressMeasureKind(smk), getMgisTangentOperatorKind(tok)}, l_, f_, mgis::behaviour::fromString(h_));
 
  // internal variables
  this->nbIvars_ = static_cast<int> (mgis::behaviour::getArraySize(mgisBehaviour_.isvs, hypothesis_));
 
  // external variables
  this->nbEvars_ = static_cast<int> (mgis::behaviour::getArraySize(mgisBehaviour_.esvs, hypothesis_));
 
  // K
  this->KSize_ = static_cast<int> (mgis::behaviour::getArraySize(mgisBehaviour_.gradients, hypothesis_) *
                                   mgis::behaviour::getArraySize(mgisBehaviour_.thermodynamic_forces, hypothesis_));
}
 
void Structural_dynamic_mesh_model::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
                                                        )
{
  mgisBehaviourData_.s0.gradients = const_cast<double *>(gradientOrStrain0);
  mgisBehaviourData_.s1.gradients = const_cast<double *>(gradientOrStrain1);
 
  mgisBehaviourData_.s0.thermodynamic_forces = stress;
  mgisBehaviourData_.s1.thermodynamic_forces = stress;
 
  mgisBehaviourData_.s0.internal_state_variables = ivar;
  mgisBehaviourData_.s1.internal_state_variables = ivar;
 
  mgisBehaviourData_.s0.external_state_variables = const_cast<double *>(evar);
  mgisBehaviourData_.s1.external_state_variables = const_cast<double *>(evar);
 
  mgisBehaviourData_.K = stiff;
 
  mgisBehaviourData_.speed_of_sound = vsound;
 
  int ierr = mgis::behaviour::integrate(mgisBehaviourData_, mgisBehaviour_);
 
  if (ierr == -1)
    {
      // Message::error("Integration of the behavior failed");
      Cerr << "Error: Integration of the behaviour with Mfront failed" << finl ;
      Process::exit() ;
    }
  else if (ierr == 0)
    {
      // Message::warning("Integration of the behavior unreliable");
      Cerr << "Warning: Integration of the behaviour with Mfront unreliable" << finl ;
    }
}
 
void Structural_dynamic_mesh_model::setLocalFields(const int elnodes[4], const int elem)
{
 
  iel_ = elem ;
 
  nbn_ = 3 ;
  if (dimension == 3) nbn_ = 4 ;
 
  for (int i=0; i<nbn_; i++)
    {
      int isom = elnodes[i] ;
      for (int j=0; j<dimension; j++)
        {
          xl_(i,j) = x(isom,j) ;
          ul_(i,j) = u(isom,j) ;
        }
    }
 
  int pos = 0 ;
  for (int j=0; j < dimension; j++)
    {
      for (int i=0; i < nbn_ ; i++)
        {
          B0l_(j,i) = B0_(iel_,pos) ;
          pos++ ;
        }
    }
 
}
 
void Structural_dynamic_mesh_model::setGlobalFields(const int elnodes[4], const double Pressure, const double VonMises)
{
 
  for (int i=0; i<nbn_; i++)
    {
      int isom = elnodes[i] ;
      for (int j=0; j<dimension; j++)
        {
          ff(isom,j) += fl_(i,j) ;
        }
    }
 
  meshPbPressure_(iel_) = Pressure ;
  meshPbVonMises_(iel_) = VonMises ;
 
}
 
double Structural_dynamic_mesh_model::computeCriticalDt(const double Vol, const double Xlong, const double cSound, const double dtMin, double& scaleMass)
{
 
  // Current density
  double currDensity = massElem_(iel_) / Vol ;
 
  // Speed of sound
  if (currDensity <= 0.)
    {
      Cerr << "Error: negative density in grid mesh motion" << finl ;
      Process::exit() ;
    }
 
  double cc = cSound ;
  scaleMass = 0. ;
  if (dtMin > 0.)
    {
      double dtc = Xlong / cc ;
      if (dtc < dtMin)
        {
          cc = Xlong / dtMin ;
          double newDensity = currDensity * pow(cSound,2) / pow(cc,2) ;
          scaleMass = (newDensity - currDensity) * Vol ;
        }
    }
 
  return Xlong / cc ;
 
}
 
void Structural_dynamic_mesh_model::computeForceFaces(const int nb_faces, const int nb_som_face, const IntTab& face_sommets)
{
 
  int i, j, s ;
  meshPbForceFace_ = 0. ;
  for (j=0; j<dimension; j++)
    {
      for (i=0; i<nb_faces; i++)
        {
          for (s=0; s<nb_som_face; s++)
            {
              meshPbForceFace_(i,j) += ff(face_sommets(i,s),j);
            }
          meshPbForceFace_(i,j)/=nb_som_face;
        }
    }
 
}
 
void Structural_dynamic_mesh_model::checkElemOrientation(int elnodes[4], const int elem)
{
 
  if (gridNStep == 0)
    {
      DoubleTab Jac(dimension, dimension) ;
      double Det, Vol, Xlong ;
      int nbn=0 ;
 
      switch (dimension)
        {
        case (2) :
          nbn = 3 ;
 
          break ;
        case (3) :
          nbn = 4 ;
 
          break ;
        }
 
      for (int i=0; i<nbn; i++)
        {
          int isom = elnodes[i] ;
          for (int j=0; j<dimension; j++)
            {
              xl_(i,j) = x(isom,j) ;
            }
        }
 
      Jac=0 ;
      for (int k=0; k < nbn; k++)
        {
          for (int j=0; j < dimension; j++)
            {
              for (int i=0; i < dimension ; i++)
                {
                  Jac(j,i) += Eta_(j,k) * xl_(k,i) ;
                }
            }
        }
 
      switch (dimension)
        {
        case (2) :
          Det = Jac(0,0) * Jac(1,1) - Jac(0,1) * Jac(1,0) ;
          triangleSurfLength_(Vol, Xlong, Det) ;
 
          break ;
        case (3) :
          Det= Jac(0,0) * (Jac(1,1) * Jac(2,2) - Jac(1,2) * Jac(2,1))
               - Jac(0,1) * (Jac(1,0) * Jac(2,2) - Jac(1,2) * Jac(2,0))
               + Jac(0,2) * (Jac(1,0) * Jac(2,1) - Jac(1,1) * Jac(2,0)) ;
          tetrahedronVolLength_(Vol, Xlong, Det) ;
 
          break ;
        }
 
      if (Vol <= 0.) invertNum_[elem] = 1 ;
 
    }
 
  if (invertNum_[elem] == 1)
    {
      // negative volume in mesh numbering, invert vertices numbering in elnodes
      int i1, i2, i3, i4 ;
 
      switch (dimension)
        {
        case (2) :
          i1 = elnodes[0] ;
          i2 = elnodes[1] ;
          i3 = elnodes[2] ;
          elnodes[0] = i1 ;
          elnodes[1] = i3 ;
          elnodes[2] = i2 ;
 
          break ;
        case (3) :
          i1 = elnodes[0] ;
          i2 = elnodes[1] ;
          i3 = elnodes[2] ;
          i4 = elnodes[3] ;
          elnodes[0] = i3 ;
          elnodes[1] = i2 ;
          elnodes[2] = i1 ;
          elnodes[3] = i4 ;
        }
    }
 
}
void Structural_dynamic_mesh_model::resumptionMesh(double tinit,DoubleTab& u_res_n, DoubleTab& v_res_n, DoubleTab& a_res_n,  DoubleTab& x_res_n, DoubleTab& B0_res_n, DoubleTab& Ft_res_n, DoubleTab& Stress_res_n)
{
  u=u_res_n;
  v=v_res_n;
  a=a_res_n;
  x=x_res_n;
 
  B0_=B0_res_n;
  Ft_=Ft_res_n;
  Stress_=Stress_res_n;
 
  gridTime=tinit;
  resumption = 1;
}
 
void Structural_dynamic_mesh_model::solveDynamicMeshProblem(const double temps, const DoubleTab& imposedVelocity, const IntVect& imposedVelocityTag,
                                                            DoubleTab& outputMeshVelocity, const int nbSom, const int nbElem, const int nbSomElem,
                                                            const IntTab& sommets, const int nbFace, const int nbSomFace, const IntTab& face_sommets, const Domaine_ALE& dom)
{
  DoubleTab x0( x) ; // Copy coordinates at the beginning of the step
  double tt =  gridTime ;
 
  double t0 = tt ;
  bool loopOnGridProblem = true ;
  int nstepCurr = 0 ;
 
  if ( configurationResetDt > 0)
    {
      if (tt >  gridResetTime)
        {
          doConfigurationReset=true ;
          gridResetTime +=  configurationResetDt ;
        }
    }
 
  // ICoCo + FSI : accelerated mode to activate within implicit iteration with subcycling
  if ( iCoCoImplicitIteration <= 0)  acceleratedSolution=0 ; // No iCoCo coupling or no possible acceleration at iteration 0
  else
    {
      if (! acceleratedSolutionEnabled)  acceleratedSolution=0 ; // Acceleration disabled after zero mesh motion at iteration 0
    }
 
  if ( acceleratedSolution == 1) Cerr << "Domaine_ALE::solveDynamicMeshProblem, accelerated solution with iCoCo implicit iteration: " <<  iCoCoImplicitIteration << finl ;
 
  while (loopOnGridProblem)
    {
      // Adjust the grid time step for smooth arrival at time = temps
      if ( gridDt >= temps - tt)
        {
          gridDt = temps - tt ;
          loopOnGridProblem = false ; // final time reached after this last step
        }
      else if ( gridDt >= 0.5 * (temps - tt)) // Avoid excessive time step variations
        {
          gridDt = 0.5 * (temps - tt) ;
        }
 
      double Dt ;
      if ( acceleratedSolution == 0)
        {
 
          Dt =  gridDt ;
 
          // Update mid-step velocities, displacements and coordinates
 
          for (int i=0; i<nbSom; i++)
            {
              int ii = dom.get_renum_som_perio(i) ; // to get imposed velocity from periodic boundaries if any
              for (int j=0; j<dimension; j++)
                {
                  vp(i,j) =  v(i,j) + 0.5 * Dt *  a(i,j) ;
                  if (imposedVelocityTag[i] == 1)  vp(i,j) = imposedVelocity(ii,j) ; // Apply imposed velocity from the boundary
 
                  double du = Dt *  vp(i,j) ;
                  u(i,j) += du ;
                  x(i,j) += du ;
                }
            }
        }
      else
        {
          // Accelerated mode
 
          Dt = temps - tt ; // One step only
          loopOnGridProblem = false ;
          x= xLast ;
          vp= vpLast ;
 
          // Update mid-step velocities, displacements and coordinates for nodes on FSI boundary only
 
          for (int i=0; i<nbSom; i++)
            {
              int ii = dom.get_renum_som_perio(i) ; // to get imposed velocity from periodic boundaries if any
              for (int j=0; j<dimension; j++)
                {
                  if (imposedVelocityTag[i] == 1)
                    {
                      vp(i,j) = imposedVelocity(ii,j) ; // Apply imposed velocity from the boundary
                      double du = Dt *  vp(i,j) ;
                      u(i,j) += du ;
                      x(i,j) = x0(i,j) + du ;
                    }
                }
            }
        }
      if ( iCoCoImplicitIteration == 0 && !loopOnGridProblem)
        {
          xLast =  x ;
          vpLast =  vp ;
          dtLast = Dt ;
        }
 
      tt += Dt ;
 
      // Compute internal forces
      ff = 0. ;
      int nbn =  getNbNodesPerElem() ;
      int elnodes[4] ;
      double volume ;
      double xlong ;
      double cSound ;
      double density =  getDensity() ;
      double dtm =  getGridDtMin() ;
      double dts ;
      double scaleMass ;
      double totalScaleMass = 0 ;
      double dtmin ;
      double pressure ;
      double vonmises ;
      double mm = 0 ;
 
      dtmin = 1.E12 ; // Initialize dtmin at <huge>
      if (dtm > 0.)  nodalScaleMass = 0. ; // Initialize nodal additional masses if needed
 
      if (! isMassBuilt)  totalMass = 0. ;
 
      for (int elem=0; elem<nbElem; elem++)
        {
          for (int i=0; i< nbn; i++) elnodes[i] = sommets(elem,i) ;
 
          if ( acceleratedSolution == 0)  skipStressComputation = false ;
          else
            {
              skipStressComputation = true ; // No stress computation by default in accelerated mode
              for (int i=0; i< nbn; i++)
                {
                  int ii=elnodes[i];
                  if (imposedVelocityTag[ii] == 1)  skipStressComputation = false ; // Compute stress for cells connected to the moving boundary
                }
            }
 
          checkElemOrientation(elnodes, elem) ; // check orientation to ensure a positive element volume
 
          setLocalFields(elnodes, elem) ; // x, u, B0 global to local + elem id
 
          computeInternalForces(volume, xlong, cSound, pressure, vonmises) ; // local force computation on element elem
 
          setGlobalFields(elnodes, pressure, vonmises) ; // ff back to global + store elem pressure and vonmises for postprocessing
 
          if (! isMassBuilt)
            {
              setMassElem(volume * density) ;
              totalMass += volume * density ;
              for (int i=0; i< nbn; i++)
                {
                  int ii = elnodes[i] ;
                  mass[ii] += volume * density / nbn ;
                }
            }
 
          // Set next time step
          dts =  computeCriticalDt(volume, xlong, cSound, dtm, scaleMass) ;
          dtmin = std::min(dtmin, dts) ;
 
          if (scaleMass > 0.)
            {
              totalScaleMass += scaleMass ;
              for (int i=0; i< nbn; i++)
                {
                  int ii = elnodes[i] ;
                  nodalScaleMass[ii] += scaleMass / nbn ;
                }
            }
        }
 
      if ( doConfigurationReset)  doConfigurationReset=false ;
 
      gridDt = mp_min(dtmin) ;
 
      MD_Vector_tools::echange_espace_virtuel( ff, MD_Vector_tools::EV_SOMME_ECHANGE) ;
      if (! isMassBuilt)
        {
          MD_Vector_tools::echange_espace_virtuel( mass, MD_Vector_tools::EV_SOMME_ECHANGE) ;
          totalMass = mp_sum( totalMass) ;
          isMassBuilt = true ;
        }
 
      if (dtm > 0.)
        {
          totalScaleMass = mp_sum(totalScaleMass) ;
          if (totalScaleMass > 0.)
            {
              MD_Vector_tools::echange_espace_virtuel( nodalScaleMass, MD_Vector_tools::EV_SOMME_ECHANGE) ;
              if ( maxAddedMassRatio > 0.)
                {
                  double r = totalScaleMass /  totalMass ;
                  if (r >  maxAddedMassRatio)  doConfigurationReset=true ;
                }
            }
        }
 
      // Compute accelerations and full step velocities
      double rhs ;
      double den ;
      double d =  getDampingCoefficient() ;
 
      double dtUpdate=Dt ;
      if ( acceleratedSolution == 1) dtUpdate= dtLast ; // Accelerated mode: acceleration and velocity update must be done with the last stable dt
 
      for (int i=0; i<nbSom; i++)
        {
          mm =  mass[i] ;
          if (dtm > 0.) mm +=  nodalScaleMass[i] ;
          for (int j=0; j<dimension; j++)
            {
              rhs = - ff(i,j) - d * mm *  vp(i,j) ;
              den = mm * (1. + 0.5 * d * dtUpdate) ;
              a(i,j) = rhs / den ;
              v(i,j) =  vp(i,j) + 0.5 * dtUpdate *  a(i,j) ;
            }
        }
 
      gridNStep += 1 ;
      applyDtCoefficient() ;
      gridTime = tt ;
      nstepCurr += 1 ;
      if (dtm > 0.)
        {
          if ( acceleratedSolution == 0) Cerr << "Grid dynamic problem, internal step: "<< nstepCurr << ", dt= " << Dt << ", dt_stab= " <<  gridDt << ", [mass scaling] added_mass_ratio= "<< totalScaleMass /  totalMass << ", time= " << tt << ", target fluid time= " << temps << finl ;
          else Cerr << "Grid dynamic problem, accelerated internal step: "<< nstepCurr << ", dt= " << Dt << ", [mass scaling] added_mass_ratio= "<< totalScaleMass /  totalMass << ", time= " << tt << ", target fluid time= " << temps << finl ;
        }
      else
        {
          if ( acceleratedSolution == 0) Cerr << "Grid dynamic problem, internal step: "<< nstepCurr << ", dt= " << Dt << ", dt_stab= " <<  gridDt << ", time= " << tt << ", target fluid time= " << temps << finl ;
          else Cerr << "Grid dynamic problem, accelerated internal step: "<< nstepCurr << ", dt= " << Dt  << ", time= " << tt << ", target fluid time= " << temps << finl ;
        }
 
    } // End time loop on grid problem
 
  if (tt > t0)
    {
      for (int i=0; i<nbSom; i++)
        {
          for (int j=0; j<dimension; j++)
            {
              outputMeshVelocity(i,j) = ( x(i,j) - x0(i,j)) / (tt - t0) ;
            }
        }
 
    }
  else
    {
      outputMeshVelocity = 0. ;
    }
 
  outputMeshVelocity.echange_espace_virtuel();
 
  // Compute forces at face centers for postprocessing only
  computeForceFaces(nbFace,nbSomFace,face_sommets) ;
}