Variable Robin boundary condition and generic boundary condition class (SimVascular#437)

* Initial implementation and test case. Works with 1 processor, but fails with multiple processors. Other improvements to make

* Put debugging statements inside debugging blocks

* Add spatially_variable_robin test to pytest suite

* Rename to VariableRobinBCData, initial VariableRobinBCData once, instead of at every Newton iteration

* Add copyright statement

* Renamed to RobinBC, partial implementation of parallelization

* Refactored distribute method to use gatherv and scatterv. Dealing with some compilation issues now.

* Passes test in parallel, fails with index error with 1 proc

* Fixed bug, now passes test in serial

* clean up debugging statements to make code more readable

* Refactor vtp file reading to accept vector of array names, more generic this way

* Split implementation into a BC base class and a RobinBC derived class

* Update README.md and add ustruct test case. ustruct result is fairly different from struct result, so could not use result_002.vtu from the struct case

* Refactor BC class to improve error handling and code clarity. Introduced custom exception classes for VTP file errors and node count mismatches. Updated point matching logic to use squared distance for efficiency. Removed redundant map clearing operations during initialization. Also, put test cases number of timesteps back to 2

* Rename BC and RobinBC classes as BoundaryCondition and RobinBoundaryCondition classes.

* Clean up constructors

* Remove extraneous comments and debugging statements

* Reformat code

* Move flag to apply in normal direction only inside RobinBoundaryCondition. Add generic flag handling in BoundaryCondition

* Split up distribute method into smaller functions

* Add copyrights

* Implement copy-and-swap idiom, add noexcept to appropriate methods

* Catch exceptions in constructors

* Move generate_load.py and generate_spatially_variable_robin.py into utilities/generate_boundary_condition_data/ directory, and update READMEs to reflect new location

* Initialize members with default values and use default default constructor

* Remove defined_ flag, just check if arrays are allocated, add new exception classes

Author: aabrown100-git

Code/Source/solver/BoundaryCondition.cpp

@@ -224,6 +224,9 @@ set(CSRCS
   SPLIT.c
 
   svZeroD_interface/LPNSolverInterface.h svZeroD_interface/LPNSolverInterface.cpp
+
+  BoundaryCondition.h BoundaryCondition.cpp
+  RobinBoundaryCondition.h RobinBoundaryCondition.cpp
 )
 
   # Set PETSc interace code.

Code/Source/solver/BoundaryCondition.h

@@ -163,4 +163,84 @@ void cmType::bcast(const CmMod& cm_mod, int* data) const
 void cmType::bcast(const CmMod& cm_mod, Vector<int>& data) const
 {
   MPI_Bcast(data.data(), data.size(), cm_mod::mpint, cm_mod.master, com());
+}
+
+/// @brief gather int array
+void cmType::gather(const CmMod& cm_mod, const int* send_data, int send_count, int* recv_data, int recv_count, int root) const
+{
+  MPI_Gather(send_data, send_count, cm_mod::mpint, recv_data, recv_count, cm_mod::mpint, root, com());
+}
+
+/// @brief gather double array
+void cmType::gather(const CmMod& cm_mod, const double* send_data, int send_count, double* recv_data, int recv_count, int root) const
+{
+  MPI_Gather(send_data, send_count, cm_mod::mpreal, recv_data, recv_count, cm_mod::mpreal, root, com());
+}
+
+/// @brief gather int Vector
+void cmType::gather(const CmMod& cm_mod, const Vector<int>& send_data, Vector<int>& recv_data, int root) const
+{
+  MPI_Gather(send_data.data(), send_data.size(), cm_mod::mpint, 
+             recv_data.data(), send_data.size(), cm_mod::mpint, root, com());
+}
+
+/// @brief gather double Vector
+void cmType::gather(const CmMod& cm_mod, const Vector<double>& send_data, Vector<double>& recv_data, int root) const
+{
+  MPI_Gather(send_data.data(), send_data.size(), cm_mod::mpreal, 
+             recv_data.data(), send_data.size(), cm_mod::mpreal, root, com());
+}
+
+/// @brief scatter int array
+void cmType::scatter(const CmMod& cm_mod, const int* send_data, int send_count, int* recv_data, int recv_count, int root) const
+{
+  MPI_Scatter(send_data, send_count, cm_mod::mpint, recv_data, recv_count, cm_mod::mpint, root, com());
+}
+
+/// @brief scatter double array
+void cmType::scatter(const CmMod& cm_mod, const double* send_data, int send_count, double* recv_data, int recv_count, int root) const
+{
+  MPI_Scatter(send_data, send_count, cm_mod::mpreal, recv_data, recv_count, cm_mod::mpreal, root, com());
+}
+
+/// @brief scatter int Vector
+void cmType::scatter(const CmMod& cm_mod, const Vector<int>& send_data, Vector<int>& recv_data, int root) const
+{
+  MPI_Scatter(send_data.data(), send_data.size() / nProcs, cm_mod::mpint, 
+              recv_data.data(), send_data.size() / nProcs, cm_mod::mpint, root, com());
+}
+
+/// @brief scatter double Vector
+void cmType::scatter(const CmMod& cm_mod, const Vector<double>& send_data, Vector<double>& recv_data, int root) const
+{
+  MPI_Scatter(send_data.data(), send_data.size() / nProcs, cm_mod::mpreal, 
+              recv_data.data(), send_data.size() / nProcs, cm_mod::mpreal, root, com());
+}
+
+/// @brief gatherv int Vector
+void cmType::gatherv(const CmMod& cm_mod, const Vector<int>& send_data, Vector<int>& recv_data, const Vector<int>& recv_counts, const Vector<int>& displs, int root) const
+{
+  MPI_Gatherv(send_data.data(), send_data.size(), cm_mod::mpint,
+              recv_data.data(), recv_counts.data(), displs.data(), cm_mod::mpint, root, com());
+}
+
+/// @brief gatherv double Vector
+void cmType::gatherv(const CmMod& cm_mod, const Vector<double>& send_data, Vector<double>& recv_data, const Vector<int>& recv_counts, const Vector<int>& displs, int root) const
+{
+  MPI_Gatherv(send_data.data(), send_data.size(), cm_mod::mpreal,
+              recv_data.data(), recv_counts.data(), displs.data(), cm_mod::mpreal, root, com());
+}
+
+/// @brief scatterv int Vector
+void cmType::scatterv(const CmMod& cm_mod, const Vector<int>& send_data, const Vector<int>& send_counts, const Vector<int>& displs, Vector<int>& recv_data, int root) const
+{
+  MPI_Scatterv(send_data.data(), send_counts.data(), displs.data(), cm_mod::mpint,
+               recv_data.data(), recv_data.size(), cm_mod::mpint, root, com());
+}
+
+/// @brief scatterv double Vector
+void cmType::scatterv(const CmMod& cm_mod, const Vector<double>& send_data, const Vector<int>& send_counts, const Vector<int>& displs, Vector<double>& recv_data, int root) const
+{
+  MPI_Scatterv(send_data.data(), send_counts.data(), displs.data(), cm_mod::mpreal,
+               recv_data.data(), recv_data.size(), cm_mod::mpreal, root, com());
 }

Code/Source/solver/CMakeLists.txt

@@ -112,6 +112,26 @@ class cmType {
 
     void bcast(const CmMod& cm_mod, Array<int>& data, const std::string& name="") const;
 
+    // Gather operations
+    void gather(const CmMod& cm_mod, const int* send_data, int send_count, int* recv_data, int recv_count, int root) const;
+    void gather(const CmMod& cm_mod, const double* send_data, int send_count, double* recv_data, int recv_count, int root) const;
+    void gather(const CmMod& cm_mod, const Vector<int>& send_data, Vector<int>& recv_data, int root) const;
+    void gather(const CmMod& cm_mod, const Vector<double>& send_data, Vector<double>& recv_data, int root) const;
+
+    // Gatherv operations
+    void gatherv(const CmMod& cm_mod, const Vector<int>& send_data, Vector<int>& recv_data, const Vector<int>& recv_counts, const Vector<int>& displs, int root) const;
+    void gatherv(const CmMod& cm_mod, const Vector<double>& send_data, Vector<double>& recv_data, const Vector<int>& recv_counts, const Vector<int>& displs, int root) const;
+
+    // Scatterv operations
+    void scatterv(const CmMod& cm_mod, const Vector<int>& send_data, const Vector<int>& send_counts, const Vector<int>& displs, Vector<int>& recv_data, int root) const;
+    void scatterv(const CmMod& cm_mod, const Vector<double>& send_data, const Vector<int>& send_counts, const Vector<int>& displs, Vector<double>& recv_data, int root) const;
+
+    // Scatter operations
+    void scatter(const CmMod& cm_mod, const int* send_data, int send_count, int* recv_data, int recv_count, int root) const;
+    void scatter(const CmMod& cm_mod, const double* send_data, int send_count, double* recv_data, int recv_count, int root) const;
+    void scatter(const CmMod& cm_mod, const Vector<int>& send_data, Vector<int>& recv_data, int root) const;
+    void scatter(const CmMod& cm_mod, const Vector<double>& send_data, Vector<double>& recv_data, int root) const;
+
     //------------
     // bcast_enum
     //------------

Code/Source/solver/CmMod.cpp

@@ -43,6 +43,7 @@
 #include "ChnlMod.h"
 #include "CmMod.h"
 #include "Parameters.h"
+#include "RobinBoundaryCondition.h"
 #include "Timer.h"
 #include "Vector.h"
 
@@ -58,6 +59,7 @@
 
 #include <array>
 #include <iostream>
+#include <memory>
 #include <string>
 #include <vector>
 #include <fstream>
@@ -150,17 +152,13 @@ class rcrType
 class bcType
 {
   public:
-
     // Strong/Weak application of Dirichlet BC
     bool weakDir = false;
 
     // Whether load vector changes with deformation
     // (Neu - struct/ustruct only)
     bool flwP = false;
 
-    // Robin: apply only in normal direction
-    bool rbnN = false;
-
     // Strong/Weak application of Dirichlet BC
     int clsFlgRis = 0;
 
@@ -191,11 +189,8 @@ class bcType
     // Neu: defined resistance
     double r = 0.0;
 
-    // Robin: stiffness
-    double k = 0.0;
-
-    // Robin: damping
-    double c = 0.0;
+    // Robin: VTP file path for per-node stiffness and damping
+    std::string robin_vtp_file = "";
 
     // RIS0D: resistance
     double resistance = 0.0;
@@ -227,6 +222,9 @@ class bcType
 
     // Neu: RCR
     rcrType RCR;
+
+    // Robin BC class
+    RobinBoundaryCondition robin_bc;
 };
 
 /// @brief Class storing data for B-Splines.

Code/Source/solver/CmMod.h

@@ -443,6 +443,7 @@ 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/ComMod.h

@@ -806,6 +806,7 @@ 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/Parameters.cpp

@@ -0,0 +1,34 @@
+/* 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 "RobinBoundaryCondition.h"
+
+#define n_debug_robin_bc
+

Code/Source/solver/Parameters.h

@@ -0,0 +1,129 @@
+/* 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.
+ */
+
+#ifndef ROBIN_BOUNDARY_CONDITION_H
+#define ROBIN_BOUNDARY_CONDITION_H
+
+#include "BoundaryCondition.h"
+#include <string>
+#include <map>
+#include <vector>
+
+/// @brief Class to handle Robin boundary condition with potentially spatially variable arrays
+/// 
+/// This class extends the generic BoundaryCondition class to handle Robin boundary conditions, which require
+/// stiffness and damping arrays. While it supports any number of named arrays through its base class,
+/// it provides specific validation and convenience methods for stiffness and damping values.
+///
+/// Example usage:
+/// ```cpp
+/// // Read multiple arrays from VTP file
+/// std::vector<std::string> array_names = {"Stiffness", "Damping"};
+/// RobinBoundaryCondition bc(vtp_file_path, array_names, face);
+///
+/// // Access values
+/// double stiffness = bc.get_stiffness(node_id);  // Convenience method
+/// double damping = bc.get_damping(node_id);      // Convenience method
+///
+/// // Create with uniform values
+/// std::map<std::string, double> uniform_values = {
+///     {"Stiffness", 1.0},
+///     {"Damping", 0.5},
+/// };
+/// RobinBoundaryCondition bc(uniform_values, face);
+/// ```
+
+#define debug_robin_bc
+class RobinBoundaryCondition : public BoundaryCondition {
+public:
+    /// @brief Default constructor - creates an empty RobinBoundaryCondition
+    RobinBoundaryCondition() : BoundaryCondition() {}
+
+    /// @brief Constructor - reads stiffness and damping from VTP file
+    /// @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
+    /// @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) {}
+
+
+    /// @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
+    /// @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) {}
+ 
+    /// @brief Apply only along normal direction (getter)
+    /// @return true if BC should be applied only along normal direction
+    /// @throws BoundaryConditionFlagException if "normal_direction_only" flag not found
+    bool normal_direction_only() const { return this->get_flag("normal_direction_only"); }
+ 
+    /// @brief Get stiffness value for a specific node (convenience method)
+    /// @param node_id Node index on the face
+    /// @return Stiffness value for the node
+    /// @throws BoundaryConditionArrayException if "Stiffness" array not found
+    /// @throws BoundaryConditionNodeIdException if node_id is out of range
+    double get_stiffness(int node_id) const {
+        return get_value("Stiffness", node_id);
+    }
+    
+    /// @brief Get damping value for a specific node (convenience method)
+    /// @param node_id Node index on the face
+    /// @return Damping value for the node
+    /// @throws BoundaryConditionArrayException if "Damping" array not found
+    /// @throws BoundaryConditionNodeIdException if node_id is out of range
+    double get_damping(int node_id) const {
+        return get_value("Damping", node_id);
+    }
+
+    /// @brief Assemble the Robin BC into the global residual vector and stiffness matrix
+    /// Currently not implemented
+    /// @return 0
+    double assemble() const { return 0; }
+
+protected:
+    /// @brief Validate array values for Robin BC
+    /// @param array_name Name of the array being validated
+    /// @param value Value to validate
+    /// @throws BoundaryConditionValidationException if validation fails
+    void validate_array_value(const std::string& array_name, double value) const override {
+        if (value < 0.0) {
+            throw BoundaryConditionValidationException(array_name, value);
+        }
+    }
+};
+
+#endif // ROBIN_BOUNDARY_CONDITION_H