Spatially variable robin bug fixes (SimVascular#446)

Author: aabrown100-git

Code/Source/solver/BoundaryCondition.cpp

@@ -40,14 +40,15 @@
 
 #define n_debug_bc
 
-BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std::vector<std::string>& array_names, const StringBoolMap& flags, const faceType& face)
+BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std::vector<std::string>& array_names, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger)
     : face_(&face)
     , global_num_nodes_(face.nNo)
     , local_num_nodes_(0)
     , array_names_(array_names)
     , spatially_variable(true)
     , vtp_file_path_(vtp_file_path)
     , flags_(flags)
+    , logger_(&logger)
 {
     try {
         global_data_ = read_data_from_vtp_file(vtp_file_path, array_names);
@@ -73,13 +74,14 @@ BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std
     }
 }
 
-BoundaryCondition::BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face)
+BoundaryCondition::BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger)
     : face_(&face)
     , global_num_nodes_(face.nNo)
     , local_num_nodes_(0)
     , spatially_variable(false)
     , vtp_file_path_("")
     , flags_(flags)
+    , logger_(&logger)
 {
     try {
         // Store array names, validate and store values
@@ -107,6 +109,7 @@ BoundaryCondition::BoundaryCondition(const BoundaryCondition& other)
     , vtp_file_path_(other.vtp_file_path_)
     , flags_(other.flags_)
     , global_node_map_(other.global_node_map_)
+    , logger_(other.logger_)
 {
     if (other.vtp_data_) {
         vtp_data_ = std::make_unique<VtkVtpData>(*other.vtp_data_);
@@ -126,6 +129,7 @@ void swap(BoundaryCondition& lhs, BoundaryCondition& rhs) noexcept {
     swap(lhs.flags_, rhs.flags_);
     swap(lhs.global_node_map_, rhs.global_node_map_);
     swap(lhs.vtp_data_, rhs.vtp_data_);
+    swap(lhs.logger_, rhs.logger_);
 }
 
 BoundaryCondition& BoundaryCondition::operator=(BoundaryCondition other) {
@@ -145,6 +149,7 @@ BoundaryCondition::BoundaryCondition(BoundaryCondition&& other) noexcept
     , flags_(std::move(other.flags_))
     , global_node_map_(std::move(other.global_node_map_))
     , vtp_data_(std::move(other.vtp_data_))
+    , logger_(other.logger_)
 {
     other.face_ = nullptr;
     other.global_num_nodes_ = 0;
@@ -203,6 +208,11 @@ BoundaryCondition::StringArrayMap BoundaryCondition::read_data_from_vtp_file(con
         #endif
     }
 
+    logger_ -> log_message("[BoundaryCondition] Loaded from VTP file");
+    logger_ -> log_message("\t File path:", vtp_file_path);
+    logger_ -> log_message("\t Arrays:", array_names);
+    logger_ -> log_message("\t Face:", face_->name);
+
     return result;
 }
 
@@ -375,11 +385,20 @@ void BoundaryCondition::distribute_spatially_variable(const ComMod& com_mod, con
         // Get VTP points for position matching
         Array<double> vtp_points = vtp_data_->get_points();
 
-         // Look up data for all nodes using point matching
+        // Get mesh scale factor from the face's mesh
+        double mesh_scale_factor = 1.0; // Default scale factor
+        if (face_ != nullptr) {
+            mesh_scale_factor = com_mod.msh[face_->iM].scF;
+            #ifdef debug_distribute_spatially_variable
+            dmsg << "Mesh scale factor: " << mesh_scale_factor << std::endl;
+            #endif
+        }
+        
+        // Look up data for all nodes using point matching
         for (const auto& array_name : array_names_) {
             all_values[array_name].resize(total_num_nodes);
             for (int i = 0; i < total_num_nodes; i++) {
-                int vtp_idx = find_vtp_point_index(all_positions(0,i), all_positions(1,i), all_positions(2,i), vtp_points);
+                int vtp_idx = find_vtp_point_index(all_positions(0,i), all_positions(1,i), all_positions(2,i), vtp_points, mesh_scale_factor);
                 all_values[array_name](i) = global_data_[array_name](vtp_idx, 0);
             }
         }
@@ -474,10 +493,13 @@ void BoundaryCondition::distribute_flags(const CmMod& cm_mod, const cmType& cm,
 }
 
 int BoundaryCondition::find_vtp_point_index(double x, double y, double z,
-                                const Array<double>& vtp_points) const
+                                const Array<double>& vtp_points, double mesh_scale_factor) const
 {
     const int num_points = vtp_points.ncols();
-    Vector<double> target_point{x, y, z};
+    
+    // Scale down the target coordinates to match the unscaled VTP coordinates
+    // The simulation coordinates are scaled by mesh_scale_factor, but VTP coordinates are not
+    Vector<double> target_point{x / mesh_scale_factor, y / mesh_scale_factor, z / mesh_scale_factor};
 
     // Simple linear search through all points in the VTP file
     for (int i = 0; i < num_points; i++) {
@@ -490,6 +512,7 @@ int BoundaryCondition::find_vtp_point_index(double x, double y, double z,
             #ifdef debug_bc_find_vtp_point_index
             DebugMsg dmsg(__func__, 0);
             dmsg << "Found VTP point index for node at position (" << x << ", " << y << ", " << z << ")" << std::endl;
+            dmsg << "Scaled target position (" << target_point(0) << ", " << target_point(1) << ", " << target_point(2) << ")" << std::endl;
             dmsg << "VTP point index: " << i << std::endl;
             #endif
 

Code/Source/solver/BoundaryCondition.h

@@ -36,6 +36,7 @@
 #include "CmMod.h"
 #include "VtkData.h"
 #include <string>
+#include "SimulationLogger.h"
 #include <memory>
 #include <map>
 #include <vector>
@@ -79,7 +80,7 @@ class BoundaryCondition {
     using StringDoubleMap = std::map<std::string, double>;
 
     /// @brief Data members for BC
-    const faceType* face_ = nullptr;         ///< Face associated with the BC (can be null)
+    const faceType* face_ = nullptr;         ///< Face associated with the BC (not owned by BoundaryCondition)
     int global_num_nodes_ = 0;               ///< Global number of nodes on the face
     int local_num_nodes_ = 0;                ///< Local number of nodes on this processor
     std::vector<std::string> array_names_;   ///< Names of arrays to read from VTP file
@@ -90,6 +91,7 @@ class BoundaryCondition {
     std::string vtp_file_path_;              ///< Path to VTP file (empty if uniform)
     std::map<int, int> global_node_map_;     ///< Maps global node IDs to local array indices
     std::unique_ptr<VtkVtpData> vtp_data_;   ///< VTP data object
+    const SimulationLogger* logger_ = nullptr;  ///< Logger for warnings/info (not owned by BoundaryCondition)
 
 public:
     /// @brief Tolerance for point matching in VTP files
@@ -102,13 +104,15 @@ class BoundaryCondition {
     /// @param vtp_file_path Path to VTP file containing arrays
     /// @param array_names Names of arrays to read from VTP file
     /// @param face Face associated with the BC
+    /// @param logger Simulation logger used to write warnings
     /// @throws std::runtime_error if file cannot be read or arrays are missing
-    BoundaryCondition(const std::string& vtp_file_path, const std::vector<std::string>& array_names, const StringBoolMap& flags, const faceType& face);
+    BoundaryCondition(const std::string& vtp_file_path, const std::vector<std::string>& array_names, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger);
 
     /// @brief Constructor for uniform values
     /// @param uniform_values Map of array names to uniform values
     /// @param face Face associated with the BC
-    BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face);
+    /// @param logger Simulation logger used to write warnings
+    BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger);
 
     /// @brief Copy constructor
     BoundaryCondition(const BoundaryCondition& other);
@@ -198,9 +202,10 @@ class BoundaryCondition {
     /// @param y Y coordinate
     /// @param z Z coordinate
     /// @param vtp_points VTP points array
+    /// @param mesh_scale_factor Scale factor applied to mesh coordinates
     /// @return Index of the matching point in the VTP array
     /// @throws std::runtime_error if no matching point is found
-    int find_vtp_point_index(double x, double y, double z, const Array<double>& vtp_points) const;
+    int find_vtp_point_index(double x, double y, double z, const Array<double>& vtp_points, double mesh_scale_factor) const;
 
     /// @brief Hook for derived classes to validate array values
     /// @param array_name Name of the array being validated

Code/Source/solver/Parameters.cpp

@@ -443,7 +443,6 @@ BoundaryConditionParameters::BoundaryConditionParameters()
   set_parameter("Spatial_profile_file_path", "", !required, spatial_profile_file_path);
   set_parameter("Spatial_values_file_path", "", !required, spatial_values_file_path);
   set_parameter("Stiffness", 1.0, !required, stiffness);
-  set_parameter("Robin_vtp_file_path", "", !required, robin_vtp_file_path);
   set_parameter("svZeroDSolver_block", "", !required, svzerod_solver_block);
 
   set_parameter("Temporal_and_spatial_values_file_path", "", !required, temporal_and_spatial_values_file_path);

Code/Source/solver/Parameters.h

@@ -806,7 +806,6 @@ class BoundaryConditionParameters : public ParameterLists
     Parameter<std::string> spatial_profile_file_path;
     Parameter<std::string> spatial_values_file_path;
     Parameter<double> stiffness;
-    Parameter<std::string> robin_vtp_file_path;
     Parameter<std::string> svzerod_solver_block;
 
     Parameter<std::string> temporal_and_spatial_values_file_path;

Code/Source/solver/RobinBoundaryCondition.cpp

@@ -29,6 +29,17 @@
  */
 
 #include "RobinBoundaryCondition.h"
+#include <iostream>
 
 #define n_debug_robin_bc
 
+RobinBoundaryCondition::RobinBoundaryCondition(double uniform_stiffness, double uniform_damping, bool normal_only, const faceType& face, SimulationLogger& logger)
+    : BoundaryCondition({{"Stiffness", uniform_stiffness}, {"Damping", uniform_damping}}, StringBoolMap{{"normal_direction_only", normal_only}}, face, logger)
+{
+    // Warning if both stiffness and damping are zero
+    if (uniform_stiffness == 0.0 && uniform_damping == 0.0) {
+        logger_ -> log_message("WARNING [RobinBoundaryCondition] Both stiffness and damping values set to zero. "
+                   "This will result in effectively no boundary condition being applied.");
+    }
+}
+

Code/Source/solver/RobinBoundaryCondition.h

@@ -70,21 +70,22 @@ class RobinBoundaryCondition : public BoundaryCondition {
     /// @param vtp_file_path Path to VTP file containing Stiffness and Damping point arrays
     /// @param normal_only Flag to apply only along normal direction
     /// @param face Face associated with the Robin BC
+    /// @param logger Simulation logger used to write warnings
     /// @throws BoundaryConditionFileException if file cannot be read
     /// @throws BoundaryConditionVtpArrayException if arrays are missing
     /// @throws BoundaryConditionValidationException if values are invalid
-    RobinBoundaryCondition(const std::string& vtp_file_path, bool normal_only, const faceType& face)
-        : BoundaryCondition(vtp_file_path, std::vector<std::string>{"Stiffness", "Damping"}, StringBoolMap{{"normal_direction_only", normal_only}}, face) {}
+    RobinBoundaryCondition(const std::string& vtp_file_path, bool normal_only, const faceType& face, SimulationLogger& logger)
+        : BoundaryCondition(vtp_file_path, std::vector<std::string>{"Stiffness", "Damping"}, StringBoolMap{{"normal_direction_only", normal_only}}, face, logger) {}
 
 
     /// @brief Constructor for uniform values
     /// @param uniform_stiffness Uniform stiffness value for all nodes
     /// @param uniform_damping Uniform damping value for all nodes
     /// @param normal_only Flag to apply only along normal direction
     /// @param face Face associated with the Robin BC
+    /// @param logger Simulation logger used to write warnings
     /// @throws BoundaryConditionValidationException if values are invalid
-    RobinBoundaryCondition(double uniform_stiffness, double uniform_damping, bool normal_only, const faceType& face)
-        : BoundaryCondition({{"Stiffness", uniform_stiffness}, {"Damping", uniform_damping}}, StringBoolMap{{"normal_direction_only", normal_only}}, face) {}
+    RobinBoundaryCondition(double uniform_stiffness, double uniform_damping, bool normal_only, const faceType& face, SimulationLogger& logger);
 
     /// @brief Apply only along normal direction (getter)
     /// @return true if BC should be applied only along normal direction

Code/Source/solver/SimulationLogger.h

@@ -47,7 +47,7 @@ class SimulationLogger {
       this->initialize(file_name, cout_write);
     }
 
-    void initialize(const std::string& file_name, bool cout_write=false) 
+    void initialize(const std::string& file_name, bool cout_write=false) const
     {
       log_file_.open(file_name);
       if (log_file_.fail()) {
@@ -58,12 +58,14 @@ class SimulationLogger {
       file_name_ = file_name;
     }
 
+    bool is_initialized() const { return log_file_.is_open(); }
+
     ~SimulationLogger() 
     {
       log_file_.close();
     }
 
-  template <class T> SimulationLogger& operator<< (const T& value)
+  template <class T> const SimulationLogger& operator<< (const T& value) const
   {
     if (file_name_ == "") { 
       return *this;
@@ -78,7 +80,28 @@ class SimulationLogger {
     return *this;
   }
 
-  SimulationLogger& operator<<(std::ostream&(*value)(std::ostream& o))
+  // Special handling for vector<string>
+  const SimulationLogger& operator<< (const std::vector<std::string>& values) const
+  {
+    if (file_name_ == "") { 
+      return *this;
+    }
+
+    bool first = true;
+    for (const auto& value : values) {
+      if (!first) {
+        log_file_ << ", ";
+        if (cout_write_) std::cout << ", ";
+      }
+      log_file_ << value;
+      if (cout_write_) std::cout << value;
+      first = false;
+    }
+
+    return *this;
+  }
+
+  const SimulationLogger& operator<<(std::ostream&(*value)(std::ostream& o)) const
   {
     if (file_name_ == "") { 
       return *this;
@@ -93,10 +116,27 @@ class SimulationLogger {
     return *this;
   };
 
+  /// @brief Log a message with automatic space separation
+  /// @param args Arguments to log
+  template<typename... Args>
+  const SimulationLogger& log_message(const Args&... args) const {
+    if (file_name_ == "") return *this;
+
+    bool first = true;
+    ((first ? (first = false, (*this << args)) : (*this << ' ' << args)), ...);
+    *this << std::endl;
+    return *this;
+  }
+
   private:
-    bool cout_write_ = false;
-    std::string file_name_;
-    std::ofstream log_file_;
+    // These members are marked mutable because they can be modified in const member functions.
+    // While logging writes to a file (modifying these members), it doesn't change the logical
+    // state of the logger - this is exactly what mutable is designed for: implementation
+    // details that need to change even when the object's public state remains constant.
+    mutable bool cout_write_ = false;
+    mutable bool initialized_ = false;
+    mutable std::string file_name_;
+    mutable std::ofstream log_file_;
 };
 
 #endif

Code/Source/solver/initialize.cpp

@@ -371,22 +371,6 @@ void initialize(Simulation* simulation, Vector<double>& timeP)
   dmsg.banner();
   #endif
 
-  // Setup logging to history file.
-  if (!simulation->com_mod.cm.slv(simulation->cm_mod)) {
-    std::string hist_file_name;
-
-    if (chnl_mod.appPath != "") { 
-      auto mkdir_arg = std::string("mkdir -p ") + chnl_mod.appPath;
-      std::system(mkdir_arg.c_str());
-      hist_file_name = chnl_mod.appPath + "/" + simulation->history_file_name;
-    } else {
-      hist_file_name = simulation->history_file_name;
-    }
-
-    bool output_to_cout = true;
-    simulation->logger.initialize(hist_file_name, output_to_cout);
-  }
-
   #ifdef debug_initialize
   dmsg << "Set time " << std::endl;
   dmsg << "com_mod.timeP[0]: " << com_mod.timeP[0] << std::endl;

Code/Source/solver/read_files.cpp

@@ -394,15 +394,17 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
   if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Robin))) { 
 
     // Read VTP file path for per-node stiffness and damping (optional)
-    if (bc_params->robin_vtp_file_path.defined()) {
-      lBc.robin_bc = RobinBoundaryCondition(bc_params->robin_vtp_file_path.value(), 
+    if (bc_params->spatial_values_file_path.defined()) {
+      lBc.robin_bc = RobinBoundaryCondition(bc_params->spatial_values_file_path.value(), 
                                             bc_params->apply_along_normal_direction.value(),
-                                            com_mod.msh[lBc.iM].fa[lBc.iFa]);
+                                            com_mod.msh[lBc.iM].fa[lBc.iFa],
+                                            simulation->logger);
     } else {
       lBc.robin_bc = RobinBoundaryCondition(bc_params->stiffness.value(), 
                                             bc_params->damping.value(),
                                             bc_params->apply_along_normal_direction.value(),
-                                            com_mod.msh[lBc.iM].fa[lBc.iFa]);
+                                            com_mod.msh[lBc.iM].fa[lBc.iFa],
+                                            simulation->logger);
     }
   }
 
@@ -1703,7 +1705,23 @@ void read_files(Simulation* simulation, const std::string& file_name)
     com_mod.stFileFlag = gen_params.continue_previous_simulation.value();
   }
 
-  // Set simulatioin and module member data from XML parameters.
+  // Setup logging to history file.
+  if (!simulation->com_mod.cm.slv(simulation->cm_mod)) {
+    std::string hist_file_name;
+
+    if (chnl_mod.appPath != "") { 
+      auto mkdir_arg = std::string("mkdir -p ") + chnl_mod.appPath;
+      std::system(mkdir_arg.c_str());
+      hist_file_name = chnl_mod.appPath + "/" + simulation->history_file_name;
+    } else {
+      hist_file_name = simulation->history_file_name;
+    }
+
+    bool output_to_cout = true;
+    simulation->logger.initialize(hist_file_name, output_to_cout);
+  }
+
+  // Set simulation and module member data from XML parameters.
   simulation->set_module_parameters();
 
   // Read mesh and BCs data.

tests/cases/struct/spatially_variable_robin/solver.xml

@@ -132,7 +132,7 @@
 
    <Add_BC name="Y0" > 
       <Type> Robin </Type> 
-      <Robin_vtp_file_path> Y0_spatially_varying_robin.vtp </Robin_vtp_file_path> 
+      <Spatial_values_file_path> Y0_spatially_varying_robin.vtp </Spatial_values_file_path> 
       <Apply_along_normal_direction> false </Apply_along_normal_direction> 
    </Add_BC>