Function input feature for time-dependent boundary conditions

src/model/FlowReferenceBC.cpp

@@ -36,6 +36,12 @@ void FlowReferenceBC::setup_dofs(DOFHandler &dofhandler) {
 
 void FlowReferenceBC::update_constant(SparseSystem &system,
                                       std::vector<double> &parameters) {
+  //std::cout << system << std::endl;
+  for (double i: parameters){
+    std::cout << i << std::endl;
+  }
+  std::cout << global_eqn_ids[0] << std::endl;
+  std::cout << global_var_ids[1] << std::endl;
   system.F.coeffRef(global_eqn_ids[0], global_var_ids[1]) = 1.0;
 }
 

src/model/FlowReferenceBC.h

@@ -88,8 +88,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

@@ -178,6 +178,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 {

src/model/Model.h

@@ -205,6 +205,14 @@ class Model {
   int add_parameter(const std::vector<double> &times,
                     const std::vector<double> &values, bool periodic = true);
 
+  /**
+   * @brief Add a time-dependent model parameter defined by a function
+   *
+   * @param expression_string String with function to determine value
+   * @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

@@ -29,6 +29,10 @@
 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "Parameter.h"
+#include <cstdio>
+#include <string>
+
+#include "exprtk.hpp"
 
 Parameter::Parameter(int id, double value) {
   this->id = id;
@@ -42,6 +46,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;
@@ -63,6 +73,12 @@ void Parameter::update(const std::vector<double> &update_times,
   }
 }
 
+
+void Parameter::update(const std::string update_string) {
+  is_function = true;
+  expression_string = update_string;
+}
+
 double Parameter::get(double time) {
   // Return the constant value if parameter is constant
   if (is_constant) {
@@ -79,6 +95,29 @@ double Parameter::get(double time) {
     rtime = time;
   }
 
+  if (is_function == true) {
+    // Adapted from example from Basic Design example at http://www.partow.net/programming/exprtk/index.html
+    typedef double T;
+
+    typedef exprtk::symbol_table<T> symbol_table_t;
+    typedef exprtk::expression<T>   expression_t;
+    typedef exprtk::parser<T>       parser_t;
+
+    T t = T(time);
+
+    symbol_table_t symbol_table;
+    symbol_table.add_variable("t",t);
+
+    expression_t expression;
+    expression.register_symbol_table(symbol_table);
+
+    parser_t parser;
+
+    parser.compile(expression_string, expression);
+    T value = expression.value();
+    return 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

@@ -39,6 +39,10 @@
 #include <iostream>
 #include <numeric>
 #include <vector>
+#include <cstdio>
+#include <string>
+
+#include "exprtk.hpp"
 
 #include "DOFHandler.h"
 
@@ -70,6 +74,13 @@ 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);
   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
@@ -80,6 +91,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;           ///< Bool value indicating if the parameter is a function
+  std::string expression_string;  ///< String with value function
+
 
   /**
    * @brief Update the parameter
@@ -97,6 +111,14 @@ class Parameter {
   void update(const std::vector<double>& times,
               const std::vector<double>& values);
 
+  /**
+   * @brief Update the parameter
+   *
+   * @param expression_string String with function
+   * @param update_value Value at time step
+   */
+  void update(const std::string update_string);
+
   /**
    * @brief Get the parameter value at the specified time.
    *
@@ -128,6 +150,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)
 
   /**
@@ -139,10 +162,11 @@ struct InputParameter {
    * @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) {}
 };