PASSMONL DEVS: stress prediction implementation

Author: cpp2cool

femutils/FemUtils.h

@@ -769,6 +769,12 @@ class Tensor2
     this->add(tb);
   }
 
+  //! Add b to this Tensor2
+  ARCCORE_HOST_DEVICE void add(const RealVector<6>& b) {
+    for (Arcane::Int32 i = 0; i < 6; ++i)
+      m_vec(i) += b(i);
+  }
+
   //! Substract b to this Tensor2
   ARCCORE_HOST_DEVICE void sub(const Tensor2& b) {
     for (Arcane::Int32 i = 0; i < 6; ++i)
@@ -780,6 +786,12 @@ class Tensor2
     this->sub(tb);
   }
 
+  //! Substract b to this Tensor2
+  ARCCORE_HOST_DEVICE void sub(const RealVector<6>& b) {
+    for (Arcane::Int32 i = 0; i < 6; ++i)
+      m_vec(i) -= b(i);
+  }
+
   //! Scalar multiplication: Tensor * scalar
   ARCCORE_HOST_DEVICE Tensor2 operator*(Real scalar) const {
     Tensor2 result;

femutils/GaussDoFsOnCells.cc

@@ -57,6 +57,7 @@ class GaussDoFsOnCells::Impl
   VariableDoFReal* m_gauss_weight = nullptr;
   VariableDoFReal* m_gauss_jacobian = nullptr;
   VariableDoFArrayReal*  m_gauss_law_param = nullptr;
+  VariableDoFArrayReal* m_gauss_law_history_param = nullptr;
 
 /*  VariableDoFArrayTensor2* m_gauss_stress = nullptr;
   VariableDoFArrayTensor2* m_gauss_strain = nullptr;
@@ -90,6 +91,7 @@ GaussDoFsOnCells::
   delete m_p->m_gauss_strain;
   delete m_p->m_gauss_strain_plastic;
   delete m_p->m_gauss_law_param;
+  delete m_p->m_gauss_law_history_param;
   delete m_p;
 }
 
@@ -171,6 +173,7 @@ initialize(IMesh* mesh, Int32 max_nb_gauss_per_cell)
   m_p->m_gauss_shape = new VariableDoFArrayReal(VariableBuildInfo(mesh, "GaussShape", "GaussCellFamily"));
   m_p->m_gauss_shapederiv = new VariableDoFArrayReal3(VariableBuildInfo(mesh, "GaussShapeDeriv", "GaussCellFamily"));
   m_p->m_gauss_law_param = new VariableDoFArrayReal(VariableBuildInfo(mesh, "GaussLawParam", "GaussCellFamily"));
+  m_p->m_gauss_law_history_param = new VariableDoFArrayReal(VariableBuildInfo(mesh, "GaussLawHistoryParam", "GaussCellFamily"));
 
 /*  m_p->m_gauss_stress = new VariableDoFArrayTensor2(VariableBuildInfo(mesh, "GaussStress", "GaussCellFamily"));
   m_p->m_gauss_strain = new VariableDoFArrayTensor2(VariableBuildInfo(mesh, "GaussStrain", "GaussCellFamily"));
@@ -236,6 +239,14 @@ VariableDoFArrayReal& GaussDoFsOnCells::
 gaussLawParam(){
   return *m_p->m_gauss_law_param;
 }
+
+/*---------------------------------------------------------------------------*/
+/*---------------------------------------------------------------------------*/
+VariableDoFArrayReal& GaussDoFsOnCells::
+gaussLawHistoryParam(){
+  return *m_p->m_gauss_law_param;
+}
+
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 

femutils/GaussDoFsOnCells.h

@@ -70,6 +70,7 @@ class GaussDoFsOnCells
   Arcane::VariableDoFReal& gaussWeight();
   Arcane::VariableDoFReal& gaussJacobian();
   Arcane::VariableDoFArrayReal& gaussLawParam();
+  Arcane::VariableDoFArrayReal& gaussLawHistoryParam();
 
 /*  Arcane::VariableDoFArrayTensor2& gaussStress();
   Arcane::VariableDoFArrayTensor2& gaussStrain();

modules/passmonl/LawDispatcher.cc

@@ -96,15 +96,15 @@ LawDispatcher::LawDispatcher(TypesNLDynamic::eLawType law_type, bool default_par
 
     switch(law_type){
 
-    case TypesNLDynamic::UNKNOWN:
-    case TypesNLDynamic::HOOKE: m_nb_law_param = 2;
+    case TypesNLDynamic::HOOKE: m_nb_law_param = 2; m_nb_law_history_param = 0;
       break;
 
     case TypesNLDynamic::DRUCKP:
-    case TypesNLDynamic::MOHRC: m_nb_law_param = 7;
+    case TypesNLDynamic::MOHRC: m_nb_law_param = 7; m_nb_law_history_param = 1;
       break;
 
-    default: m_nb_law_param = 2;
+    case TypesNLDynamic::UNKNOWN:
+    default: m_nb_law_param = 2; m_nb_law_history_param = 0;
       break;
     }
 }
@@ -171,10 +171,11 @@ Tensor4 LawDispatcher::computeTangentTensor(const Tensor2& sig) {
 
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
-bool LawDispatcher::initState(const Tensor2& sig, RealConstArrayView& history_vars)
+bool LawDispatcher::initState(const Tensor2& sig)
 {
     if (sig == Tensor2::zero()) return true;
 
+    ConstArrayView<Real> history_vars = m_history_vars.constView();
     m_history_vars = initHistoryVars(history_vars);
 
     auto f = m_init_state[m_law_type];
@@ -260,4 +261,5 @@ void ReadLawBlock(istream& is, Integer nblock) {
         stop = true;
     }
   } while(!stop);
-}
+}
+

modules/passmonl/LawDispatcher.h

@@ -45,7 +45,7 @@ class LawDispatcher {
 public:
     Tensor4 computeElastTensor(const Tensor2& sig);
     Tensor4 computeTangentTensor(const Tensor2& sig);
-    bool initState(const Tensor2& sig, RealConstArrayView& history_vars);
+    bool initState(const Tensor2& sig);
     RealUniqueArray initHistoryVars(RealConstArrayView& history_vars);
     void computeStress(bool isRef);
     RealUniqueArray initConsts(RealConstArrayView& law_params);
@@ -69,6 +69,7 @@ class LawDispatcher {
     void	setStrainIncrement(const Tensor2&);
 
     [[nodiscard]] Integer getNbLawParam() const { return m_nb_law_param; }
+    [[nodiscard]] Integer getNbLawHistoryParam() const { return m_nb_law_history_param; }
 
     void  setLambda(Real lambda) { m_Lambda = lambda; }
     void  setMu(Real mu) { m_Mu = mu; }
@@ -77,9 +78,12 @@ class LawDispatcher {
     void setLawParams(const RealUniqueArray& lawparams) {
       m_law_params = lawparams.clone();
     }
+    void setLawHistoryParams(const RealUniqueArray& lawhistparams) {
+      m_history_vars = lawhistparams.clone();
+    }
 
 private:
- std::function<Tensor4(RealConstArrayView& law_params, RealArrayView& history_vars, Tensor2& sig, Tensor2& eps, Tensor2& epsp, Tensor2& dsig,
+    std::function<Tensor4(RealConstArrayView& law_params, RealArrayView& history_vars, Tensor2& sig, Tensor2& eps, Tensor2& epsp, Tensor2& dsig,
                        const Tensor2& deps, bool isRef)> m_compute_stress[NB_LAW_TYPE];
     std::function<Tensor4(RealConstArrayView& law_params, const Tensor2& sig)> m_compute_elast_tensor[NB_LAW_TYPE];
     std::function<Tensor4(RealConstArrayView& law_params, RealArrayView& history_vars,
@@ -106,6 +110,7 @@ class LawDispatcher {
     String m_name{};
     bool m_default{true};
     Integer m_nb_law_param{2};
+    Integer m_nb_law_history_param{0};
 };
 
 #endif //PASSMO_LAWDISPATCHER_H

modules/passmonl/NLDynamic.axl

@@ -505,6 +505,9 @@
         <simple name="nb-law-param" type="integer" default = "2" optional="true">
           <description>Number of law parameters for the selected constitutive model (including the 2 elastic coef.)</description>
         </simple>
+        <simple name="nb-law-hist-param" type="integer" default = "0" optional="true">
+          <description>Number of law history parameters for the selected constitutive model</description>
+        </simple>
         <simple name="i-law-param" type="integer" default = "0" optional="true">
           <description>ith law parameters block in the constitutive models file</description>
         </simple>