Variable Robin boundary condition and generic boundary condition class

Code/Source/solver/BoundaryCondition.h

@@ -0,0 +1,273 @@
+/* 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 BOUNDARY_CONDITION_H
+#define BOUNDARY_CONDITION_H
+
+#include "Array.h"
+#include "Vector.h"
+#include "CmMod.h"
+#include "VtkData.h"
+#include <string>
+#include <memory>
+#include <map>
+#include <vector>
+#include <stdexcept>
+
+// Forward declarations. These are needed because including ComMod.h causes a 
+// circular header dependency.
+class faceType;
+class ComMod;
+
+/// @brief Base class for boundary conditions with spatially variable arrays
+/// 
+/// This class provides common functionality for boundary conditions that need to
+/// read and manage arrays of values from VTP files or uniform values. It handles
+/// distribution of data across processes and provides efficient access to values.
+///
+/// This class is intended to be subclassed by specific boundary condition types.
+///
+/// Development note: this class is intended to eventually be an object-oriented
+/// replacement of the existing bcType, although it is not yet complete.
+///
+/// Example usage:
+/// ```cpp
+/// class MyBoundaryCondition : public BoundaryCondition {
+/// public:
+///     MyBoundaryCondition(const std::string& vtp_file_path, const std::vector<std::string>& array_names, const faceType& face)
+///         : BoundaryCondition(vtp_file_path, array_names, face) {}
+///
+///     MyBoundaryCondition(const std::map<std::string, double>& uniform_values, const faceType& face)
+///         : BoundaryCondition(uniform_values, face) {}
+///
+///     // Add BC-specific functionality here
+/// };
+/// ```
+class BoundaryCondition {
+protected:
+    /// @brief Type alias for map of array names to array data
+    using StringArrayMap = std::map<std::string, Array<double>>;
+    using StringBoolMap = std::map<std::string, bool>;
+    using StringDoubleMap = std::map<std::string, double>;
+
+    /// @brief Data members for BC
+    const faceType* face_;                   ///< Face associated with the BC (can be null)
+    int global_num_nodes_;                   ///< Global number of nodes on the face
+    int local_num_nodes_;                    ///< Local number of nodes on this processor
+    std::vector<std::string> array_names_;   ///< Names of arrays to read from VTP file
+    StringArrayMap local_data_;              ///< Local array values for each node on this processor
+    StringArrayMap global_data_;             ///< Global array values (only populated on master)
+    StringBoolMap flags_;                    ///< Named boolean flags for BC behavior
+    bool spatially_variable;                 ///< Flag indicating if data is from VTP file
+    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
+    bool defined_;                           ///< Whether this BC has been properly initialized
+
+public:
+    /// @brief Tolerance for point matching in VTP files
+    static constexpr double POINT_MATCH_TOLERANCE = 1e-12;
+
+    /// @brief Default constructor - creates an empty BC
+    BoundaryCondition() : face_(nullptr), global_num_nodes_(0), local_num_nodes_(0), spatially_variable(false), defined_(false) {}
+
+    /// @brief Constructor - reads data from VTP file
+    /// @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
+    /// @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);
+
+    /// @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);
+
+    /// @brief Copy constructor
+    BoundaryCondition(const BoundaryCondition& other);
+
+    /// @brief Copy assignment operator
+    BoundaryCondition& operator=(const BoundaryCondition& other);
+
+    /// @brief Move constructor
+    BoundaryCondition(BoundaryCondition&& other) noexcept;
+
+    /// @brief Move assignment operator
+    BoundaryCondition& operator=(BoundaryCondition&& other) noexcept;
+
+    /// @brief Virtual destructor
+    virtual ~BoundaryCondition() = default;
+
+
+    /// @brief Get value for a specific array and node
+    /// @param array_name Name of the array
+    /// @param node_id Node index on the face
+    /// @return Value for the array at the specified node
+    double get_value(const std::string& array_name, int node_id) const;
+
+    /// @brief Get a boolean flag by name
+    /// @param name Name of the flag
+    /// @return Value of the flag
+    /// @throws std::runtime_error if flag is not found
+    bool get_flag(const std::string& name) const;
+
+    /// @brief Get a string representation of the flags
+    /// @return String representation of the flags
+    std::string flags_to_string() const;
+
+    /// @brief Get global number of nodes
+    /// @return Global number of nodes on the face
+    int get_global_num_nodes() const {
+        return global_num_nodes_;
+    }
+
+    /// @brief Get local number of nodes
+    /// @return Local number of nodes on the face on this processor
+    int get_local_num_nodes() const {
+        return local_num_nodes_;
+    }
+
+    /// @brief Get local array index for a global node ID
+    /// @param global_node_id The global node ID defined on the face
+    /// @return Local array index for data arrays
+    /// @throws std::runtime_error if global_node_id is not found in the map
+    int get_local_index(int global_node_id) const;
+
+    /// @brief Check if data is loaded from VTP file
+    /// @return true if loaded from VTP, false if using uniform values
+    bool is_from_vtp() const {
+        return spatially_variable;
+    }
+
+    /// @brief Get the VTP file path (empty if using uniform values)
+    /// @return VTP file path
+    const std::string& get_vtp_path() const {
+        return vtp_file_path_;
+    }
+
+    /// @brief Check if this BC is properly defined
+    /// @return true if BC has been initialized with either VTP data or uniform values
+    bool is_defined() const {
+        return defined_;
+    }
+
+    /// @brief Distribute BC data from the master process to the slave processes
+    /// @param com_mod Reference to ComMod object for global coordinates
+    /// @param cm_mod Reference to CmMod object for MPI communication
+    /// @param cm Reference to cmType object for MPI communication
+    /// @param face Face associated with the BC
+    void distribute(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const faceType& face);
+
+protected:
+    /// @brief Read data from VTP file
+    /// @param vtp_file_path Path to VTP file
+    /// @param array_names Names of arrays to read
+    /// @return Map of array names to array data
+    StringArrayMap read_data_from_vtp_file(const std::string& vtp_file_path, const std::vector<std::string>& array_names);
+
+    /// @brief Find index of a point in the VTP points array
+    /// @param x X coordinate
+    /// @param y Y coordinate
+    /// @param z Z coordinate
+    /// @param vtp_points VTP points array
+    /// @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;
+
+    /// @brief Hook for derived classes to validate array values
+    /// @param array_name Name of the array being validated
+    /// @param value Value to validate
+    /// @throws std::runtime_error if validation fails
+    virtual void validate_array_value(const std::string& array_name, double value) const {}
+
+     // ---- distribute helpers ----
+     void distribute_metadata(const CmMod& cm_mod, const cmType& cm, bool is_slave);
+     void distribute_spatially_variable(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bool is_slave);
+     void distribute_uniform(const CmMod& cm_mod, const cmType& cm, bool is_slave);
+     void distribute_flags(const CmMod& cm_mod, const cmType& cm, bool is_slave);
+
+};
+
+/// @brief Base exception class for BC errors
+class BoundaryConditionBaseException : public std::exception {
+public:
+    /// @brief Constructor
+    /// @param msg Error message
+    explicit BoundaryConditionBaseException(const std::string& msg) : message_(msg) {}
+
+    /// @brief Get error message
+    /// @return Error message
+    const char* what() const noexcept override {
+        return message_.c_str();
+    }
+
+private:
+    std::string message_;
+};
+
+/// @brief Exception thrown when VTP file cannot be read or is invalid
+class BoundaryConditionFileException : public BoundaryConditionBaseException {
+public:
+    /// @brief Constructor
+    /// @param file Path to VTP file
+    explicit BoundaryConditionFileException(const std::string& file)
+        : BoundaryConditionBaseException("Failed to open or read the VTP file '" + file + "'") {}
+};
+
+/// @brief Exception thrown when node count mismatch between VTP and face
+class BoundaryConditionNodeCountException : public BoundaryConditionBaseException {
+public:
+    /// @brief Constructor
+    /// @param vtp_file VTP file path
+    /// @param face_name Face name
+    explicit BoundaryConditionNodeCountException(const std::string& vtp_file, const std::string& face_name)
+        : BoundaryConditionBaseException("Number of nodes in VTP file '" + vtp_file +
+                                       "' does not match number of nodes on face '" + face_name + "'") {}
+};
+
+/// @brief Exception thrown when array validation fails
+class BoundaryConditionValidationException : public BoundaryConditionBaseException {
+public:
+    /// @brief Constructor
+    /// @param array_name Name of array that failed validation
+    /// @param value Value that failed validation
+    explicit BoundaryConditionValidationException(const std::string& array_name, double value) 
+        : BoundaryConditionBaseException("Invalid value " + std::to_string(value) + 
+                                       " for array '" + array_name + "'") {}
+};
+
+/// @brief Exception thrown when a requested flag is not defined
+class BoundaryConditionFlagException : public BoundaryConditionBaseException {
+public:
+    explicit BoundaryConditionFlagException(const std::string& flag_name)
+        : BoundaryConditionBaseException("BoundaryCondition flag not found: '" + flag_name + "'") {}
+};
+
+#endif // BOUNDARY_CONDITION_H

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
     //------------