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

Author: aabrown100-git

Code/Source/solver/CMakeLists.txt

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

Code/Source/solver/ComMod.h

@@ -197,6 +197,9 @@ class bcType
     // 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;
 

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/RobinBCData.cpp

@@ -0,0 +1,194 @@
+#include "RobinBCData.h"
+#include "ComMod.h"
+#include <stdexcept>
+#include <iostream>
+
+RobinBCData::RobinBCData(const std::string& vtp_file_path, const faceType& face)
+{
+    std::cerr << "[RobinBCData] Constructor with VTP file path" << std::endl;
+    num_nodes_ = face.nNo;
+    from_vtp_ = true;
+    vtp_file_path_ = vtp_file_path;
+
+    //std::cout << "[RobinBCData] Loading Robin BC from VTP file: " << vtp_file_path << std::endl;
+    load_from_vtp(vtp_file_path, face);
+}
+
+RobinBCData::RobinBCData(double uniform_stiffness, double uniform_damping, int num_nodes)
+{
+    std::cout << "[RobinBCData] Constructor with uniform values" << std::endl;
+    num_nodes_ = num_nodes;
+    from_vtp_ = false;
+    vtp_file_path_ = "";
+
+    std::cout << "[RobinBCData] Initializing uniform Robin BC values" << std::endl;
+    initialize_uniform(uniform_stiffness, uniform_damping, num_nodes);
+}
+
+void RobinBCData::load_from_vtp(const std::string& vtp_file_path, const faceType& face)
+{
+    //std::cout << "[load_from_vtp] Loading Robin BC from VTP file: " << vtp_file_path << std::endl;
+    // Check if file exists
+    if (FILE *file = fopen(vtp_file_path.c_str(), "r")) {
+        fclose(file);
+        //std::cout << "[load_from_vtp] File exists and is readable" << std::endl;
+    } else {
+        throw std::runtime_error("Robin BC VTP file '" + vtp_file_path + "' cannot be read.");
+    }
+    
+    // Read the VTP file
+    //std::cout << "[load_from_vtp] About to create VtkVtpData object" << std::endl;
+    VtkVtpData vtp_data(vtp_file_path, true);
+    //std::cout << "[load_from_vtp] VtkVtpData object created successfully" << std::endl;
+    
+    // Validate the VTP data matches the face geometry
+    validate_vtp_data(vtp_data, face);
+    
+    // Check for required point arrays
+    if (!vtp_data.has_point_data("Stiffness")) {
+        throw std::runtime_error("VTP file '" + vtp_file_path + "' does not contain 'Stiffness' point array.");
+    }
+    
+    if (!vtp_data.has_point_data("Damping")) {
+        throw std::runtime_error("VTP file '" + vtp_file_path + "' does not contain 'Damping' point array.");
+    }
+    
+    // Load stiffness array
+    auto stiffness_data = vtp_data.get_point_data("Stiffness");
+    if (stiffness_data.nrows() != num_nodes_ || stiffness_data.ncols() != 1) {
+        throw std::runtime_error("'Stiffness' array in VTP file '" + vtp_file_path + 
+                                "' has incorrect dimensions. Expected " + std::to_string(num_nodes_) + 
+                                " x 1, got " + std::to_string(stiffness_data.nrows()) + " x " + 
+                                std::to_string(stiffness_data.ncols()) + ".");
+    }
+    
+    // Load damping array
+    auto damping_data = vtp_data.get_point_data("Damping");
+    if (damping_data.nrows() != num_nodes_ || damping_data.ncols() != 1) {
+        throw std::runtime_error("'Damping' array in VTP file '" + vtp_file_path + 
+                                "' has incorrect dimensions. Expected " + std::to_string(num_nodes_) + 
+                                " x 1, got " + std::to_string(damping_data.nrows()) + " x " + 
+                                std::to_string(damping_data.ncols()) + ".");
+    }
+    
+    // Copy data to internal arrays and build global-to-local map
+    stiffness_array_.resize(num_nodes_);
+    damping_array_.resize(num_nodes_);
+    global_to_local_map_.clear();
+    
+    for (int i = 0; i < num_nodes_; i++) {
+        stiffness_array_(i) = stiffness_data(i, 0);
+        damping_array_(i) = damping_data(i, 0);
+        
+        // Store mapping from global node ID to local index
+        global_to_local_map_[face.gN(i)] = i;
+        
+        // Validate that values are non-negative
+        if (stiffness_array_(i) < 0.0) {
+            throw std::runtime_error("Negative stiffness value at node " + std::to_string(i) + 
+                                    " in VTP file '" + vtp_file_path + "'.");
+        }
+        if (damping_array_(i) < 0.0) {
+            throw std::runtime_error("Negative damping value at node " + std::to_string(i) + 
+                                    " in VTP file '" + vtp_file_path + "'.");
+        }
+    }
+}
+
+void RobinBCData::initialize_uniform(double stiffness, double damping, int num_nodes)
+{
+    if (stiffness < 0.0) {
+        throw std::runtime_error("Uniform stiffness value cannot be negative.");
+    }
+    if (damping < 0.0) {
+        throw std::runtime_error("Uniform damping value cannot be negative.");
+    }
+    if (num_nodes <= 0) {
+        throw std::runtime_error("Number of nodes must be positive.");
+    }
+    
+    stiffness_array_.resize(num_nodes);
+    damping_array_.resize(num_nodes);
+    global_to_local_map_.clear();
+    
+    for (int i = 0; i < num_nodes; i++) {
+        stiffness_array_(i) = stiffness;
+        damping_array_(i) = damping;
+        // For uniform values, we assume global IDs are sequential starting from 0
+        global_to_local_map_[i] = i;
+    }
+}
+
+void RobinBCData::validate_vtp_data(const VtkVtpData& vtp_data, const faceType& face)
+{
+    //std::cout << "[validate_vtp_data] Validating VTP data" << std::endl;
+    int vtp_num_points = vtp_data.num_points();
+    if (vtp_num_points != face.nNo) {
+        throw std::runtime_error("Number of points in VTP file (" + std::to_string(vtp_num_points) + 
+                                ") does not match number of nodes on face '" + face.name + 
+                                "' (" + std::to_string(face.nNo) + ").");
+    }
+    
+    int vtp_num_elems = vtp_data.num_elems();
+    if (vtp_num_elems != face.nEl) {
+        throw std::runtime_error("Number of elements in VTP file (" + std::to_string(vtp_num_elems) + 
+                                ") does not match number of elements on face '" + face.name + 
+                                "' (" + std::to_string(face.nEl) + ").");
+    }
+
+    // Get points from VTP file, (3,N) array
+    Array<double> vtp_points = vtp_data.get_points();
+    
+    // Compare each point's coordinates
+    const double tolerance = 1e-12;  // Small tolerance for floating point comparison
+    for (int i = 0; i < face.nNo; i++) {
+        double face_x = face.x(0,i);
+        double face_y = face.x(1,i);
+        double face_z = face.x(2,i);
+        double vtp_x = vtp_points(0,i);
+        double vtp_y = vtp_points(1,i);
+        double vtp_z = vtp_points(2,i);
+        double dx = std::abs(vtp_x - face_x);
+        double dy = std::abs(vtp_y - face_y);
+        double dz = std::abs(vtp_z - face_z);
+        
+        if (dx > tolerance || dy > tolerance || dz > tolerance) {
+            throw std::runtime_error("Node " + std::to_string(i) + " in VTP file is not colocated with " +
+                                   "corresponding node in face '" + face.name + "'. " +
+                                   "VTP: (" + std::to_string(vtp_x) + ", " + 
+                                   std::to_string(vtp_y) + ", " + 
+                                   std::to_string(vtp_z) + ") " +
+                                   "Face: (" + std::to_string(face_x) + ", " + 
+                                   std::to_string(face_y) + ", " + 
+                                   std::to_string(face_z) + ")");
+        }
+    }
+}
+
+double RobinBCData::get_stiffness(int node_id) const
+{
+    if (node_id < 0 || node_id >= num_nodes_) {
+        throw std::runtime_error("Node ID " + std::to_string(node_id) + 
+                                " is out of range [0, " + std::to_string(num_nodes_ - 1) + "].");
+    }
+    return stiffness_array_(node_id);
+}
+
+double RobinBCData::get_damping(int node_id) const
+{
+    if (node_id < 0 || node_id >= num_nodes_) {
+        throw std::runtime_error("Node ID " + std::to_string(node_id) + 
+                                " is out of range [0, " + std::to_string(num_nodes_ - 1) + "].");
+    }
+    return damping_array_(node_id);
+}
+
+int RobinBCData::get_local_index(int global_node_id) const
+{
+    auto it = global_to_local_map_.find(global_node_id);
+    if (it == global_to_local_map_.end()) {
+        throw std::runtime_error("Global node ID " + std::to_string(global_node_id) + 
+                               " not found in global-to-local map.");
+    }
+    return it->second;
+}

Code/Source/solver/RobinBCData.h

@@ -0,0 +1,86 @@
+#ifndef ROBIN_BC_DATA_H
+#define ROBIN_BC_DATA_H
+
+#include "Array.h"
+#include "Vector.h"
+#include "VtkData.h"
+#include "ComMod.h"
+#include <string>
+#include <memory>
+
+
+/// @brief Class to handle Robin boundary condition data with per-node stiffness and damping
+/// 
+/// This class reads stiffness and damping arrays from a VTP file and provides
+/// efficient access to per-node values during boundary condition assembly.
+class RobinBCData {
+public:
+    /// @brief Constructor - reads data from VTP file
+    /// @param vtp_file_path Path to VTP file containing Stiffness and Damping point arrays
+    /// @param face Reference to the face for validation
+    /// @throws std::runtime_error if file cannot be read or arrays are missing
+    RobinBCData(const std::string& vtp_file_path, const faceType& face);
+    
+    /// @brief Default constructor for uniform values
+    /// @param uniform_stiffness Uniform stiffness value for all nodes
+    /// @param uniform_damping Uniform damping value for all nodes
+    /// @param num_nodes Number of nodes on the face
+    RobinBCData(double uniform_stiffness, double uniform_damping, int num_nodes);
+    
+    /// @brief Destructor
+    ~RobinBCData() = default;
+    
+    /// @brief Get stiffness value for a specific node
+    /// @param node_id Node index on the face
+    /// @return Stiffness value for the node
+    double get_stiffness(int node_id) const;
+    
+    /// @brief Get damping value for a specific node
+    /// @param node_id Node index on the face
+    /// @return Damping value for the node
+    double get_damping(int node_id) const;
+    
+    /// @brief Get number of nodes
+    /// @return Number of nodes on the face
+    int get_num_nodes() const { return 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 stiffness_array_ and damping_array_
+    /// @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 from_vtp_; }
+    
+    /// @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_; }
+
+private:
+    /// @brief Load stiffness and damping arrays from VTP file
+    /// @param vtp_file_path Path to VTP file
+    /// @param face Reference to face for validation
+    void load_from_vtp(const std::string& vtp_file_path, const faceType& face);
+    
+    /// @brief Initialize uniform values
+    /// @param stiffness Uniform stiffness value
+    /// @param damping Uniform damping value
+    /// @param num_nodes Number of nodes
+    void initialize_uniform(double stiffness, double damping, int num_nodes);
+    
+    /// @brief Validate that VTP data matches face geometry
+    /// @param vtp_data VTP data object
+    /// @param face Reference to face
+    void validate_vtp_data(const VtkVtpData& vtp_data, const faceType& face);
+    
+    int num_nodes_;                    ///< Number of nodes on the face
+    Vector<double> stiffness_array_;   ///< Stiffness values for each node
+    Vector<double> damping_array_;     ///< Damping values for each node
+    bool from_vtp_;                    ///< 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_to_local_map_; ///< Maps global node IDs to local array indices
+};
+
+#endif // ROBIN_BC_DATA_H