coup_wss as material parameter

Code/Source/svFSI/ComMod.h

@@ -389,6 +389,7 @@ class grModelType
  int n_t_pre = 0;
  int n_t_end = 0;
  int example = 0;
+ bool coup_wss = false;
  double KsKi = 0.0;
  double curve = 0.0;
  double mult = 0.0;
@@ -1523,9 +1524,6 @@ class ComMod {
  /// @brief Number of internal growth and remodeling variables
  int nGrInt = 0;
 
- /// @brief Is growth and remodeling coupled to fluids?
- bool gr_coup_wss = 0;
-
  //----- double members -----//
 
  /// @brief Time step size

Code/Source/svFSI/Parameters.cpp

@@ -736,6 +736,7 @@ GREquilibratedParameters::GREquilibratedParameters()
  set_parameter("n_t_pre", 0, required, n_t_pre);
  set_parameter("n_t_end", 0, required, n_t_end);
  set_parameter("example", 0, required, example);
+ set_parameter("coup_wss", false, required, coup_wss);
  set_parameter("KsKi", 0.0, required, KsKi);
  set_parameter("curve", 0.0, required, curve);
  set_parameter("mult", 0.0, required, mult);

Code/Source/svFSI/Parameters.h

@@ -514,6 +514,7 @@ class GREquilibratedParameters : public ParameterLists
  Parameter<int> n_t_pre;
  Parameter<int> n_t_end;
  Parameter<int> example;
+ Parameter<bool> coup_wss;
  Parameter<double> KsKi;
  Parameter<double> curve;
  Parameter<double> mult;

Code/Source/svFSI/distribute.cpp

@@ -1220,6 +1220,7 @@ void dist_gr_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm
 {
  using namespace consts;
 
+ cm.bcast(cm_mod, &grM.coup_wss);
  cm.bcast(cm_mod, &grM.KsKi);
  cm.bcast(cm_mod, &grM.curve);
  cm.bcast(cm_mod, &grM.mult);

Code/Source/svFSI/gr_equilibrated.cpp

@@ -38,8 +38,8 @@ namespace gr_equilibrated_ns {
 void stress_tangent_(const grModelType &grM, const double Fe[3][3],
  const double time, const Vector<double> &eVWP,
  Vector<double> &grInt, double S_out[3][3],
- double CC_out[3][3][3][3], const bool coup_wss,
- const bool eval_s, const bool eval_cc) {
+ double CC_out[3][3][3][3], const bool eval_s,
+ const bool eval_cc) {
  // convert deformation gradient to FEBio format
  mat3d F(Fe[0][0], Fe[0][1], Fe[0][2], Fe[1][0], Fe[1][1], Fe[1][2], Fe[2][0],
  Fe[2][1], Fe[2][2]);
@@ -74,7 +74,7 @@ void stress_tangent_(const grModelType &grM, const double Fe[3][3],
 
  // get current time
  double t;
- if (coup_wss)
+ if (grM.coup_wss)
  // time from partitioned coupling (passed through input file)
  t = eVWP(7);
  else
@@ -614,7 +614,7 @@ void stress_tangent_(const grModelType &grM, const double Fe[3][3],
  ro / rIo * lt -
  (ro - rIo) / rIo * lr; // rIrIo -> rIorIo = 1 for F -> Fo
 
- if (coup_wss)
+ if (grM.coup_wss)
  tau_ratio = tau / tauo;
  else
  tau_ratio = pow(rIrIo, -3);
@@ -785,7 +785,7 @@ void stress_tangent_(const grModelType &grM, const double Fe[3][3],
  (IxIss - 2.0 * IoIss);
 
  // wss linearization
- if (!coup_wss) {
+ if (!grM.coup_wss) {
  css += 3.0 * pow(rIrIo, -4) * phim * 2.0 * Cratio * CS * EPS *
  exp(-Cratio * Cratio) / (1.0 - exp(-CB * CB)) *
  (ro / rIo / lt * saoxntt - (ro - rIo) / rIo / lr * saoxnrr);
@@ -967,4 +967,58 @@ void stress_tangent_(const grModelType &grM, const double Fe[3][3],
  grInt(k + 11) = phic_gp;
  }
 }
+
+void stress_tangent_stvk(const grModelType &grM, const double Fe[3][3],
+ const double time, const Vector<double> &eVWP,
+ Vector<double> &grInt, double S_out[3][3],
+ double CC_out[3][3][3][3], const bool coup_wss,
+ const bool eval_s, const bool eval_cc) {
+ // convert deformation gradient to FEBio format
+ mat3d F(Fe[0][0], Fe[0][1], Fe[0][2], Fe[1][0], Fe[1][1], Fe[1][2], Fe[2][0],
+ Fe[2][1], Fe[2][2]);
+
+ // material parameters
+ const double young = 240.56596E6;
+ const double nu = 0.4;
+
+ // lame parameters
+ const double mu = young / (2.0 * (1.0 + nu));
+ const double lambda = nu * young / ((1.0 + nu) * (1.0 - 2.0 * nu));
+
+ // define identity tensor and some useful dyadic products of the identity
+ // tensor
+ const mat3dd I(1.0);
+ tens4ds IxI = dyad1s(I);
+ tens4ds IoI = dyad4s(I);
+ tens4dmm cIxI = tens4dmm(IxI);
+ tens4dmm cIoI = tens4dmm(IoI);
+
+ // green-lagrange strain
+ mat3ds E = 0.5 * ((F.transpose() * F).sym() - I);
+
+ // stress
+ mat3ds S = lambda * E.tr() * I + 2.0 * mu * E;
+
+ // tangent
+ tens4dmm css = lambda * cIxI + 2.0 * mu * cIoI;
+
+ // convert to vector for FORTRAN
+ typedef double(*ten2)[3];
+ typedef double(*ten4)[3][3][3];
+
+ if (eval_s) {
+ for (int i = 0; i < 3; i++)
+ for (int j = 0; j < 3; j++)
+ S_out[i][j] = S(i, j);
+ }
+
+ if (eval_cc) {
+ for (int i = 0; i < 3; i++)
+ for (int j = 0; j < 3; j++)
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ CC_out[i][j][k][l] = css(i, j, k, l);
+ }
+}
+
 } // namespace gr_equilibrated_ns

Code/Source/svFSI/gr_equilibrated.h

@@ -47,20 +47,26 @@
 #include <vec3d.h>
 
 #include "Array.h"
-#include "Tensor4.h"
 #include "ComMod.h"
+#include "Tensor4.h"
 
 #ifndef GR_EQUILIBRATED
 #define GR_EQUILIBRATED
 
 namespace gr_equilibrated_ns {
 
-void stress_tangent_(const grModelType &grM, const double Fe[3][3], const double time,
- const Vector<double> &eVWP, Vector<double> &grInt,
- double S_out[3][3], double CC_out[3][3][3][3],
- const bool coup_wss, const bool eval_s = true,
+void stress_tangent_(const grModelType &grM, const double Fe[3][3],
+ const double time, const Vector<double> &eVWP,
+ Vector<double> &grInt, double S_out[3][3],
+ double CC_out[3][3][3][3], const bool eval_s = true,
  const bool eval_cc = true);
 
-};
+void stress_tangent_stvk(const grModelType &grM, const double Fe[3][3],
+ const double time, const Vector<double> &eVWP,
+ Vector<double> &grInt, double S_out[3][3],
+ double CC_out[3][3][3][3], const bool eval_s = true,
+ const bool eval_cc = true);
+
+}; // namespace gr_equilibrated_ns
 
 #endif

Code/Source/svFSI/gr_struct.h

@@ -122,14 +122,13 @@ void get_pk2cc(const ComMod &com_mod, const dmnType &lDmn, const double F[N][N],
  const auto &stM = lDmn.stM;
  const auto &grM = lDmn.grM;
  const int nsd = com_mod.nsd;
- const bool coup_wss = com_mod.gr_coup_wss;
 
  double CC[N][N][N][N];
 
  switch (stM.isoType) {
  case ConstitutiveModelType::GR_equi: {
  gr_equilibrated_ns::stress_tangent_(grM, F, com_mod.time, gr_props, gr_int,
- S, CC, coup_wss, eval_s, eval_cc);
+ S, CC, eval_s, eval_cc);
  } break;
 
  default:

Code/Source/svFSI/main.cpp

@@ -174,7 +174,6 @@ void iterate_solution(Simulation* simulation)
 
  bool exit_now = false;
  double elapsed_time = 0.0;
- com_mod.gr_coup_wss = false;
 
  // Uncomment these two lines to enable writting values to a file.
  //Array<double>::write_enabled = true;

Code/Source/svFSI/pic.cpp

@@ -652,11 +652,18 @@ void picp(Simulation* simulation)
  // Update the solution vectors? (new load step in FSG coupling)
  auto & lM = com_mod.msh[0];
  const bool update = std::fabs(lM.gr_props(12, 0) - 1.0) < 1.0e-6;
+
+ // Check if coupled to fluid
+ int cEq = com_mod.cEq;
+ auto &eq = com_mod.eq[cEq];
+ int cDmn = com_mod.cDmn;
+ auto &dmn = eq.dmn[cDmn];
+ const bool coup_wss = dmn.grM.coup_wss;
 
  // Predict new load step
  for (int i = s; i <= e; i++) {
  for (int j = 0; j < Do.ncols(); j++) {
- if (com_mod.gr_coup_wss && !update) {
+ if (coup_wss && !update) {
  Dn(i,j) = Do(i,j);
  } else {
  An(i,j) = Yo(i,j);

Code/Source/svFSI/set_equation_props.h

@@ -311,7 +311,7 @@ SetEquationPropertiesMapType set_equation_props = {
 
  // Set the possible equations for fsi: fluid (required), struct/ustruct/lElas
  EquationPhys phys { EquationType::phys_fluid, EquationType::phys_struct, EquationType::phys_ustruct, 
- EquationType::phys_lElas, EquationType::phys_gr };
+ EquationType::phys_lElas };
 
  // Set fluid properties.
  int n = 0;

Code/Source/svFSI/set_material_props.h

@@ -183,6 +183,7 @@ SeMaterialPropertiesMapType set_material_props = {
  lDmn.grM.n_t_pre = params.n_t_pre.value();
  lDmn.grM.n_t_end = params.n_t_end.value();
  lDmn.grM.example = params.example.value();
+ lDmn.grM.coup_wss = params.coup_wss.value();
  lDmn.grM.KsKi = params.KsKi.value();
  lDmn.grM.curve = params.curve.value();
  lDmn.grM.mult = params.mult.value();

tests/cases/gr/equilibrated/svFSI.xml

@@ -58,6 +58,7 @@
  <example> 0 </example>
  <n_t_pre> 1 </n_t_pre>
  <n_t_end> 4 </n_t_end>
+ <coup_wss> false </coup_wss>
  <KsKi> 1.0 </KsKi>
  <mult> 1.0 </mult>
  <rIo> 0.6468 </rIo>