Add simple blocks

src/algebra/SparseSystem.cpp

@@ -101,6 +101,7 @@ void SparseSystem::update_jacobian(double time_coeff_ydot,
 void SparseSystem::solve() {
   solver->factorize(jacobian);
   if (solver->info() != Eigen::Success) {
+    std::cout<<"Jacobian: "<<std::endl<<jacobian<<std::endl;
     throw std::runtime_error(
         "System is singular. Check your model (connections, boundary "
         "conditions, parameters).");

src/model/Block.h

@@ -50,15 +50,15 @@
  * to the global system.
  */
 struct TripletsContributions {
-  TripletsContributions(){};
+  TripletsContributions() {};
   /**
    * @brief Set the number of triplets that the element contributes
    * to the global system.
    * @param F Contributions to F matrix
    * @param E Contributions to E matrix
    * @param D Contributions to dC/dy matrix
    */
-  TripletsContributions(int F, int E, int D) : F(F), E(E), D(D){};
+  TripletsContributions(int F, int E, int D) : F(F), E(E), D(D) {};
   /**
    * @brief Set the number of triplets that the element contributes
    * to the global system.
@@ -113,6 +113,7 @@ class Block {
 
   bool steady = false;             ///< Toggle steady behavior
   bool input_params_list = false;  ///< Are input parameters given as a list?
+  bool is_vessel = false;          ///< Is this block a blood vessel?
 
   /**
    * @brief Construct a new Block object

src/model/BlockType.h

@@ -54,7 +54,14 @@ enum class BlockType {
   closed_loop_rcr_bc = 11,
   closed_loop_heart_pulmonary = 12,
   valve_tanh = 13,
-  chamber_elastance_inductor = 14
+  chamber_elastance_inductor = 14,
+  resistance = 15,
+  capacitance = 16,
+  inductance = 17,
+  parallel_rc = 18,
+  parallel_rl = 19,
+  serial_rc = 20,
+  four_element_windkessel_bc = 21
 };
 
 /**

src/model/BloodVessel.h

@@ -168,7 +168,9 @@ class BloodVessel : public Block {
               {{"R_poiseuille", InputParameter()},
                {"C", InputParameter(true)},
                {"L", InputParameter(true)},
-               {"stenosis_coefficient", InputParameter(true)}}) {}
+               {"stenosis_coefficient", InputParameter(true)}}) {
+    is_vessel = true;
+  }
 
   /**
    * @brief Set up the degrees of freedom (DOF) of the block

src/model/CMakeLists.txt

@@ -36,6 +36,7 @@ set(CXXSRCS
   Block.cpp
   BloodVessel.cpp
   BloodVesselJunction.cpp
+  Capacitance.cpp
   ChamberElastanceInductor.cpp
   ClosedLoopCoronaryBC.cpp
   ClosedLoopCoronaryLeftBC.cpp
@@ -44,14 +45,20 @@ set(CXXSRCS
   ClosedLoopRCRBC.cpp
   DOFHandler.cpp
   FlowReferenceBC.cpp
+  FourElementWindkesselBC.cpp
+  Inductance.cpp
   Junction.cpp
   Model.cpp
   Node.cpp
   OpenLoopCoronaryBC.cpp
   Parameter.cpp
+  ParallelRC.cpp
+  ParallelRL.cpp
   PressureReferenceBC.cpp
+  Resistance.cpp
   ResistanceBC.cpp
   ResistiveJunction.cpp
+  SerialRC.cpp
   ValveTanh.cpp
   WindkesselBC.cpp
 )
@@ -61,6 +68,7 @@ set(HDRS
   BlockType.h
   BloodVessel.h
   BloodVesselJunction.h
+  Capacitance.h
   ChamberElastanceInductor.h
   ClosedLoopCoronaryBC.h
   ClosedLoopCoronaryLeftBC.h
@@ -69,14 +77,20 @@ set(HDRS
   ClosedLoopRCRBC.h
   DOFHandler.h
   FlowReferenceBC.h
+  FourElementWindkesselBC.h
+  Inductance.h
   Junction.h
   Model.h
   Node.h
   OpenLoopCoronaryBC.h
   Parameter.h
+  ParallelRC.h
+  ParallelRL.h
   PressureReferenceBC.h
+  Resistance.h
   ResistanceBC.h
   ResistiveJunction.h
+  SerialRC.h
   ValveTanh.h
   WindkesselBC.h
 )

src/model/Capacitance.cpp

@@ -0,0 +1,51 @@
+// 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 "Capacitance.h"
+
+void Capacitance::setup_dofs(DOFHandler &dofhandler) {
+  Block::setup_dofs_(dofhandler, 2, {});
+}
+
+void Capacitance::update_constant(SparseSystem &system,
+                                  std::vector<double> &parameters) {
+  double capacitance = parameters[global_param_ids[ParamId::CAPACITANCE]];
+
+  // unknowns: y = [P_in, Q_in, P_out, Q_out]
+
+  // Q_{in} - Q_{out} - C * \dot{P}_{in} = 0
+  system.E.coeffRef(global_eqn_ids[0], global_var_ids[0]) = -capacitance;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = 1.0;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[3]) = -1.0;
+
+  // 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;
+}

src/model/Capacitance.h

@@ -0,0 +1,150 @@
+// 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.
+/**
+ * @file Capacitance.h
+ * @brief model::Capacitance source file
+ */
+#ifndef SVZERODSOLVER_MODEL_CAPACITANCE_HPP_
+#define SVZERODSOLVER_MODEL_CAPACITANCE_HPP_
+
+#include <math.h>
+
+#include "Block.h"
+#include "SparseSystem.h"
+
+/**
+ * @brief Simple capacitance element for 0D blood flow modeling
+ *
+ * Models a pure capacitance element relating flow to change in pressure.
+ *
+ * \f[
+ * \begin{circuitikz} \draw
+ *   node[left] {$Q_{in}$} [-latex] (0,0) -- (0.8,0);
+ *   (1,0) node[anchor=south]{$P_{in}$}
+ *   to [short, -*] (1,0)
+ *   to [C, l=$C$] (5,0)
+ *   to [short, -*] (7,0)
+ *   node[anchor=south]{$P_{out}$}
+ *   (7.2,0) -- (8,0) node[right] {$Q_{out}$};
+ * \end{circuitikz}
+ * \f]
+ *
+ * ### Governing equations
+ *
+ * \f[
+ * P_\text{in}-P_\text{out} = 0 \f]
+ *
+ * \f[
+ * Q_\text{in} - C (\dot{P}_\text{in} - \dot{P}_\text{out}) = 0 \f]
+ *
+ * ### Local contributions
+ *
+ * \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}
+ * 0 & 1 &  0 & 0 \\
+ * 1 & 0 & -1 & 0
+ * \end{array}\right]
+ * \f]
+ *
+ * \f[
+ * \mathbf{E}^{e}=\left[\begin{array}{cccc}
+ * -C &  0 & C & 0 \\
+ *  0 &  0 & 0 & 0
+ * \end{array}\right]
+ * \f]
+ *
+ * ### Parameters
+ *
+ * Parameter sequence for constructing this block
+ *
+ * * `0` Capacitance value
+ *
+ * ### Internal variables
+ *
+ * This block has no internal variables.
+ *
+ */
+class Capacitance : public Block {
+ public:
+  /**
+   * @brief Local IDs of the parameters
+   *
+   */
+  enum ParamId {
+    CAPACITANCE = 0,
+  };
+
+  /**
+   * @brief Construct a new Capacitance object
+   *
+   * @param id Global ID of the block
+   * @param model The model to which the block belongs
+   */
+  Capacitance(int id, Model *model)
+      : Block(id, model, BlockType::capacitance, BlockClass::vessel,
+              {{"C", InputParameter()}}) {
+    is_vessel = true;
+  }
+
+  /**
+   * @brief Set up the degrees of freedom (DOF) of the block
+   *
+   * Set \ref global_var_ids and \ref global_eqn_ids of the element based on the
+   * number of equations and the number of internal variables of the
+   * element.
+   *
+   * @param dofhandler Degree-of-freedom handler to register variables and
+   * equations at
+   */
+  void setup_dofs(DOFHandler &dofhandler);
+
+  /**
+   * @brief Update the constant contributions of the element in a sparse
+   * system
+   *
+   * @param system System to update contributions at
+   * @param parameters Parameters of the model
+   */
+  void update_constant(SparseSystem &system, std::vector<double> &parameters);
+
+  /**
+   * @brief Number of triplets of element
+   *
+   * Number of triplets that the element contributes to the global system
+   * (relevant for sparse memory reservation)
+   */
+  TripletsContributions num_triplets{4, 1, 0};
+};
+
+#endif  // SVZERODSOLVER_MODEL_CAPACITANCE_HPP_

src/model/FourElementWindkesselBC.cpp

@@ -0,0 +1,69 @@
+// 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 "FourElementWindkesselBC.h"
+
+void FourElementWindkesselBC::setup_dofs(DOFHandler &dofhandler) {
+  Block::setup_dofs_(dofhandler, 3, {"pressure_c", "pressure_L"});
+}
+
+void FourElementWindkesselBC::update_constant(SparseSystem &system,
+                                            std::vector<double> &parameters) {
+  // P_in - P_L = 0 (constant part)
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0;
+  system.F.coeffRef(global_eqn_ids[0], global_var_ids[3]) = -1.0;
+
+  // P_L - P_c - R_p*Q_in = 0 (constant part)
+  system.F.coeffRef(global_eqn_ids[1], global_var_ids[2]) = -1.0;
+  system.F.coeffRef(global_eqn_ids[1], global_var_ids[3]) = 1.0;
+
+  // R_d*Q_in - P_c + P_ref = 0 (constant part)
+  system.F.coeffRef(global_eqn_ids[2], global_var_ids[2]) = -1.0;
+}
+
+void FourElementWindkesselBC::update_time(SparseSystem &system,
+                                        std::vector<double> &parameters) {
+  double inductance = parameters[global_param_ids[0]];
+  double resistance_proximal = parameters[global_param_ids[1]];
+  double capacitance = parameters[global_param_ids[2]];
+  double resistance_distal = parameters[global_param_ids[3]];
+  double pressure_distal = parameters[global_param_ids[4]];
+
+  // P_in - P_L - L*dQ_in/dt = 0 (time-dependent part)
+  system.E.coeffRef(global_eqn_ids[0], global_var_ids[1]) = -inductance;
+
+  // P_L - P_c - R_p*Q_in = 0 (time-dependent part)
+  system.F.coeffRef(global_eqn_ids[1], global_var_ids[1]) = -resistance_proximal;
+
+  // R_d*Q_in - P_c + P_ref - R_d*C*dP_c/dt = 0 (time-dependent part)
+  system.F.coeffRef(global_eqn_ids[2], global_var_ids[1]) = resistance_distal;
+  system.E.coeffRef(global_eqn_ids[2], global_var_ids[2]) = -resistance_distal * capacitance;
+  system.C(global_eqn_ids[2]) = pressure_distal;
+}