Variable Robin boundary condition and generic boundary condition class

Code/Source/solver/CMakeLists.txt

@@ -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/CmMod.cpp

@@ -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/CmMod.h

@@ -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/ComMod.h

@@ -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/Parameters.cpp

@@ -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/Parameters.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/RobinBoundaryCondition.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/RobinBoundaryCondition.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