Add time-dependent function input for boundary conditions

CMakeLists.txt

@@ -133,6 +133,25 @@ FetchContent_Declare(
 
 FetchContent_MakeAvailable(pybind11_json)
 
+# -----------------------------------------------------------------------------
+# Fetch ExprTk
+# -----------------------------------------------------------------------------
+# ExprTk is a header-only mathematical expression parser used by Parameter.
+#
+FetchContent_Declare(
+    exprtk
+    GIT_REPOSITORY https://github.com/ArashPartow/exprtk.git
+    GIT_TAG        master
+    GIT_SHALLOW    TRUE
+    GIT_PROGRESS   TRUE
+)
+
+FetchContent_MakeAvailable(exprtk)
+
+add_library(exprtk INTERFACE)
+add_library(exprtk::exprtk ALIAS exprtk)
+target_include_directories(exprtk INTERFACE ${exprtk_SOURCE_DIR})
+
 # -----------------------------------------------------------------------------
 # Set executables for each application.
 # -----------------------------------------------------------------------------

src/model/Block.h

@@ -25,15 +25,15 @@
  * to the global system.
  */
 struct TripletsContributions {
-  TripletsContributions() {};
+  TripletsContributions() = default;
   /**
    * @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.

src/model/CMakeLists.txt

@@ -72,6 +72,11 @@ set(HDRS
 
 add_library(${lib} OBJECT ${CXXSRCS} ${HDRS})
 
+# exprtk.hpp generates too many COFF sections for the default MSVC object format
+if(MSVC)
+  set_source_files_properties(Parameter.cpp PROPERTIES COMPILE_FLAGS "/bigobj")
+endif()
+
 target_include_directories( ${lib} PUBLIC
   ${CMAKE_SOURCE_DIR}/src/algebra
   ${CMAKE_SOURCE_DIR}/src/model
@@ -80,4 +85,4 @@ target_include_directories( ${lib} PUBLIC
 
 target_link_libraries( ${lib} Eigen3::Eigen )
 target_link_libraries( ${lib} nlohmann_json::nlohmann_json )
-
+target_link_libraries( ${lib} exprtk::exprtk )

src/model/FlowReferenceBC.h

@@ -84,8 +84,9 @@ class FlowReferenceBC : public Block {
    */
   FlowReferenceBC(int id, Model* model)
       : Block(id, model, BlockType::flow_bc, BlockClass::boundary_condition,
-              {{"t", InputParameter(false, true)},
-               {"Q", InputParameter(false, true)}}) {}
+              {{"t", InputParameter(true, true)},
+               {"Q", InputParameter(true, true)},
+               {"fn", InputParameter(true, false, false, true)}}) {}
 
   /**
    * @brief Set up the degrees of freedom (DOF) of the block

src/model/Model.cpp

@@ -158,6 +158,13 @@ int Model::add_parameter(const std::vector<double>& times,
   return parameter_count++;
 }
 
+int Model::add_parameter(const std::string expression_string) {
+  auto param = Parameter(parameter_count, expression_string);
+  parameter_values.push_back(param.get(0.0));
+  parameters.push_back(std::move(param));
+  return parameter_count++;
+}
+
 Parameter* Model::get_parameter(int param_id) { return &parameters[param_id]; }
 
 double Model::get_parameter_value(int param_id) const {
@@ -181,6 +188,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

@@ -186,6 +186,14 @@ class Model {
   int add_parameter(const std::vector<double>& times,
                     const std::vector<double>& values, bool periodic = true);
 
+  /**
+   * @brief Add a function-based model parameter
+   *
+   * @param expression_string String with time-dependent function expression
+   * @return int Global ID of the parameter
+   */
+  int add_parameter(const std::string expression_string);
+
   /**
    * @brief Get a parameter by its global ID
    *

src/model/Parameter.cpp

@@ -2,6 +2,36 @@
 // University of California, and others. SPDX-License-Identifier: BSD-3-Clause
 #include "Parameter.h"
 
+#include <cstdio>
+#include <string>
+
+// Suppress -Wshadow warnings from the exprtk third-party header
+#if defined(__clang__)
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wshadow"
+#elif defined(__GNUC__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wshadow"
+#endif
+#include <exprtk.hpp>
+#if defined(__clang__)
+#pragma clang diagnostic pop
+#elif defined(__GNUC__)
+#pragma GCC diagnostic pop
+#endif
+
+/// PIMPL implementation holding the exprtk expression state
+struct Parameter::ExprtkImpl {
+  exprtk::symbol_table<double>
+      symbol_table;  ///< Symbol table mapping 't' to current time
+  exprtk::expression<double> expression;  ///< Compiled expression object
+  double* time_value = nullptr;  ///< Pointer into the symbol table for 't'
+};
+
+Parameter::~Parameter() = default;
+Parameter::Parameter(Parameter&&) noexcept = default;
+Parameter& Parameter::operator=(Parameter&&) noexcept = default;
+
 Parameter::Parameter(int id, double value) {
   this->id = id;
   update(value);
@@ -14,6 +44,12 @@ Parameter::Parameter(int id, const std::vector<double>& times,
   update(times, values);
 }
 
+Parameter::Parameter(int id, const std::string expression_string) {
+  this->id = id;
+  this->expression_string = expression_string;
+  update(expression_string);
+}
+
 void Parameter::update(double update_value) {
   is_constant = true;
   is_periodic = true;
@@ -35,6 +71,48 @@ void Parameter::update(const std::vector<double>& update_times,
   }
 }
 
+void Parameter::update(const std::string update_string) {
+  is_function = true;
+  is_constant = false;
+  expression_string = update_string;
+
+  exprtk_impl_ = std::make_unique<ExprtkImpl>();
+
+  if (!exprtk_impl_->symbol_table.create_variable("t")) {
+    throw std::runtime_error(
+        "Error failed to create time_value in symbol_table.");
+  }
+
+  exprtk_impl_->symbol_table.add_constants();
+
+  exprtk_impl_->time_value =
+      &exprtk_impl_->symbol_table.get_variable("t")->ref();
+  exprtk_impl_->expression.register_symbol_table(exprtk_impl_->symbol_table);
+
+  exprtk::parser<double> parser;
+
+  if (!parser.compile(expression_string, exprtk_impl_->expression)) {
+    is_function = false;
+    exprtk_impl_.reset();
+
+    printf("Error: %s\tExpression: %s\n", parser.error().c_str(),
+           expression_string.c_str());
+
+    for (std::size_t i = 0; i < parser.error_count(); ++i) {
+      const auto error = parser.get_error(i);
+
+      printf(
+          "Error: %02d  Position: %02d Type: [%14s] Msg: %s\tExpression: %s\n",
+          static_cast<unsigned int>(i),
+          static_cast<unsigned int>(error.token.position),
+          exprtk::parser_error::to_str(error.mode).c_str(),
+          error.diagnostic.c_str(), expression_string.c_str());
+    }
+    throw std::runtime_error(
+        "Error when compiling the function provided in 'fn'.");
+  }
+}
+
 double Parameter::get(double time) {
   // Return the constant value if parameter is constant
   if (is_constant) {
@@ -51,6 +129,13 @@ double Parameter::get(double time) {
     rtime = time;
   }
 
+  if (is_function) {
+    assert(exprtk_impl_ != nullptr);
+    assert(exprtk_impl_->time_value != nullptr);
+    *exprtk_impl_->time_value = time;
+    return exprtk_impl_->expression.value();
+  }
+
   // Determine the lower and upper element for interpolation
   auto i = lower_bound(times.begin(), times.end(), rtime);
   int k = i - times.begin();

src/model/Parameter.h

@@ -10,7 +10,9 @@
 #include <algorithm>
 #include <cmath>
 #include <iostream>
+#include <memory>
 #include <numeric>
+#include <string>
 #include <vector>
 
 #include "DOFHandler.h"
@@ -43,6 +45,30 @@ class Parameter {
   Parameter(int id, const std::vector<double>& times,
             const std::vector<double>& values, bool periodic = true);
 
+  /**
+   * @brief Construct a new Parameter object
+   *
+   * @param id Global ID of the parameter
+   * @param expression_string string with time-dependent function to get values
+   */
+  Parameter(int id, const std::string expression_string);
+
+  /**
+   * @brief Destructor
+   */
+  ~Parameter();
+
+  /**
+   * @brief Move constructor
+   */
+  Parameter(Parameter&&) noexcept;
+
+  /**
+   * @brief Move assignment operator
+   * @return Reference to this object
+   */
+  Parameter& operator=(Parameter&&) noexcept;
+
   int id;                      ///< Global ID of the parameter
   std::vector<double> times;   ///< Time steps if parameter is time-dependent
   std::vector<double> values;  ///< Values if parameter is time-dependent
@@ -53,6 +79,9 @@ class Parameter {
   bool is_constant;  ///< Bool value indicating if the parameter is constant
   bool is_periodic;  ///< Bool value indicating if the parameter is periodic
                      ///< with the cardiac cycle
+  bool is_function = false;  ///< Bool value indicating if the parameter is a
+                             ///< function
+  std::string expression_string;  ///< String with value function
 
   /**
    * @brief Update the parameter
@@ -70,6 +99,13 @@ class Parameter {
   void update(const std::vector<double>& times,
               const std::vector<double>& values);
 
+  /**
+   * @brief Update the parameter
+   *
+   * @param update_string String with function
+   */
+  void update(const std::string update_string);
+
   /**
    * @brief Get the parameter value at the specified time.
    *
@@ -92,6 +128,8 @@ class Parameter {
 
  private:
   bool steady_converted = false;
+  struct ExprtkImpl;                         ///< Forward declared PIMPL struct
+  std::unique_ptr<ExprtkImpl> exprtk_impl_;  ///< PIMPL for exprtk types
 };
 
 /**
@@ -101,6 +139,7 @@ struct InputParameter {
   bool is_optional;    ///< Is this parameter optional?
   bool is_array;       ///< Is this parameter an array?
   bool is_number;      ///< Is this parameter a number?
+  bool is_function;    ///< Is this parameter a time-dependent function?
   double default_val;  ///< Default value (if parameter is optional)
 
   /**
@@ -109,13 +148,16 @@ struct InputParameter {
    * @param is_optional Is this parameter optional?
    * @param is_array Is this parameter an array?
    * @param is_number Is this parameter a number?
+   * @param is_function Is this parameter a time-dependent function?
    * @param default_val Default value (if parameter is optional)
    */
   InputParameter(bool is_optional = false, bool is_array = false,
-                 bool is_number = true, double default_val = 0.0)
+                 bool is_number = true, bool is_function = false,
+                 double default_val = 0.0)
       : is_optional(is_optional),
         is_array(is_array),
         is_number(is_number),
+        is_function(is_function),
         default_val(default_val) {}
 };
 

src/model/PressureReferenceBC.h

@@ -85,8 +85,9 @@ class PressureReferenceBC : public Block {
    */
   PressureReferenceBC(int id, Model* model)
       : Block(id, model, BlockType::pressure_bc, BlockClass::boundary_condition,
-              {{"t", InputParameter(false, true)},
-               {"P", InputParameter(false, true)}}) {}
+              {{"t", InputParameter(true, true)},
+               {"P", InputParameter(true, true)},
+               {"fn", InputParameter(true, false, false, true)}}) {}
 
   /**
    * @brief Set up the degrees of freedom (DOF) of the block