Holzapfel Ogden Modified Anisotropy (HO-ma) Constitutive Model

Code/Source/svFSI/ComMod.h

@@ -333,6 +333,9 @@ class fibStrsType
     // Constant steady value
     double g = 0.0;
 
+    // Cross fiber stress parameter
+    double eta_s = 0.0;
+
     // Unsteady time-dependent values
     fcType gt;
 };

Code/Source/svFSI/Parameters.cpp

@@ -467,7 +467,8 @@ const std::string ConstitutiveModelParameters::xml_element_name_ = "Constitutive
 // [TODO] Should use the types defined in consts.h.
 const std::string ConstitutiveModelParameters::GUCCIONE_MODEL = "Guccione";
 const std::string ConstitutiveModelParameters::HGO_MODEL = "HGO";
-const std::string ConstitutiveModelParameters::HOLZAPFEL_MODEL = "Holzapfel";
+const std::string ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL = "HolzapfelOgden";
+const std::string ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL = "HolzapfelOgden-ModifiedAnisotropy";
 const std::string ConstitutiveModelParameters::LEE_SACKS = "Lee-Sacks";
 const std::string ConstitutiveModelParameters::NEOHOOKEAN_MODEL = "neoHookean";
 const std::string ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL = "stVenantKirchhoff";
@@ -479,7 +480,9 @@ const std::map<std::string, std::string> ConstitutiveModelParameters::constituti
 
   {ConstitutiveModelParameters::HGO_MODEL, ConstitutiveModelParameters::HGO_MODEL},
 
-  {ConstitutiveModelParameters::HOLZAPFEL_MODEL, ConstitutiveModelParameters::HOLZAPFEL_MODEL},
+  {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL, ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL},
+
+  {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL, ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL},
 
   {ConstitutiveModelParameters::LEE_SACKS, ConstitutiveModelParameters::LEE_SACKS},
 
@@ -498,7 +501,8 @@ using SetConstitutiveModelParamMapType = std::map<std::string, std::function<voi
 SetConstitutiveModelParamMapType SetConstitutiveModelParamMap = {
   {ConstitutiveModelParameters::GUCCIONE_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->guccione.set_values(params);}},
   {ConstitutiveModelParameters::HGO_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel_gasser_ogden.set_values(params);}},
-  {ConstitutiveModelParameters::HOLZAPFEL_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel.set_values(params);}},
+  {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel.set_values(params);}},
+  {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel.set_values(params);}},
   {ConstitutiveModelParameters::LEE_SACKS, [](CmpType cp, CmpXmlType params) -> void {cp->lee_sacks.set_values(params);}},
   {ConstitutiveModelParameters::NEOHOOKEAN_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->neo_hookean.set_values(params);}},
   {ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->stvenant_kirchhoff.set_values(params);}},
@@ -510,6 +514,8 @@ using PrintConstitutiveModelParamMapType = std::map<std::string, std::function<v
 PrintConstitutiveModelParamMapType PrintConstitutiveModelParamMap = {
   {ConstitutiveModelParameters::GUCCIONE_MODEL, [](CmpType cp) -> void {cp->guccione.print_parameters();}},
   {ConstitutiveModelParameters::HGO_MODEL, [](CmpType cp) -> void {cp->holzapfel_gasser_ogden.print_parameters();}},
+  {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MODEL, [](CmpType cp) -> void {cp->holzapfel.print_parameters();}},
+  {ConstitutiveModelParameters::HOLZAPFEL_OGDEN_MA_MODEL, [](CmpType cp) -> void {cp->holzapfel.print_parameters();}},
   {ConstitutiveModelParameters::LEE_SACKS, [](CmpType cp) -> void {cp->lee_sacks.print_parameters();}},
   {ConstitutiveModelParameters::NEOHOOKEAN_MODEL, [](CmpType cp) -> void {cp->neo_hookean.print_parameters();}},
   {ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL, [](CmpType cp) -> void {cp->stvenant_kirchhoff.print_parameters();}},

Code/Source/svFSI/Parameters.h

@@ -517,7 +517,8 @@ class ConstitutiveModelParameters : public ParameterLists
     // Model types supported.
     static const std::string GUCCIONE_MODEL;
     static const std::string HGO_MODEL;
-    static const std::string HOLZAPFEL_MODEL;
+    static const std::string HOLZAPFEL_OGDEN_MODEL;
+    static const std::string HOLZAPFEL_OGDEN_MA_MODEL;
     static const std::string LEE_SACKS;
     static const std::string NEOHOOKEAN_MODEL;
     static const std::string STVENANT_KIRCHHOFF_MODEL;

Code/Source/svFSI/VtkData.cpp

@@ -926,6 +926,24 @@ void VtkVtuData::copy_point_data(const std::string& data_name, Vector<int>& mesh
   }
 }
 
+/// @brief Copy points into the given array.
+//
+void VtkVtuData::copy_points(Array<double>& points)
+{ 
+  auto vtk_points = impl->vtk_ugrid->GetPoints();
+  auto num_points = vtk_points->GetNumberOfPoints();
+  Array<double> points_array(3, num_points);
+
+  double point[3];
+  for (int i = 0; i < num_points; i++) {
+    vtk_points->GetPoint(i, point);
+    points(0,i) = point[0];
+    points(1,i) = point[1];
+    points(2,i) = point[2];
+  } 
+
+  return;
+}
 
 bool VtkVtuData::has_point_data(const std::string& data_name)
 {

Code/Source/svFSI/VtkData.h

@@ -60,6 +60,7 @@ class VtkData {
 
     virtual bool has_point_data(const std::string& data_name) = 0;
 
+    virtual void copy_points(Array<double>& points) = 0;
     virtual void copy_point_data(const std::string& data_name, Array<double>& mesh_data) = 0;
     virtual void copy_point_data(const std::string& data_name, Vector<double>& mesh_data) = 0;
     virtual void write() = 0;
@@ -119,6 +120,7 @@ class VtkVtuData : public VtkData {
     virtual int num_points();
     virtual void read_file(const std::string& file_name);
 
+    void copy_points(Array<double>& points);
     void copy_point_data(const std::string& data_name, Array<double>& mesh_data);
     void copy_point_data(const std::string& data_name, Vector<double>& mesh_data);
     void copy_point_data(const std::string& data_name, Vector<int>& mesh_data);

Code/Source/svFSI/consts.cpp

@@ -86,7 +86,10 @@ const std::map<std::string,ConstitutiveModelType> constitutive_model_name_to_typ
   {"Gucci", ConstitutiveModelType::stIso_Gucci},
 
   {"HO", ConstitutiveModelType::stIso_HO}, 
-  {"Holzapfel", ConstitutiveModelType::stIso_HO},
+  {"HolzapfelOgden", ConstitutiveModelType::stIso_HO},
+
+  {"HO_ma", ConstitutiveModelType::stIso_HO_ma}, 
+  {"HolzapfelOgden-ModifiedAnisotropy", ConstitutiveModelType::stIso_HO_ma},
 
   {"quad", ConstitutiveModelType::stVol_Quad},
   {"Quad", ConstitutiveModelType::stVol_Quad},

Code/Source/svFSI/distribute.cpp

@@ -1220,7 +1220,7 @@ void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c
    cm.bcast(cm_mod, lStM.Tf.gt.r, "lStM.Tf.gt.r");
    cm.bcast(cm_mod, lStM.Tf.gt.i, "lStM.Tf.gt.i");
   }
-
+  cm.bcast(cm_mod, &lStM.Tf.eta_s);
 }
 
 

Code/Source/svFSI/eq_assem.cpp

@@ -102,7 +102,7 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector<double>&
 
     for (int g = 0; g < lFa.nG; g++) {
       Vector<double> nV(nsd);
-      auto Nx = lFa.Nx.slice(g);
+      auto Nx = lFa.Nx.rslice(g);
       nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV);
       double Jac = sqrt(utils::norm(nV));
       nV = nV / Jac;
@@ -271,7 +271,7 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const
 
       // Get surface normal vector
       Vector<double> nV(nsd);
-      auto Nx_g = lFa.Nx.slice(g);
+      auto Nx_g = lFa.Nx.rslice(g);
       nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx_g, nV);
       Jac = sqrt(utils::norm(nV));
       nV = nV / Jac;
@@ -302,6 +302,74 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const
   }
 }
 
+/// @brief Update the surface integral involved in the coupled/resistance BC
+/// contribution to the stiffness matrix to reflect deformed geometry, if using
+/// a follower pressure load.
+/// The value of this integral is stored in lhs%face%val.
+/// This integral is sV = int_Gammat (Na * n_i) (See Brown et al. 2024, Eq. 56)
+/// where Na is the shape function and n_i is the normal vector.
+///
+/// This function updates the variable lhs%face%val with the new value, which
+/// is eventually used in ADDBCMUL() in the linear solver to add the contribution
+/// from the resistance BC to the matrix-vector product of the tangent matrix and
+/// an arbitrary vector.
+void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa)
+{
+  using namespace consts;
+  using namespace utils;
+  using namespace fsi_linear_solver;
+
+  #define n_debug_fsi_ls_upd
+  #ifdef debug_fsi_ls_upd
+  DebugMsg dmsg(__func__, com_mod.cm.idcm());
+  dmsg.banner();
+  dmsg << "lFa.name: " << lFa.name;
+  #endif
+
+  auto& cm = com_mod.cm;
+  int nsd = com_mod.nsd;
+  int tnNo = com_mod.tnNo;
+
+  int iM = lFa.iM;
+  int nNo = lFa.nNo;
+
+  Array<double> sVl(nsd,nNo); 
+  Array<double> sV(nsd,tnNo); 
+
+  // Updating the value of the surface integral of the normal vector
+  // using the deformed configuration ('n' = new = timestep n+1)
+  sV = 0.0;
+  for (int e = 0; e < lFa.nEl; e++) {
+    if (lFa.eType == ElementType::NRB) {
+      // CALL NRBNNXB(msh(iM),lFa,e)
+    }
+    for (int g = 0; g < lFa.nG; g++) {
+      Vector<double> n(nsd);
+      auto Nx = lFa.Nx.rslice(g);
+
+      auto cfg = MechanicalConfigurationType::new_timestep;
+
+      nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, cfg);
+      // 
+      for (int a = 0; a < lFa.eNoN; a++) {
+        int Ac = lFa.IEN(a,e);
+        for (int i = 0; i < nsd; i++) {
+          sV(i,Ac) = sV(i,Ac) + lFa.N(a,g)*lFa.w(g)*n(i);
+        }
+      }
+    }
+  }
+
+  if (sVl.size() != 0) { 
+    for (int a = 0; a < lFa.nNo; a++) {
+      int Ac = lFa.gN(a);
+      sVl.set_col(a, sV.col(Ac));
+    }
+  }
+  // Update lhs.face(i).val with the new value of the surface integral
+  fsils_bc_update(com_mod.lhs, lBc.lsPtr, lFa.nNo, nsd, sVl); 
+};
+
 /// @brief This routine assembles the equation on a given mesh.
 ///
 /// Ag(tDof,tnNo), Yg(tDof,tnNo), Dg(tDof,tnNo)

Code/Source/svFSI/mat_fun_carray.h

@@ -458,6 +458,15 @@ double norm(const Vector<double>& u, const double v[N])
   return norm;
 }
 
+template <size_t N>
+double norm(const Vector<double>& u, const Vector<double>& v)
+{
+  double norm = 0.0;
+  for (int i = 0; i < N; i++) {
+    norm += u[i]*v[i];
+  }
+  return norm;
+}
 
 template <size_t N>
 void mat_symm(const double A[N][N], double S[N][N])