Piecewise Chamber and Valve

docs/references.bib

@@ -125,3 +125,18 @@ @article{caruel13
  volume = {13},
  year = {2013}
 }
+
+@article{Regazzoni2022,
+title = {A cardiac electromechanical model coupled with a lumped-parameter model for closed-loop blood circulation},
+journal = {Journal of Computational Physics},
+volume = {457},
+pages = {111083},
+year = {2022},
+issn = {0021-9991},
+doi = {https://doi.org/10.1016/j.jcp.2022.111083},
+url = {https://www.sciencedirect.com/science/article/pii/S0021999122001450},
+author = {F. Regazzoni and M. Salvador and P.C. Africa and M. Fedele and L. Dedè and A. Quarteroni},
+keywords = {Mathematical modeling, Multiphysics models, Multiscale models, Cardiac modeling, Cardiac electromechanics},
+abstract = {We propose a novel mathematical and numerical model for cardiac electromechanics, wherein biophysically detailed core models describe the different physical processes concurring to the cardiac function. The core models, once suitably approximated, are coupled by a computationally efficient strategy. Our model is based on: (1) the combination of implicit-explicit (IMEX) schemes to solve the different core cardiac models, (2) an Artificial Neural Network based model, that surrogates a biophysically detailed but computationally demanding microscale model of active force generation and (3) appropriate partitioned schemes to couple the different models in this multiphysics setting. We employ a flexible and scalable intergrid transfer operator, which allows to interpolate Finite Element functions between nested meshes and, possibly, among arbitrary Finite Element spaces for the different core models. Our core 3D electromechanical model of the left ventricle is coupled with a closed-loop 0D model of the vascular network (and the other cardiac chambers) by an approach that is energy preserving. More precisely, we derive a balance law for the mechanical energy of the whole circulatory network. This provides a quantitative insight into the energy utilization, dissipation and transfer among the different compartments of the cardiovascular network and during different stages of the heartbeat. On this ground, a new tool is proposed to validate some energy indicators adopted in the daily clinical practice. A further contribution of this paper is the proposition of a robust algorithm for the reconstruction of the stress-free reference configuration. This feature is fundamental to correctly initialize our electromechanical simulations. As a matter of fact, the geometry acquired from medical imaging typically refers to a configuration affected by residual internal stresses, whereas the elastodynamics equations that govern the mechanics core model are related to a stress-free configuration. To prove the biophysical accuracy of our computational model, we address different scenarios of clinical interest, namely by varying preload, afterload and contractility.}
+}
+

src/model/BlockType.h

@@ -28,7 +28,9 @@ enum class BlockType {
   closed_loop_heart_pulmonary = 12,
   valve_tanh = 13,
   chamber_elastance_inductor = 14,
-  chamber_sphere = 15
+  chamber_sphere = 15,
+  piecewise_cosine_chamber = 17,
+  piecewise_valve = 18
 };
 
 /**

src/model/CMakeLists.txt

@@ -28,6 +28,8 @@ set(CXXSRCS
   ResistiveJunction.cpp
   ValveTanh.cpp
   WindkesselBC.cpp
+  PiecewiseCosineChamber.cpp
+  PiecewiseValve.cpp
 )
 
 set(HDRS 
@@ -54,6 +56,8 @@ set(HDRS
   ResistiveJunction.h
   ValveTanh.h
   WindkesselBC.h
+  PiecewiseCosineChamber.h
+  PiecewiseValve.h
 )
 
 add_library(${lib} OBJECT ${CXXSRCS} ${HDRS})

src/model/Model.cpp

@@ -28,7 +28,9 @@ Model::Model() {
       {"RESISTANCE", block_factory<ResistanceBC>()},
       {"resistive_junction", block_factory<ResistiveJunction>()},
       {"ValveTanh", block_factory<ValveTanh>()},
-      {"ChamberElastanceInductor", block_factory<ChamberElastanceInductor>()}};
+      {"ChamberElastanceInductor", block_factory<ChamberElastanceInductor>()},
+      {"PiecewiseValve", block_factory<PiecewiseValve>()},
+      {"PiecewiseCosineChamber", block_factory<PiecewiseCosineChamber>()}};
 }
 
 Model::~Model() {}
@@ -116,7 +118,7 @@ std::string Model::get_block_name(int block_id) const {
 int Model::add_node(const std::vector<Block*>& inlet_eles,
                     const std::vector<Block*>& outlet_eles,
                     const std::string_view& name) {
-  // DEBUG_MSG("Adding node " << name);
+  DEBUG_MSG("Adding node " << name);
   auto node = std::shared_ptr<Node>(
       new Node(node_count, inlet_eles, outlet_eles, this));
   nodes.push_back(node);
@@ -174,6 +176,10 @@ void Model::finalize() {
   for (auto& block : blocks) {
     block->setup_model_dependent_params();
   }
+
+  if (cardiac_cycle_period < 0.0) {
+    cardiac_cycle_period = 1.0;
+  }
 }
 
 int Model::get_num_blocks(bool internal) const {

src/model/Model.h

@@ -31,6 +31,8 @@
 #include "Node.h"
 #include "OpenLoopCoronaryBC.h"
 #include "Parameter.h"
+#include "PiecewiseCosineChamber.h"
+#include "PiecewiseValve.h"
 #include "PressureReferenceBC.h"
 #include "ResistanceBC.h"
 #include "ResistiveJunction.h"

src/model/OpenLoopCoronaryBC.h

@@ -73,7 +73,7 @@
  * Parameter sequence for constructing this block
  *
  * * `0` Ra: Small artery resistance
- * * `1` Ram: Microvascualar resistance
+ * * `1` Ram: Microvascualr resistance
  * * `2` Rv: Venous resistance
  * * `3` Ca: Small artery capacitance
  * * `4` Cim: Intramyocardial capacitance

src/model/PiecewiseCosineChamber.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 "PiecewiseCosineChamber.h"
+
+void PiecewiseCosineChamber::setup_dofs(DOFHandler& dofhandler) {
+  // Internal variable is chamber volume
+  Block::setup_dofs_(dofhandler, 3, {"Vc"});
+}
+
+void PiecewiseCosineChamber::update_constant(SparseSystem& system,
+                                             std::vector<double>& parameters) {
+  // 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 = 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;
+
+  // 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 PiecewiseCosineChamber::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 * parameters[global_param_ids[ParamId::VREST]];
+}
+
+void PiecewiseCosineChamber::get_elastance_values(
+    std::vector<double>& parameters) {
+  double Emax = parameters[global_param_ids[ParamId::EMAX]];
+  double Epass = parameters[global_param_ids[ParamId::EPASS]];
+  double Vrest = parameters[global_param_ids[ParamId::VREST]];
+  double contract_start = parameters[global_param_ids[ParamId::CSTART]];
+  double relax_start = parameters[global_param_ids[ParamId::RSTART]];
+  double contract_duration = parameters[global_param_ids[ParamId::CDUR]];
+  double relax_duration = parameters[global_param_ids[ParamId::RDUR]];
+
+  auto T_HB = model->cardiac_cycle_period;
+
+  double phi = 0;
+
+  auto piecewise_condition = fmod(model->time - contract_start, T_HB);
+
+  if (0 <= piecewise_condition && piecewise_condition < contract_duration) {
+    phi = 0.5 * (1 - cos((M_PI * piecewise_condition) / contract_duration));
+  } else {
+    piecewise_condition = fmod(model->time - relax_start, T_HB);
+    if (0 <= piecewise_condition && piecewise_condition < relax_duration) {
+      phi = 0.5 * (1 + cos((M_PI * piecewise_condition) / relax_duration));
+    }
+  }
+
+  Elas = Epass + Emax * phi;
+}