PASSMO-NL: stress prediction at Gauss points

Author: cpp2cool

femutils/ArcaneFemFunctions.h

@@ -28,7 +28,7 @@
 
 using namespace Arcane;
 using namespace Arcane::FemUtils;
-
+Real REL_PREC{1.0e-15};
 /*---------------------------------------------------------------------------*/
 /**
  * @brief Contains various functions & operations related to FEM calculations.
@@ -1436,7 +1436,7 @@ class ArcaneFemFunctions
   /**
    * @brief Provides methods based on the Dispatcher mechanism available in Arcane,
    * allowing to compute FEM methods without declaring the FE entity type
-   * (coming fro PASSMO).
+   * (coming from PASSMO).
    */
   /*---------------------------------------------------------------------------*/
   class CellFEMDispatcher
@@ -1464,7 +1464,7 @@ class ArcaneFemFunctions
    * This class includes static methods for computing shape functions and their
    * derivatives depending on finite element types. These methods are used within
    * the Dispatcher mechanism available in Arcane through the class
-   * CellFEMDispatcher (coming fro PASSMO).
+   * CellFEMDispatcher (coming from PASSMO).
    */
   /*---------------------------------------------------------------------------*/
   class FemShapeMethods
@@ -2560,7 +2560,7 @@ class ArcaneFemFunctions
    * @brief Provides methods for Gauss quadrature.
    *
    * This class includes static methods for computing Gauss-Legendre integration
-   * depending on finite element types (coming fro PASSMO).
+   * depending on finite element types (coming from PASSMO).
    */
   /*---------------------------------------------------------------------------*/
   class FemGaussQuadrature

femutils/FemUtils.h

@@ -30,7 +30,6 @@
 #include <arcane/MeshVariableArrayRef.h>
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
-extern Arcane::Real REL_PREC;
 
 struct Real4
 {
@@ -631,8 +630,9 @@ class Tensor2
 
   //! Define the == operator
   ARCCORE_HOST_DEVICE bool operator==(const Tensor2& vec) {
+    Real eps{1.0e-15};
     for (Arcane::Int32 i = 0; i < 6; ++i) {
-      if (fabs(m_vec(i) - vec(i)) > REL_PREC)
+      if (fabs(m_vec(i) - vec(i)) > eps)
         return false;
     }
     return true;
@@ -792,6 +792,14 @@ class Tensor2
     return sqrt(dot(m_vec,m_vec));
   }
 
+  //! Function to convert Real3x3 matrix (symmetric) to Tensor
+  ARCCORE_HOST_DEVICE void fromReal3x3ToTensor2(const Real3x3& mat) {
+    for (Arcane::Int32 i = 0; i < 3; i++) (*this)(i) = mat[i][i];
+    for (Arcane::Int32 i = 3; i < 5; i++) (*this)(i) = mat[0][i - 2];
+    (*this)(5) = mat[1][2];
+  }
+
+
   //! Friend function for scalar multiplication: scalar * Tensor
   ARCCORE_HOST_DEVICE friend Tensor2 operator*(Real scalar, const Tensor2& vector) {
     Tensor2 result;

modules/passmo/ElastodynamicModule.cc

@@ -29,6 +29,7 @@
 /*---------------------------------------------------------------------------*/
 using namespace Arcane;
 using namespace Arcane::FemUtils;
+extern Real REL_PREC;
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 ElastodynamicModule::ElastodynamicModule(const ModuleBuildInfo& mbi)
@@ -521,14 +522,16 @@ compute(){
   // Update the nodal variable according to the integration scheme (e.g. Newmark)
   _updateNewmark();
 
-
+/*
   if (t < tf) {
     if (t + dt > tf) {
       dt = tf - t;
       m_global_deltat = dt;
     }
   }
   else {
+*/
+  if (t >= tf) {
     // Save/Check results
     _checkResultFile();
     subDomain()->timeLoopMng()->stopComputeLoop(true);
@@ -546,7 +549,7 @@ _checkResultFile()
   if (filename.empty())
     return;
   const double epsilon = 1.0e-4;
-  const double min_value_to_test = 1.0e-10;
+  const double min_value_to_test = 1.0e-16;
   Arcane::FemUtils::checkNodeResultFile(traceMng(), filename, m_displ, epsilon, min_value_to_test);
 }
 /*---------------------------------------------------------------------------*/

modules/passmonl/CMakeLists.txt

@@ -4,14 +4,14 @@ add_executable(Passmonl
   NLDynamicModule.h
   NLDynamicModule.cc
   main.cc
+  utilFEM.h
+  utilFEM.cc
+  LawDispatcher.h
+  LawDispatcher.cc
   analytic_func.cc
   analytical_func.h
-  utilFEM.h
   hooke.cc
   druckerp.cc
-  LawDispatcher.h
-  LawDispatcher.cc
-  utilFEM.cc
 )
 
 arcane_generate_axl(NLDynamic)

modules/passmonl/LawDispatcher.cc

@@ -93,6 +93,20 @@ LawDispatcher::LawDispatcher(TypesNLDynamic::eLawType law_type, bool default_par
     //! Initialize intern useful constants allowing to the material constitutive model type
     m_init_consts[TypesNLDynamic::HOOKE] = HookeInitConsts;
     m_init_consts[TypesNLDynamic::DRUCKP] = DruckPInitConsts;
+
+    switch(law_type){
+
+    case TypesNLDynamic::UNKNOWN:
+    case TypesNLDynamic::HOOKE: m_nb_law_param = 2;
+      break;
+
+    case TypesNLDynamic::DRUCKP:
+    case TypesNLDynamic::MOHRC: m_nb_law_param = 7;
+      break;
+
+    default: m_nb_law_param = 2;
+      break;
+    }
 }
 
 /////////////////////////////////////////////////////////////////////////////
@@ -228,24 +242,22 @@ RealUniqueArray LawDispatcher::initConsts(RealConstArrayView& law_params)
 
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
-void ReadLawBlock(istream& is, Integer nblock){
+void ReadLawBlock(istream& is, Integer nblock) {
+
   char c[500];
   String str;
-  bool stop { false};
-  String s_end_b{"</"};
-
-  for (Integer i = 0; i < nblock; ++i){
-
-    if (!i) {
-      is >> str;
-      is.getline(c, 500); // read until end of line "\n"
-      stop = (str.contains(s_end_b));
-    }
-
-    while(!stop) {
-      is >> str;
-      is.getline(c, 500); // read until end of line "\n"
-      stop = (str.contains(s_end_b));
+  bool stop{ false };
+  String str_nblock = String::fromNumber(nblock);
+  String sblock = String("<") + str_nblock + String(">");
+
+  do {
+
+    is.getline(c, 500); // read whole line
+    std::istringstream iss(c);
+    String token;
+    while (iss >> token) {
+      if (token.contains(sblock))
+        stop = true;
     }
-  }
+  } while(!stop);
 }

modules/passmonl/LawDispatcher.h

@@ -68,10 +68,15 @@ class LawDispatcher {
     [[nodiscard]] Tensor2	getStrainIncrement() const;
     void	setStrainIncrement(const Tensor2&);
 
+    [[nodiscard]] Integer getNbLawParam() const { return m_nb_law_param; }
+
     void  setLambda(Real lambda) { m_Lambda = lambda; }
     void  setMu(Real mu) { m_Mu = mu; }
     void  setName(const String& name) { m_name = name; }
     void  setDefault(bool is_default) { m_default = is_default; }
+    void setLawParams(const RealUniqueArray& lawparams) {
+      m_law_params = lawparams.clone();
+    }
 
 private:
     std::function<Tensor4(RealConstArrayView& law_params, RealArrayView& history_vars, Tensor2& sig, Tensor2& eps, Tensor2& epsp, Tensor2& dsig,
@@ -84,6 +89,7 @@ class LawDispatcher {
     std::function<RealUniqueArray(RealConstArrayView& law_params)> m_init_consts[NB_LAW_TYPE];
 
     Tensor2 m_sig{};
+    Tensor2 m_sign{};
     Tensor2 m_eps{};
     Tensor2 m_dsig{};
     Tensor2 m_epsp{};
@@ -97,6 +103,7 @@ class LawDispatcher {
     TypesNLDynamic::eLawType m_law_type{TypesNLDynamic::HOOKE};
     String m_name{};
     bool m_default{true};
+    Integer m_nb_law_param{2};
 };
 
 #endif //PASSMO_LAWDISPATCHER_H

modules/passmonl/NLDynamicModule.cc

@@ -24,6 +24,7 @@
 #include "IDoFLinearSystemFactory.h"
 #include "ArcaneFemFunctions.h"
 #include "NLDynamicModule.h"
+#include "LawDispatcher.h"
 
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
@@ -455,9 +456,9 @@ _applyInitialCellConditions(){
   for (Integer i = 0, nb = options()->lawModel().size(); i < nb; ++i){
 
     CellGroup cell_group = options()->lawModel[i]->cellGroup();
-    Integer ilaw = options()->lawModel[i]->iLawParam();
-    TypesNLDynamic::eLawType lawtyp = options()->lawModel[i]->lawType();
-    Integer nblaw = options()->lawModel[i]->nbLawParam();
+    auto ilaw = options()->lawModel[i]->iLawParam();
+    auto lawtyp = options()->lawModel[i]->lawType();
+    auto nblaw = options()->lawModel[i]->nbLawParam();
     m_nb_law_param = math::max(m_nb_law_param, nblaw);
 
     bool is_default = (m_law_param_file.empty() || nblaw == 2);
@@ -552,6 +553,7 @@ _applyInitialCellConditions(){
 void NLDynamicModule::
 _initGaussStep()
 {
+  Real eps{1.0e-15};
   auto gauss_point(m_gauss_on_cells.gaussCellConnectivityView());
 
   VariableDoFReal& gauss_jacobian(m_gauss_on_cells.gaussJacobian());
@@ -597,7 +599,7 @@ _initGaussStep()
       else
         jacobian = ArcaneFemFunctions::MeshOperation::computeLengthEdge2(cell, m_node_coord) / 2.;
 
-      if (fabs(jacobian) < REL_PREC) {
+      if (fabs(jacobian) < eps) {
         ARCANE_FATAL("Cell jacobian is null");
       }
       gauss_jacobian[gauss_pti] = jacobian;
@@ -1583,6 +1585,7 @@ Real3x3 NLDynamicModule::
 _computeJacobian(const ItemWithNodes& cell,const Int32& ig, const RealUniqueArray& vec, Real& jacobian) {
 
   auto	n = cell.nbNode();
+  Real eps{1.0e-15};
 
   // Jacobian matrix computed at the integration point
   Real3x3	jac;
@@ -1602,17 +1605,16 @@ _computeJacobian(const ItemWithNodes& cell,const Int32& ig, const RealUniqueArra
   }
 
   Int32 ndim = ArcaneFemFunctions::MeshOperation::getGeomDimension(cell);
-  //  if (NDIM == 3)
+
   if (ndim == 3)
     jacobian = math::matrixDeterminant(jac);
 
-  //  else if (NDIM == 2) {
   else if (ndim == 2)
     jacobian = jac.x.x * jac.y.y - jac.x.y * jac.y.x;
   else
     jacobian = ArcaneFemFunctions::MeshOperation::computeLengthEdge2(cell, m_node_coord) / 2.;
 
-  if (fabs(jacobian) < REL_PREC) {
+  if (fabs(jacobian) < eps) {
     ARCANE_FATAL("Cell jacobian is null");
   }
   return jac;
@@ -2343,19 +2345,73 @@ _assembleLinearRHS(){
 /*---------------------------------------------------------------------------*/
 // ! Stress prediction
 //
-bool NLDynamicModule::stress_prediction(bool init, bool isRef)
+void NLDynamicModule::stress_prediction(bool init, bool isRef)
 {
-  bool	stop = false;
-  return stop;
+ auto gauss_point(m_gauss_on_cells.gaussCellConnectivityView());
+
+  VariableDoFArrayReal3x3& gauss_stress(m_gauss_on_cells.gaussStress());
+  VariableDoFArrayReal3x3& gauss_strain(m_gauss_on_cells.gaussStrain());
+  VariableDoFArrayReal3x3& gauss_strain_plastic(m_gauss_on_cells.gaussStrainPlastic());
+  VariableDoFArrayReal& gauss_lawparam(m_gauss_on_cells.gaussLawParam());
+
+  ENUMERATE_CELL (icell, allCells()) {
+    const Cell& cell = *icell;
+    auto cell_type = cell.type();
+    auto cell_nbnod = cell.nbNode();
+    Int32 numcell = cell.localId();
+
+    auto is_default = m_default_law[cell];
+    auto lawtyp = static_cast<TypesNLDynamic::eLawType>(m_law[cell]);
+    auto ilaw = m_iparam_law[cell];
+    LawDispatcher cell_law(lawtyp, is_default);
+    auto nblaw = cell_law.getNbLawParam();
+    RealUniqueArray lawparams(nblaw);
+
+    auto cell_nbgauss = ArcaneFemFunctions::FemGaussQuadrature::getNbGaussPointsfromOrder(cell_type, ninteg);
+
+    for (Int32 ig = 0; ig < cell_nbgauss; ++ig) {
+      DoFLocalId gauss_pti = gauss_point.dofId(cell, ig);
+      Int32 gaussnum = gauss_pti.localId();
+
+      Tensor2 epsn;
+      epsn.fromReal3x3ToTensor2(gauss_strain[gauss_pti][1]);
+      Tensor2 sign;
+      sign.fromReal3x3ToTensor2(gauss_stress[gauss_pti][1]);
+      Tensor2 epspn;
+      epspn.fromReal3x3ToTensor2(gauss_strain_plastic[gauss_pti][1]);
+      Tensor2 deps;
+      {
+        // calcul deps à faire avec unodes
+      }
+
+      for (Int32 ip = 0; ip < nblaw; ++ip) {
+        lawparams[ip] = gauss_lawparam[gauss_pti][ip];
+      }
+      cell_law.setLawParams(lawparams);
+      cell_law.setStrain(epsn);
+      cell_law.setStress(sign);
+      cell_law.setPlasticStrain(epspn);
+      cell_law.setStrainIncrement(deps);
+
+      // bool isRef = false at prediction
+      // => it avoids computing tangent stiffness operator
+      cell_law.computeStress(isRef);
+
+      // Current plastic strains and stresses have been updated by the law
+      gauss_strain_plastic[gauss_pti][2] = fromTensor2Real3x3(cell_law.getPlasticStrain());
+      gauss_stress[gauss_pti][2] = fromTensor2Real3x3(cell_law.getStress());
+      gauss_strain[gauss_pti][2] = fromTensor2Real3x3(epsn + deps);
+
+      }
+    }
 }
+
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 // ! Stress correction
 //
-bool NLDynamicModule::stress_correction(bool is_converge)
+void NLDynamicModule::stress_correction(bool isRef)
 {
-  bool stop = false;
-  return stop;
 }
 
 /*---------------------------------------------------------------------------*/