Different material models for ChamberSphere

applications/svZeroDGUI/static/index.html

@@ -4,9 +4,8 @@
       <title>Vasculature Simulation </title>
       <link rel="stylesheet" href="static/css/style.css">
       <meta name="viewport" content="width=device-width, user-scalable=no, initial-scale=1, maximum-scale=1">
-      <script src="https://unpkg.com/cytoscape/dist/cytoscape.min.js"></script>
+      <script src="/static/js/cytoscape.min.js"></script>
 
-    <!-- for testing with local version of cytoscape.js -->
       <script src="https://cdnjs.cloudflare.com/ajax/libs/lodash.js/4.17.10/lodash.js"></script>
       <script src="/static/js/cytoscape-edgehandles.js"></script>
       <script src="/static/js/app.js"></script>

scripts/ChamberSphere_expmat.yaml

@@ -0,0 +1,61 @@
+variables:
+  - Pin
+  - Qin
+  - Pout
+  - Qout
+  - radius
+  - velo
+  - stress
+  - tau
+  - volume
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dPout_dt
+  - dQout_dt
+  - dradius_dt
+  - dvelo_dt
+  - dstress_dt
+  - dtau_dt
+  - dvolume_dt
+
+constants:
+  - rho
+  - thick0
+  - radius0
+  - C0
+  - C1
+  - C2
+  - C3
+  - eta
+  - act
+  - act_plus
+  - sigma_max
+
+time_dependent:
+  - act
+  - act_plus
+
+helper_functions: |
+  def CG(radius):
+      return (1 + (radius / radius0)) ** 2
+
+  def dCG(radius, dradius_dt):
+      return 2 * (1 + (radius / radius0)) * (1 / radius0) * dradius_dt
+
+  def dW1(radius): 
+      return C0 * exp(C1 * (((1 + (radius / radius0))**(-4) + 2*(1 + (radius / radius0))**2) - 3)**2) * 2 * C1 * (((1 + (radius / radius0))**(-4) + 2*(1 + (radius / radius0))**2) - 3)
+
+  def dW4(radius):
+      return C2 * exp(C3 * ((((1 + (radius / radius0)) ** 2)-1)**2)) * 2 * C3 * (((1 + (radius / radius0)) ** 2) - 1)
+
+
+residuals:
+  - rho * thick0 * dvelo_dt + (thick0 / radius0) * (1 + (radius / radius0)) * stress - Pout * CG(radius)
+  - -stress + 4 * dW1(radius) *(1 - CG(radius)**(-3)) + 2 * dW4(radius)  + eta * dCG(radius,dradius_dt) * (1 + 2  * (CG(radius)**(-6))) + tau 
+  - 4 * pi * radius0 ** 2 * CG(radius) * velo - dvolume_dt
+  - dtau_dt + act * tau - sigma_max * act_plus
+  - dradius_dt - velo
+  - Qin - Qout - dvolume_dt
+  - Pin - Pout

scripts/jacobian.py

@@ -1,4 +1,4 @@
-from sympy import symbols, Matrix, simplify, pi, Abs
+from sympy import symbols, Matrix, simplify, pi, Abs, exp
 from sympy.printing import ccode
 import re
 import pdb
@@ -16,6 +16,7 @@ def load_model(filepath):
     context = {str(s): s for s in list(variables) + list(derivatives) + list(constants)}
     context['pi'] = pi
     context['abs'] = Abs
+    context['exp'] = exp
 
     if 'helper_functions' in data:
         exec(data['helper_functions'], context, context)

src/model/Block.h

@@ -15,6 +15,7 @@
 
 #include "ActivationFunction.h"
 #include "BlockType.h"
+#include "SphereMaterial.h"
 #include "DOFHandler.h"
 #include "Parameter.h"
 #include "SparseSystem.h"
@@ -297,6 +298,18 @@ class Block {
   virtual void set_activation_function(std::unique_ptr<ActivationFunction> af) {
     (void)af;  // Included to avoid unused parameter warning
   }
+
+  /**
+   * @brief Set wall material (for ChamberSphere blocks).
+   *
+   * Default no-op. Overridden by ChamberSphere to take ownership of the
+   * material.
+   *
+   * @param m Unique pointer to the material (caller transfers ownership)
+   */
+  virtual void set_material(std::unique_ptr<SphereMaterial> m) {
+    (void)m;  // Included to avoid unused parameter warning
+  }
 };
 
 #endif

src/model/CMakeLists.txt

@@ -5,8 +5,9 @@
 
 set(lib svzero_model_library)
 
-set(CXXSRCS 
+set(CXXSRCS
   ActivationFunction.cpp
+  SphereMaterial.cpp
   BloodVesselRC.cpp
   Block.cpp
   BloodVessel.cpp
@@ -37,8 +38,9 @@ set(CXXSRCS
   PiecewiseValve.cpp
 )
 
-set(HDRS 
+set(HDRS
   ActivationFunction.h
+  SphereMaterial.h
   BloodVesselRC.h
   Block.h
   BlockType.h

src/model/ChamberSphere.cpp

@@ -53,13 +53,10 @@ void ChamberSphere::update_solution(
     SparseSystem& system, std::vector<double>& parameters,
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
-  const double W1 = parameters[global_param_ids[ParamId::W1]];
-  const double W2 = parameters[global_param_ids[ParamId::W2]];
-  const double eta = parameters[global_param_ids[ParamId::eta]];
   const double thick0 = parameters[global_param_ids[ParamId::thick0]];
   const double sigma_max = parameters[global_param_ids[ParamId::sigma_max]];
-
   const double radius0 = parameters[global_param_ids[ParamId::radius0]];
+
   const double velo = y[global_var_ids[5]];
   const double Pout = y[global_var_ids[2]];
   const double stress = y[global_var_ids[6]];
@@ -77,23 +74,13 @@ void ChamberSphere::update_solution(
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[6]) =
       thick0 * (radius + radius0) / pow(radius0, 2);
 
-  // spherical stress
-  system.C.coeffRef(global_eqn_ids[1]) =
-      2 *
-      (dradius_dt * eta * (2 * pow(radius0, 12) + pow(radius + radius0, 12)) +
-       2 * pow(radius + radius0, 5) *
-           (-pow(radius0, 6) + pow(radius + radius0, 6)) *
-           (W1 * pow(radius0, 2) + W2 * pow(radius + radius0, 2))) /
-      (pow(radius0, 2) * pow(radius + radius0, 11));
+  // spherical stress (material-dependent)
+  const auto mat = material_->compute(radius, radius0, dradius_dt);
+  system.C.coeffRef(global_eqn_ids[1]) = mat.C_val;
   system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[4]) =
-      24 * W1 * pow(radius0, 6) / pow(radius + radius0, 7) +
-      8 * W2 * radius / pow(radius0, 2) +
-      16 * W2 * pow(radius0, 4) / pow(radius + radius0, 5) + 8 * W2 / radius0 -
-      44 * dradius_dt * eta * pow(radius0, 10) / pow(radius + radius0, 12) +
-      2 * dradius_dt * eta / pow(radius0, 2);
+      mat.dC_dy_radius;
   system.dC_dydot.coeffRef(global_eqn_ids[1], global_var_ids[4]) =
-      2 * eta * (2 * pow(radius0, 12) + pow(radius + radius0, 12)) /
-      (pow(radius0, 2) * pow(radius + radius0, 11));
+      mat.dC_dydot_radius;
 
   // volume change
   system.C.coeffRef(global_eqn_ids[2]) =
@@ -107,6 +94,10 @@ void ChamberSphere::update_solution(
   system.C.coeffRef(global_eqn_ids[3]) = -act_plus * sigma_max;
 }
 
+void ChamberSphere::set_material(std::unique_ptr<SphereMaterial> m) {
+  material_ = std::move(m);
+}
+
 void ChamberSphere::get_elastance_values(std::vector<double>& parameters) {
   const double alpha_max = parameters[global_param_ids[ParamId::alpha_max]];
   const double alpha_min = parameters[global_param_ids[ParamId::alpha_min]];

src/model/ChamberSphere.h

@@ -9,7 +9,10 @@
 
 #include <math.h>
 
+#include <memory>
+
 #include "Block.h"
+#include "SphereMaterial.h"
 #include "SparseSystem.h"
 
 /**
@@ -83,17 +86,14 @@
  * * `rho` - Density \f$\rho\f$
  * * `thick0` - Wall thickness \f$d_0\f$
  * * `radius0` - Reference radius \f$r_0\f$
- * * `W1` - Material constant \f$W_1\f$
- * * `W2` - Material constant \f$W_2\f$
- * * `eta` - Viscosity parameter \f$\eta\f$
  * * `sigma_max` - Maximum active stress \f$\sigma_\text{max}\f$
  * * `alpha_max` - Maximum activation parameter \f$\alpha_\text{max}\f$
  * * `alpha_min` - Minimum activation parameter \f$\alpha_\text{min}\f$
  * * `tsys` - Systole timing parameter \f$t_\text{sys}\f$
  * * `tdias` - Diastole timing parameter \f$t_\text{dias}\f$
  * * `steepness` - Activation steepness parameter \f$\gamma\f$
  *
- * ### Usage in json configuration file
+ * ### Usage in json configuration file (with example material)
  *
  *     "vessels": [
  *        {
@@ -106,16 +106,19 @@
  *                "rho" : 1e3,
  *                "thick0" : 0.01,
  *                "radius0" : 0.05,
- *                "W1" : 10e3,
- *                "W2" : 40,
- *                "eta" : 10.0,
  *                "sigma_max" : 185e3,
  *                "alpha_max": 30.0,
  *                "alpha_min": -30.0,
  *                "tsys": 0.170,
  *                "tdias": 0.484,
  *                "steepness": 0.005
  *            }
+ *            "material": {
+ *                "type": "mooney_rivlin",
+ *                "W1": 10e3,
+ *                "W2": 40,
+ *                "eta": 10.0
+ *            }
  *        }
  *     ]
  *
@@ -140,15 +143,12 @@ class ChamberSphere : public Block {
     rho = 0,
     thick0 = 1,
     radius0 = 2,
-    W1 = 3,
-    W2 = 4,
-    eta = 5,
-    sigma_max = 6,
-    alpha_max = 7,
-    alpha_min = 8,
-    tsys = 9,
-    tdias = 10,
-    steepness = 11
+    sigma_max = 3,
+    alpha_max = 4,
+    alpha_min = 5,
+    tsys = 6,
+    tdias = 7,
+    steepness = 8
   };
 
   /**
@@ -162,9 +162,6 @@ class ChamberSphere : public Block {
               {{"rho", InputParameter()},
                {"thick0", InputParameter()},
                {"radius0", InputParameter()},
-               {"W1", InputParameter()},
-               {"W2", InputParameter()},
-               {"eta", InputParameter()},
                {"sigma_max", InputParameter()},
                {"alpha_max", InputParameter()},
                {"alpha_min", InputParameter()},
@@ -225,10 +222,14 @@ class ChamberSphere : public Block {
    */
   void get_elastance_values(std::vector<double>& parameters);
 
+  void set_material(std::unique_ptr<SphereMaterial> m) override;
+
  private:
   double act = 0.0;       // activation function
   double act_plus = 0.0;  // act_plus = max(act, 0)
 
+  std::unique_ptr<SphereMaterial> material_;
+
   /**
    * @brief Number of triplets of element
    *