Cardiac chamber model (#87)

Co-authored-by: kharold23 <[email protected]> Co-authored-by: Karthik Menon <[email protected]> Co-authored-by: Karthik Menon <[email protected]> Co-authored-by: Karthik Menon <[email protected]> Co-authored-by: Karthik Menon <[email protected]> Co-authored-by: Karthik Menon <[email protected]>
SimVascular · Dec 11, 2023 · 8a42c12 · 8a42c12
1 parent 75e2c88
commit 8a42c12
Show file tree

Hide file tree

Showing 23 changed files with 903 additions and 53 deletions.
diff --git a/docs/references.bib b/docs/references.bib
@@ -101,3 +101,13 @@ @article{sankaran2012patient
   year={2012},
   publisher={Springer}
 }
+
+@article{kerckhoffs2007coupling,
+  title={Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation},
+  author={Kerckhoffs, Roy CP and Neal, Maxwell L and Gu, Quan and Bassingthwaighte, James B and Omens, Jeff H and McCulloch, Andrew D},
+  journal={Annals of biomedical engineering},
+  volume={35},
+  pages={1--18},
+  year={2007},
+  publisher={Springer}
+}
diff --git a/src/algebra/SparseSystem.h b/src/algebra/SparseSystem.h
@@ -68,11 +68,11 @@ class Model;
  * Integrator. We then use the Newton-Raphson method to iteratively solve
  *
  * \f[
- * \mathbf{K}^{i} \cdot \Delta\mathbf{y}^{i} = - \mathbf{r}^{i},
+ * \mathbf{K}^{i} \cdot \Delta\dot{\mathbf{y}}^{i} = - \mathbf{r}^{i},
  * \f]
  *
- * with solution increment \f$\Delta\mathbf{y}\f$ in iteration \f$i\f$. The
- * linearization of the time-discretized system is
+ * with solution increment \f$\Delta\dot{\mathbf{y}}^{i}\f$ in iteration
+ * \f$i\f$. The linearization of the time-discretized system is
  *
  * \f[
  * \mathbf{K} =

diff --git a/src/model/BlockType.h b/src/model/BlockType.h
@@ -53,7 +53,8 @@ enum class BlockType {
   closed_loop_coronary_right_bc = 10,
   closed_loop_rcr_bc = 11,
   closed_loop_heart_pulmonary = 12,
-  valve_tanh = 13
+  valve_tanh = 13,
+  chamber_elastance_inductor = 14
 };
 
 /**
@@ -66,7 +67,8 @@ enum class BlockClass {
   boundary_condition = 2,
   closed_loop = 3,
   external = 4,
-  valve = 5
+  valve = 5,
+  chamber = 6
 };
 
 #endif
diff --git a/src/model/BloodVessel.h b/src/model/BloodVessel.h
@@ -61,18 +61,18 @@
  * ### Governing equations
  *
  * \f[
- * P_\text{in}^{e}-P_\text{out}^{e} - (R + S|Q_\text{in}|) Q_\text{in}^{e}-L
- * \dot{Q}_\text{out}^{e}=0 \f]
+ * P_\text{in}-P_\text{out} - (R + S|Q_\text{in}|) Q_\text{in}-L
+ * \dot{Q}_\text{out}=0 \f]
  *
  * \f[
- * Q_\text{in}^{e}-Q_\text{out}^{e} - C \dot{P}_\text{in}^{e}+C(R +
- * 2S|Q_\text{in}|) \dot{Q}_{in}^{e}=0 \f]
+ * Q_\text{in}-Q_\text{out} - C \dot{P}_\text{in}+C(R +
+ * 2S|Q_\text{in}|) \dot{Q}_{in}=0 \f]
  *
  * ### Local contributions
  *
  * \f[
- * \mathbf{y}^{e}=\left[\begin{array}{llll}P_{i n}^{e} & Q_{in}^{e} &
- * P_{out}^{e} & Q_{out}^{e}\end{array}\right]^\text{T} \f]
+ * \mathbf{y}^{e}=\left[\begin{array}{llll}P_{i n} & Q_{in} &
+ * P_{out} & Q_{out}\end{array}\right]^\text{T} \f]
  *
  * \f[
  * \mathbf{F}^{e}=\left[\begin{array}{cccc}

diff --git a/src/model/BloodVesselJunction.h b/src/model/BloodVesselJunction.h
@@ -88,9 +88,9 @@
  * ### Local contributions
  *
  * \f[
- * \mathbf{y}^{e}=\left[\begin{array}{lllllll}P_\text{in}^{e} & Q_\text{in}^{e}
- * & P_{out, 1}^{e} & Q_{out, 1}^{e} &
- * \dots & P_{out, i}^{e} & Q_{out, i}^{e}\end{array}\right] \f]
+ * \mathbf{y}^{e}=\left[\begin{array}{lllllll}P_\text{in} & Q_\text{in}
+ * & P_{out, 1} & Q_{out, 1} &
+ * \dots & P_{out, i} & Q_{out, i}\end{array}\right] \f]
  *
  * \f[
  * \mathbf{F}^{e} = \left[\begin{array}{ccccccc}

diff --git a/src/model/CMakeLists.txt b/src/model/CMakeLists.txt
@@ -36,6 +36,7 @@ set(CXXSRCS
   Block.cpp
   BloodVessel.cpp
   BloodVesselJunction.cpp
+  ChamberElastanceInductor.cpp
   ClosedLoopCoronaryBC.cpp
   ClosedLoopCoronaryLeftBC.cpp
   ClosedLoopCoronaryRightBC.cpp
@@ -60,6 +61,7 @@ set(HDRS
   BlockType.h
   BloodVessel.h
   BloodVesselJunction.h
+  ChamberElastanceInductor.h
   ClosedLoopCoronaryBC.h
   ClosedLoopCoronaryLeftBC.h
   ClosedLoopCoronaryRightBC.h

diff --git a/src/model/ChamberElastanceInductor.cpp b/src/model/ChamberElastanceInductor.cpp
@@ -0,0 +1,89 @@
+// Copyright (c) Stanford University, The Regents of the University of
+//               California, and others.
+//
+// All Rights Reserved.
+//
+// See Copyright-SimVascular.txt for additional details.
+//
+// Permission is hereby granted, free of charge, to any person obtaining
+// a copy of this software and associated documentation files (the
+// "Software"), to deal in the Software without restriction, including
+// without limitation the rights to use, copy, modify, merge, publish,
+// distribute, sublicense, and/or sell copies of the Software, and to
+// permit persons to whom the Software is furnished to do so, subject
+// to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included
+// in all copies or substantial portions of the Software.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+// OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+#include "ChamberElastanceInductor.h"
+
+void ChamberElastanceInductor::setup_dofs(DOFHandler &dofhandler) {
+  // Internal variable is chamber volume
+  Block::setup_dofs_(dofhandler, 3, {"Vc"});
+}
+
+void ChamberElastanceInductor::update_constant(
+    SparseSystem &system, std::vector<double> &parameters) {
+  double L = parameters[global_param_ids[ParamId::IMPEDANCE]];
+
+  // Eq 0: P_in - E(t)(Vc - Vrest) = 0
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0;
+
+  // Eq 1: P_in - P_out - L*dQ_out = 0
+  system.F.coeffRef(global_eqn_ids[1], global_var_ids[0]) = 1.0;
+  system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -1.0;
+  system.E.coeffRef(global_eqn_ids[1], global_var_ids[3]) = -L;
+
+  // Eq 2: Q_in - Q_out - dVc = 0
+  system.F.coeffRef(global_eqn_ids[2], global_var_ids[1]) = 1.0;
+  system.F.coeffRef(global_eqn_ids[2], global_var_ids[3]) = -1.0;
+  system.E.coeffRef(global_eqn_ids[2], global_var_ids[4]) = -1.0;
+}
+
+void ChamberElastanceInductor::update_time(SparseSystem &system,
+                                           std::vector<double> &parameters) {
+  get_elastance_values(parameters);
+
+  // Eq 0: P_in - E(t)(Vc - Vrest) = P_in - E(t)*Vc + E(t)*Vrest = 0
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[4]) = -1 * Elas;
+  system.C.coeffRef(global_eqn_ids[0]) = Elas * Vrest;
+}
+
+void ChamberElastanceInductor::get_elastance_values(
+    std::vector<double> &parameters) {
+  double Emax = parameters[global_param_ids[ParamId::EMAX]];
+  double Emin = parameters[global_param_ids[ParamId::EMIN]];
+  double Vrd = parameters[global_param_ids[ParamId::VRD]];
+  double Vrs = parameters[global_param_ids[ParamId::VRS]];
+  double t_active = parameters[global_param_ids[ParamId::TACTIVE]];
+  double t_twitch = parameters[global_param_ids[ParamId::TTWITCH]];
+
+  auto T_cardiac = model->cardiac_cycle_period;
+  auto t_in_cycle = fmod(model->time, T_cardiac);
+
+  double t_contract = 0;
+  if (t_in_cycle >= t_active) {
+    t_contract = t_in_cycle - t_active;
+  }
+
+  double act = 0;
+  if (t_contract <= t_twitch) {
+    act = -0.5 * cos(2 * M_PI * t_contract / t_twitch) + 0.5;
+  }
+
+  Vrest = (1.0 - act) * (Vrd - Vrs) + Vrs;
+  Elas = (Emax - Emin) * act + Emin;
+}