diff --git a/Code/Source/solver/BoundaryCondition.cpp b/Code/Source/solver/BoundaryCondition.cpp new file mode 100644 index 000000000..736f40382 --- /dev/null +++ b/Code/Source/solver/BoundaryCondition.cpp @@ -0,0 +1,509 @@ +/* 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 "BoundaryCondition.h" +#include "ComMod.h" +#include "DebugMsg.h" +#include "Vector.h" +#include +#include +#include +#include +#include + +#define n_debug_bc + +BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std::vector& array_names, const StringBoolMap& flags, const faceType& face) + : 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) +{ + try { + global_data_ = read_data_from_vtp_file(vtp_file_path, array_names); + + // Validate values + for (const auto& [name, data] : global_data_) { + for (int i = 0; i < global_num_nodes_; i++) { + validate_array_value(name, data(i, 0)); + } + } + + // In case we are running sequentially, we need to fill the local arrays + // and the global node map as well, because distribute is not called in sequential mode. + local_data_ = global_data_; + + global_node_map_.clear(); + for (int i = 0; i < global_num_nodes_; i++) { + global_node_map_[face_->gN(i)] = i; + } + } catch (const std::exception& e) { + // Constructor failed - resources are automatically cleaned up + throw; // Re-throw the exception + } +} + +BoundaryCondition::BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face) + : face_(&face) + , global_num_nodes_(face.nNo) + , local_num_nodes_(0) + , spatially_variable(false) + , vtp_file_path_("") + , flags_(flags) +{ + try { + // Store array names, validate and store values + array_names_.clear(); + for (const auto& [name, value] : uniform_values) { + array_names_.push_back(name); + validate_array_value(name, value); + local_data_[name] = Array(1, 1); + local_data_[name](0, 0) = value; + } + } catch (const std::exception& e) { + // Constructor failed - resources are automatically cleaned up + throw; // Re-throw the exception + } +} + +BoundaryCondition::BoundaryCondition(const BoundaryCondition& other) + : face_(other.face_) + , global_num_nodes_(other.global_num_nodes_) + , local_num_nodes_(other.local_num_nodes_) + , array_names_(other.array_names_) + , local_data_(other.local_data_) + , global_data_(other.global_data_) + , spatially_variable(other.spatially_variable) + , vtp_file_path_(other.vtp_file_path_) + , flags_(other.flags_) + , global_node_map_(other.global_node_map_) +{ + if (other.vtp_data_) { + vtp_data_ = std::make_unique(*other.vtp_data_); + } +} + +void swap(BoundaryCondition& lhs, BoundaryCondition& rhs) noexcept { + using std::swap; + swap(lhs.face_, rhs.face_); + swap(lhs.global_num_nodes_, rhs.global_num_nodes_); + swap(lhs.local_num_nodes_, rhs.local_num_nodes_); + swap(lhs.array_names_, rhs.array_names_); + swap(lhs.local_data_, rhs.local_data_); + swap(lhs.global_data_, rhs.global_data_); + swap(lhs.spatially_variable, rhs.spatially_variable); + swap(lhs.vtp_file_path_, rhs.vtp_file_path_); + swap(lhs.flags_, rhs.flags_); + swap(lhs.global_node_map_, rhs.global_node_map_); + swap(lhs.vtp_data_, rhs.vtp_data_); +} + +BoundaryCondition& BoundaryCondition::operator=(BoundaryCondition other) { + swap(*this, other); + return *this; +} + +BoundaryCondition::BoundaryCondition(BoundaryCondition&& other) noexcept + : face_(other.face_) + , global_num_nodes_(other.global_num_nodes_) + , local_num_nodes_(other.local_num_nodes_) + , array_names_(std::move(other.array_names_)) + , local_data_(std::move(other.local_data_)) + , global_data_(std::move(other.global_data_)) + , spatially_variable(other.spatially_variable) + , vtp_file_path_(std::move(other.vtp_file_path_)) + , flags_(std::move(other.flags_)) + , global_node_map_(std::move(other.global_node_map_)) + , vtp_data_(std::move(other.vtp_data_)) +{ + other.face_ = nullptr; + other.global_num_nodes_ = 0; + other.local_num_nodes_ = 0; + other.spatially_variable = false; +} + + +BoundaryCondition::StringArrayMap BoundaryCondition::read_data_from_vtp_file(const std::string& vtp_file_path, const std::vector& array_names) +{ + #ifdef debug_bc + DebugMsg dmsg(__func__, 0); + dmsg << "Loading data from VTP file: " << vtp_file_path << std::endl; + dmsg << "Array names: " << array_names[0] << " and " << array_names[1] << std::endl; + #endif + + // Check if file exists + if (FILE *file = fopen(vtp_file_path.c_str(), "r")) { + fclose(file); + } else { + throw BoundaryConditionFileException(vtp_file_path); + } + + // Read the VTP file + try { + vtp_data_ = std::make_unique(vtp_file_path, true); + } catch (const std::exception& e) { + throw BoundaryConditionFileException(vtp_file_path); + } + + if (global_num_nodes_ != face_->nNo) { + throw BoundaryConditionNodeCountException(vtp_file_path, face_->name); + } + + // Read in the data from the VTP file + StringArrayMap result; + for (const auto& array_name : array_names) { + if (!vtp_data_->has_point_data(array_name)) { + throw BoundaryConditionVtpArrayException(vtp_file_path, array_name); + } + + auto array_data = vtp_data_->get_point_data(array_name); + + if (array_data.nrows() != global_num_nodes_ || array_data.ncols() != 1) { + throw BoundaryConditionVtpArrayDimensionException(vtp_file_path, array_name, + global_num_nodes_, 1, + array_data.nrows(), array_data.ncols()); + } + + // Store array in result map + result[array_name] = array_data; + + #ifdef debug_bc + dmsg << "Successfully loaded " << array_name << " data" << std::endl; + dmsg << array_name << " data size: " << array_data.nrows() << " x " << array_data.ncols() << std::endl; + #endif + } + + return result; +} + +double BoundaryCondition::get_value(const std::string& array_name, int node_id) const +{ + auto it = local_data_.find(array_name); + if (it == local_data_.end()) { + throw BoundaryConditionArrayException(array_name); + } + + if (node_id < 0 || node_id >= global_num_nodes_) { + throw BoundaryConditionNodeIdException(node_id, global_num_nodes_); + } + + // Return value + if (spatially_variable) { + return it->second(node_id, 0); + } else { + return it->second(0, 0); + } +} + +bool BoundaryCondition::get_flag(const std::string& name) const +{ + auto it = flags_.find(name); + if (it == flags_.end()) { + #ifdef debug_bc + DebugMsg dmsg(__func__, 0); + dmsg << "Flag '" << name << "' not found. Available flags: " << flags_to_string() << std::endl; + #endif + throw BoundaryConditionFlagException(name); + } + return it->second; +} + +int BoundaryCondition::get_local_index(int global_node_id) const +{ + if (spatially_variable) { + auto it = global_node_map_.find(global_node_id); + if (it == global_node_map_.end()) { + throw BoundaryConditionGlobalNodeIdException(global_node_id); + } + return it->second; + } else { + return 0; + } +} + +void BoundaryCondition::distribute(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const faceType& face) +{ + #define n_debug_distribute + #ifdef debug_distribute + DebugMsg dmsg(__func__, cm.idcm()); + dmsg << "Distributing BC data" << std::endl; + #endif + + // Before this point, only the master process had a face_ pointer, which contained + // the entire (global) face. Within the distribute::distribute() function, + // the global face was partitioned and distributed among all processes. Each + // local face contains only a portion of the global face, corresponding to + // the portion of the volume mesh owned by this process. Here, we update the + // face_ pointer to the local face. + face_ = &face; + + // Number of nodes on the face on this processor + local_num_nodes_ = face_->nNo; + const bool is_slave = cm.slv(cm_mod); + distribute_metadata(cm_mod, cm, is_slave); + if (spatially_variable) { + distribute_spatially_variable(com_mod, cm_mod, cm, is_slave); + } else { + distribute_uniform(cm_mod, cm, is_slave); + } + distribute_flags(cm_mod, cm, is_slave); + + #ifdef debug_distribute + dmsg << "Finished distributing BC data" << std::endl; + dmsg << "Number of face nodes on this processor: " << local_num_nodes_ << std::endl; + #endif +} + +void BoundaryCondition::distribute_metadata(const CmMod& cm_mod, const cmType& cm, bool is_slave) +{ + cm.bcast(cm_mod, &spatially_variable); + + // Not necessary, but we do it for consistency + if (spatially_variable) { + cm.bcast(cm_mod, vtp_file_path_); + } + + // Not necessary, but we do it for consistency + cm.bcast(cm_mod, &global_num_nodes_); + + // Communicate array names + int num_arrays = array_names_.size(); + cm.bcast(cm_mod, &num_arrays); + if (is_slave) { + array_names_.resize(num_arrays); + } + for (int i = 0; i < num_arrays; i++) { + if (!is_slave) { + std::string& array_name = array_names_[i]; + cm.bcast(cm_mod, array_name); + } else { + std::string array_name; + cm.bcast(cm_mod, array_name); + array_names_[i] = array_name; + } + } +} + +void BoundaryCondition::distribute_spatially_variable(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bool is_slave) +{ + #define n_debug_distribute_spatially_variable + #ifdef debug_distribute_spatially_variable + DebugMsg dmsg(__func__, 0); + #endif + + if (face_ == nullptr) { + throw BoundaryConditionNullFaceException(); + } + // Each processor collects the global node IDs and nodal positions of its + // associated face portion + Vector local_global_ids = face_->gN; + Array local_positions(3, local_num_nodes_); + for (int i = 0; i < local_num_nodes_; i++) { + local_positions.set_col(i, com_mod.x.col(face_->gN(i))); + } + + #ifdef debug_distribute_spatially_variable + dmsg << "Number of face nodes on this processor: " << local_num_nodes_ << std::endl; + dmsg << "Local global IDs: " << local_global_ids << std::endl; + dmsg << "Local positions: " << local_positions << std::endl; + #endif + // Gather number of face nodes from each processor to master + Vector proc_num_nodes(cm.np()); + cm.gather(cm_mod, &local_num_nodes_, 1, proc_num_nodes.data(), 1, 0); + + // Calculate displacements for gatherv/scatterv and compute total number of nodes + // total_num_nodes is the total number of face nodes across all processors. + Vector displs(cm.np()); + int total_num_nodes = 0; + for (int i = 0; i < cm.np(); i++) { + displs(i) = total_num_nodes; + total_num_nodes += proc_num_nodes(i); + } + + // Master process: gather the nodal positions of face nodes from all processors, + // get the corresponding array values by matching the positions to the VTP points, + // and scatter the data back to all processors. + Array all_positions; + std::map> all_values; + if (!is_slave) { + // Resize receive buffers based on total number of nodes + all_positions.resize(3, total_num_nodes); + + // Gather all positions to master using gatherv + for (int d = 0; d < 3; d++) { + Vector local_pos_d(local_num_nodes_); + Vector all_pos_d(total_num_nodes); + for (int i = 0; i < local_num_nodes_; i++) { + local_pos_d(i) = local_positions(d,i); + } + cm.gatherv(cm_mod, local_pos_d, all_pos_d, proc_num_nodes, displs, 0); + for (int i = 0; i < total_num_nodes; i++) { + all_positions(d,i) = all_pos_d(i); + } + } + + // Get VTP points for position matching + Array vtp_points = vtp_data_->get_points(); + + // 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); + all_values[array_name](i) = global_data_[array_name](vtp_idx, 0); + } + } + + // Clear global data to save memory + global_data_.clear(); + } else { + // Slave processes: send node positions to master + for (int d = 0; d < 3; d++) { + Vector local_pos_d(local_num_nodes_); + for (int i = 0; i < local_num_nodes_; i++) { + local_pos_d(i) = local_positions(d,i); + } + Vector dummy_recv(total_num_nodes); + cm.gatherv(cm_mod, local_pos_d, dummy_recv, proc_num_nodes, displs, 0); + } + } + + // Scatter data back to all processes using scatterv + local_data_.clear(); + for (const auto& array_name : array_names_) { + Vector local_values(local_num_nodes_); + cm.scatterv(cm_mod, all_values[array_name], proc_num_nodes, displs, local_values, 0); + local_data_[array_name] = Array(local_num_nodes_, 1); + local_data_[array_name].set_col(0, local_values); + } + + // Build mapping from face global node IDs to local array indices so we can + // get data from a global node ID + global_node_map_.clear(); + for (int i = 0; i < local_num_nodes_; i++) { + global_node_map_[local_global_ids(i)] = i; + } + + #ifdef debug_distribute_spatially_variable + dmsg << "Checking if local arrays and node positions are consistent" << std::endl; + for (int i = 0; i < local_num_nodes_; i++) { + dmsg << "Local global ID: " << local_global_ids(i) << std::endl; + dmsg << "Local index: " << get_local_index(local_global_ids(i)) << std::endl; + dmsg << "Local position: " << com_mod.x.col(local_global_ids(i)) << std::endl; + for (const auto& array_name : array_names_) { + dmsg << "Local " << array_name << ": " << local_data_[array_name](i, 0) << std::endl; + } + } + #endif +} + +void BoundaryCondition::distribute_uniform(const CmMod& cm_mod, const cmType& cm, bool is_slave) +{ + if (!is_slave) { + for (const auto& array_name : array_names_) { + double uniform_value = local_data_[array_name](0, 0); + cm.bcast(cm_mod, &uniform_value); + } + } else { + local_data_.clear(); + for (const auto& array_name : array_names_) { + double uniform_value; + cm.bcast(cm_mod, &uniform_value); + local_data_[array_name] = Array(1, 1); + local_data_[array_name](0, 0) = uniform_value; + } + } +} + +void BoundaryCondition::distribute_flags(const CmMod& cm_mod, const cmType& cm, bool is_slave) +{ + if (cm.seq()) return; + int num_flags = 0; + if (!is_slave) { + num_flags = static_cast(flags_.size()); + } + cm.bcast(cm_mod, &num_flags); + if (is_slave) { + flags_.clear(); + } + for (int i = 0; i < num_flags; i++) { + std::string key; + bool val = false; + if (!is_slave) { + auto it = std::next(flags_.begin(), i); + key = it->first; + val = it->second; + cm.bcast(cm_mod, key); + cm.bcast(cm_mod, &val); + } else { + cm.bcast(cm_mod, key); + cm.bcast(cm_mod, &val); + flags_[key] = val; + } + } +} + +int BoundaryCondition::find_vtp_point_index(double x, double y, double z, + const Array& vtp_points) const +{ + const int num_points = vtp_points.ncols(); + Vector target_point{x, y, z}; + + // Simple linear search through all points in the VTP file + for (int i = 0; i < num_points; i++) { + auto vtp_point = vtp_points.col(i); + auto diff = vtp_point - target_point; + double distance = sqrt(diff.dot(diff)); + + if (distance <= POINT_MATCH_TOLERANCE) { + #define n_debug_bc_find_vtp_point_index + #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 << "VTP point index: " << i << std::endl; + #endif + + return i; + } + } + + throw BoundaryConditionPointNotFoundException(x, y, z); +} + +std::string BoundaryCondition::flags_to_string() const { + std::string result; + for (const auto& [name, value] : flags_) { + result += name + ": " + (value ? "true" : "false") + ", "; + } + return result; +} \ No newline at end of file diff --git a/Code/Source/solver/BoundaryCondition.h b/Code/Source/solver/BoundaryCondition.h new file mode 100644 index 000000000..7105c8e0e --- /dev/null +++ b/Code/Source/solver/BoundaryCondition.h @@ -0,0 +1,340 @@ +/* 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 +#include +#include +#include +#include +#include + +// 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& array_names, const faceType& face) +/// : BoundaryCondition(vtp_file_path, array_names, face) {} +/// +/// MyBoundaryCondition(const std::map& 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>; + using StringBoolMap = std::map; + using StringDoubleMap = std::map; + + /// @brief Data members for BC + const faceType* face_ = nullptr; ///< Face associated with the BC (can be null) + 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 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 = false; ///< Flag indicating if data is from VTP file + std::string vtp_file_path_; ///< Path to VTP file (empty if uniform) + std::map global_node_map_; ///< Maps global node IDs to local array indices + std::unique_ptr vtp_data_; ///< VTP data object + +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() = default; + + /// @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& 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 Unified assignment operator (handles both copy and move) + BoundaryCondition& operator=(BoundaryCondition other); + + /// @brief Move constructor + BoundaryCondition(BoundaryCondition&& other) noexcept; + + /// @brief Virtual destructor + virtual ~BoundaryCondition() noexcept = default; + + /// @brief Swap function for copy-and-swap idiom (friend function) + friend void swap(BoundaryCondition& lhs, BoundaryCondition& rhs) noexcept; + + + /// @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 + /// @throws std::runtime_error if array_name is not found + 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 noexcept { + 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 noexcept { + 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 noexcept { + 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 noexcept { + return vtp_file_path_; + } + + /// @brief Check if this BC has been properly initialized with data + /// @return true if BC has data (either global or local arrays are populated) + bool is_initialized() const noexcept { + return !global_data_.empty() || !local_data_.empty(); + } + + /// @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& 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& 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 +/// +/// These exceptions indicate fatal errors that should terminate the solver. +/// They should not be caught and handled gracefully - the solver should exit. +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 + "'") {} +}; + +/// @brief Exception thrown when a requested array is not found +class BoundaryConditionArrayException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionArrayException(const std::string& array_name) + : BoundaryConditionBaseException("BoundaryCondition array not found: '" + array_name + "'") {} +}; + +/// @brief Exception thrown when BoundaryCondition is not properly initialized +class BoundaryConditionNotInitializedException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionNotInitializedException() + : BoundaryConditionBaseException("BoundaryCondition not properly initialized - no data available") {} +}; + +/// @brief Exception thrown when a node ID is out of range +class BoundaryConditionNodeIdException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionNodeIdException(int node_id, int max_node_id) + : BoundaryConditionBaseException("Node ID " + std::to_string(node_id) + + " is out of range [0, " + std::to_string(max_node_id - 1) + "]") {} +}; + +/// @brief Exception thrown when a global node ID is not found in the global-to-local map +class BoundaryConditionGlobalNodeIdException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionGlobalNodeIdException(int global_node_id) + : BoundaryConditionBaseException("Global node ID " + std::to_string(global_node_id) + + " not found in global-to-local map") {} +}; + +/// @brief Exception thrown when a VTP file doesn't contain a required array +class BoundaryConditionVtpArrayException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionVtpArrayException(const std::string& vtp_file, const std::string& array_name) + : BoundaryConditionBaseException("VTP file '" + vtp_file + "' does not contain '" + array_name + "' point array") {} +}; + +/// @brief Exception thrown when a VTP array has incorrect dimensions +class BoundaryConditionVtpArrayDimensionException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionVtpArrayDimensionException(const std::string& vtp_file, const std::string& array_name, + int expected_rows, int expected_cols, int actual_rows, int actual_cols) + : BoundaryConditionBaseException("'" + array_name + "' array in VTP file '" + vtp_file + + "' has incorrect dimensions. Expected " + std::to_string(expected_rows) + + " x " + std::to_string(expected_cols) + ", got " + std::to_string(actual_rows) + + " x " + std::to_string(actual_cols)) {} +}; + +/// @brief Exception thrown when face_ is nullptr during distribute +class BoundaryConditionNullFaceException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionNullFaceException() + : BoundaryConditionBaseException("face_ is nullptr during distribute") {} +}; + +/// @brief Exception thrown when a point cannot be found in VTP file +class BoundaryConditionPointNotFoundException : public BoundaryConditionBaseException { +public: + explicit BoundaryConditionPointNotFoundException(double x, double y, double z) + : BoundaryConditionBaseException("Could not find matching point in VTP file for node at position (" + + std::to_string(x) + ", " + std::to_string(y) + ", " + std::to_string(z) + ")") {} +}; + +#endif // BOUNDARY_CONDITION_H diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 85da4dd30..e30adb391 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/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. diff --git a/Code/Source/solver/CmMod.cpp b/Code/Source/solver/CmMod.cpp index d5a97b1a3..25052be32 100644 --- a/Code/Source/solver/CmMod.cpp +++ b/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& 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& send_data, Vector& 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& send_data, Vector& 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& send_data, Vector& 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& send_data, Vector& 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& send_data, Vector& recv_data, const Vector& recv_counts, const Vector& 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& send_data, Vector& recv_data, const Vector& recv_counts, const Vector& 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& send_data, const Vector& send_counts, const Vector& displs, Vector& 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& send_data, const Vector& send_counts, const Vector& displs, Vector& 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()); } \ No newline at end of file diff --git a/Code/Source/solver/CmMod.h b/Code/Source/solver/CmMod.h index 1c00bd082..b04942101 100644 --- a/Code/Source/solver/CmMod.h +++ b/Code/Source/solver/CmMod.h @@ -112,6 +112,26 @@ class cmType { void bcast(const CmMod& cm_mod, Array& 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& send_data, Vector& recv_data, int root) const; + void gather(const CmMod& cm_mod, const Vector& send_data, Vector& recv_data, int root) const; + + // Gatherv operations + void gatherv(const CmMod& cm_mod, const Vector& send_data, Vector& recv_data, const Vector& recv_counts, const Vector& displs, int root) const; + void gatherv(const CmMod& cm_mod, const Vector& send_data, Vector& recv_data, const Vector& recv_counts, const Vector& displs, int root) const; + + // Scatterv operations + void scatterv(const CmMod& cm_mod, const Vector& send_data, const Vector& send_counts, const Vector& displs, Vector& recv_data, int root) const; + void scatterv(const CmMod& cm_mod, const Vector& send_data, const Vector& send_counts, const Vector& displs, Vector& 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& send_data, Vector& recv_data, int root) const; + void scatter(const CmMod& cm_mod, const Vector& send_data, Vector& recv_data, int root) const; + //------------ // bcast_enum //------------ diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index da2ff5f90..d44835a35 100644 --- a/Code/Source/solver/ComMod.h +++ b/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 #include +#include #include #include #include @@ -150,7 +152,6 @@ class rcrType class bcType { public: - // Strong/Weak application of Dirichlet BC bool weakDir = false; @@ -158,9 +159,6 @@ class bcType // (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. diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index a818ab9c6..02cffa0a4 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/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); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b70b72246..e9fd424d1 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -806,6 +806,7 @@ class BoundaryConditionParameters : public ParameterLists Parameter spatial_profile_file_path; Parameter spatial_values_file_path; Parameter stiffness; + Parameter robin_vtp_file_path; Parameter svzerod_solver_block; Parameter temporal_and_spatial_values_file_path; diff --git a/Code/Source/solver/RobinBoundaryCondition.cpp b/Code/Source/solver/RobinBoundaryCondition.cpp new file mode 100644 index 000000000..dd16a9f2b --- /dev/null +++ b/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 + diff --git a/Code/Source/solver/RobinBoundaryCondition.h b/Code/Source/solver/RobinBoundaryCondition.h new file mode 100644 index 000000000..13b53241c --- /dev/null +++ b/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 +#include +#include + +/// @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 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 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{"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 diff --git a/Code/Source/solver/VtkData.cpp b/Code/Source/solver/VtkData.cpp index 207306f0a..29ba340ec 100644 --- a/Code/Source/solver/VtkData.cpp +++ b/Code/Source/solver/VtkData.cpp @@ -30,6 +30,7 @@ #include "VtkData.h" #include "Array.h" +#include "DebugMsg.h" #include #include "vtkCellArray.h" @@ -59,6 +60,8 @@ class VtkVtpData::VtkVtpDataImpl { public: VtkVtpDataImpl(); + VtkVtpDataImpl(const VtkVtpDataImpl& other); + VtkVtpDataImpl& operator=(const VtkVtpDataImpl& other); void read_file(const std::string& file_name); void set_connectivity(const int nsd, const Array& conn, const int pid); void set_point_data(const std::string& data_name, const Vector& data); @@ -77,6 +80,29 @@ VtkVtpData::VtkVtpDataImpl::VtkVtpDataImpl() vtk_polydata = vtkSmartPointer::New(); } +VtkVtpData::VtkVtpDataImpl::VtkVtpDataImpl(const VtkVtpDataImpl& other) + : elem_type(other.elem_type) + , num_elems(other.num_elems) + , np_elem(other.np_elem) + , num_points(other.num_points) +{ + vtk_polydata = vtkSmartPointer::New(); + vtk_polydata->DeepCopy(other.vtk_polydata); +} + +VtkVtpData::VtkVtpDataImpl& VtkVtpData::VtkVtpDataImpl::operator=(const VtkVtpDataImpl& other) +{ + if (this != &other) { + elem_type = other.elem_type; + num_elems = other.num_elems; + np_elem = other.np_elem; + num_points = other.num_points; + vtk_polydata = vtkSmartPointer::New(); + vtk_polydata->DeepCopy(other.vtk_polydata); + } + return *this; +} + void VtkVtpData::VtkVtpDataImpl::read_file(const std::string& file_name) { auto reader = vtkSmartPointer::New(); @@ -98,14 +124,21 @@ void VtkVtpData::VtkVtpDataImpl::read_file(const std::string& file_name) void VtkVtpData::VtkVtpDataImpl::set_connectivity(const int nsd, const Array& conn, const int pid) { - //std::cout << "[VtkVtpData.set_connectivity] " << std::endl; - //std::cout << "[VtkVtpData.set_connectivity] vtk_polydata: " << vtk_polydata << std::endl; - //std::cout << "[VtkVtpData.set_connectivity] nsd: " << nsd << std::endl; + #define n_debug_vtk_data + #ifdef debug_vtk_data + DebugMsg dmsg(__func__, 0); + dmsg << "[VtkVtpData.set_connectivity] vtk_polydata: " << vtk_polydata; + dmsg << "[VtkVtpData.set_connectivity] nsd: " << nsd; + #endif + int num_elems = conn.ncols(); int np_elem = conn.nrows(); unsigned char vtk_cell_type; - //std::cout << "[VtkVtpData.set_connectivity] num_elems: " << num_elems << std::endl; - //std::cout << "[VtkVtpData.set_connectivity] np_elem: " << np_elem << std::endl; + + #ifdef debug_vtk_data + dmsg << "[VtkVtpData.set_connectivity] num_elems: " << num_elems; + dmsg << "[VtkVtpData.set_connectivity] np_elem: " << np_elem; + #endif if (nsd == 2) { if (np_elem == 4) { @@ -130,7 +163,9 @@ void VtkVtpData::VtkVtpDataImpl::set_connectivity(const int nsd, const Array element_cells = vtkSmartPointer::New(); for (int i = 0; i < num_elems; i++) { - //std::cout << "[VtkVtpData.set_connectivity] ---------- i " << i << std::endl; + #ifdef debug_vtk_data + dmsg << "[VtkVtpData.set_connectivity] ---------- i " << i; + #endif for (int j = 0; j < np_elem; j++) { - //std::cout << "[VtkVtpData.set_connectivity] ----- j " << j << std::endl; - //std::cout << "[VtkVtpData.set_connectivity] conn(j,i): " << conn(j,i) << std::endl; + #ifdef debug_vtk_data + dmsg << "[VtkVtpData.set_connectivity] ----- j " << j; + dmsg << "[VtkVtpData.set_connectivity] conn(j,i): " << conn(j,i); + #endif elem_nodes->SetId(j, conn(j,i)); } element_cells->InsertNextCell(elem_nodes); @@ -201,12 +240,19 @@ void VtkVtpData::VtkVtpDataImpl::set_point_data(const std::string& data_name, co // void VtkVtpData::VtkVtpDataImpl::set_points(const Array& points) { - //std::cout << "[VtkVtpData.set_points] vtk_polydata: " << vtk_polydata << std::endl; + #ifdef debug_vtk_data + DebugMsg dmsg(__func__, 0); + dmsg << "[VtkVtpData.set_points] vtk_polydata: " << vtk_polydata; + #endif + int num_coords = points.ncols(); if (num_coords == 0) { throw std::runtime_error("Error in VTK VTP set_points: the number of points is zero."); } - //std::cout << "[VtkVtpData.set_points] num_coords: " << num_coords << std::endl; + + #ifdef debug_vtk_data + dmsg << "[VtkVtpData.set_points] num_coords: " << num_coords; + #endif auto node_coords = vtkSmartPointer::New(); node_coords->Allocate(num_coords, 1000); @@ -229,7 +275,10 @@ void VtkVtpData::VtkVtpDataImpl::set_points(const Array& points) void VtkVtpData::VtkVtpDataImpl::write(const std::string& file_name) { - //std::cout << "[VtkVtpData.write] file_name: " << file_name << std::endl; + #ifdef debug_vtk_data + DebugMsg dmsg(__func__, 0); + dmsg << "[VtkVtpData.write] file_name: " << file_name; + #endif auto writer = vtkSmartPointer::New(); writer->SetInputDataObject(vtk_polydata); writer->SetFileName(file_name.c_str()); @@ -572,7 +621,23 @@ VtkVtpData::~VtkVtpData() delete impl; } -Array VtkVtpData::get_connectivity() +VtkVtpData::VtkVtpData(const VtkVtpData& other) +{ + impl = new VtkVtpDataImpl(*other.impl); + file_name = other.file_name; +} + +VtkVtpData& VtkVtpData::operator=(const VtkVtpData& other) +{ + if (this != &other) { + delete impl; + impl = new VtkVtpDataImpl(*other.impl); + file_name = other.file_name; + } + return *this; +} + +Array VtkVtpData::get_connectivity() const { int num_elems = impl->num_elems; int np_elem = impl->np_elem; @@ -719,7 +784,7 @@ std::vector VtkVtpData::get_point_data_names() /// @brief Get an array of point data from an unstructured grid. // -Array VtkVtpData::get_points() +Array VtkVtpData::get_points() const { auto vtk_points = impl->vtk_polydata->GetPoints(); auto num_points = vtk_points->GetNumberOfPoints(); @@ -749,22 +814,22 @@ bool VtkVtpData::has_point_data(const std::string& data_name) return false; } -int VtkVtpData::elem_type() +int VtkVtpData::elem_type() const { return impl->elem_type; } -int VtkVtpData::num_elems() +int VtkVtpData::num_elems() const { return impl->num_elems; } -int VtkVtpData::np_elem() +int VtkVtpData::np_elem() const { return impl->np_elem; } -int VtkVtpData::num_points() +int VtkVtpData::num_points() const { return impl->num_points; } @@ -840,7 +905,7 @@ VtkVtuData::~VtkVtuData() delete impl; } -Array VtkVtuData::get_connectivity() +Array VtkVtuData::get_connectivity() const { int num_elems = impl->num_elems; int np_elem = impl->np_elem; @@ -999,7 +1064,7 @@ Array VtkVtuData::get_point_data(const std::string& data_name) return data; } -Array VtkVtuData::get_points() +Array VtkVtuData::get_points() const { auto vtk_points = impl->vtk_ugrid->GetPoints(); auto num_points = vtk_points->GetNumberOfPoints(); @@ -1016,22 +1081,22 @@ Array VtkVtuData::get_points() return points_array; } -int VtkVtuData::elem_type() +int VtkVtuData::elem_type() const { return impl->elem_type; } -int VtkVtuData::num_elems() +int VtkVtuData::num_elems() const { return impl->num_elems; } -int VtkVtuData::np_elem() +int VtkVtuData::np_elem() const { return impl->np_elem; } -int VtkVtuData::num_points() +int VtkVtuData::num_points() const { return impl->num_points; } diff --git a/Code/Source/solver/VtkData.h b/Code/Source/solver/VtkData.h index 41a240948..c8c0a956d 100644 --- a/Code/Source/solver/VtkData.h +++ b/Code/Source/solver/VtkData.h @@ -41,12 +41,12 @@ class VtkData { VtkData(); virtual ~VtkData(); - virtual Array get_connectivity() = 0; - virtual Array get_points() = 0; - virtual int num_elems() = 0; - virtual int elem_type() = 0; - virtual int np_elem() = 0; - virtual int num_points() = 0; + virtual Array get_connectivity() const = 0; + virtual Array get_points() const = 0; + virtual int num_elems() const = 0; + virtual int elem_type() const = 0; + virtual int np_elem() const = 0; + virtual int num_points() const = 0; virtual void read_file(const std::string& file_name) = 0; virtual void set_element_data(const std::string& data_name, const Array& data) = 0; @@ -78,32 +78,37 @@ class VtkVtpData : public VtkData { VtkVtpData(const std::string& file_name, bool reader=true); ~VtkVtpData(); - virtual Array get_connectivity(); - virtual Array get_points(); - virtual int elem_type(); - virtual int num_elems(); - virtual int np_elem(); - virtual int num_points(); - virtual void read_file(const std::string& file_name); - - void copy_points(Array& points); - void copy_point_data(const std::string& data_name, Array& mesh_data); - void copy_point_data(const std::string& data_name, Vector& mesh_data); + // Copy constructor + VtkVtpData(const VtkVtpData& other); + // Copy assignment operator + VtkVtpData& operator=(const VtkVtpData& other); + + virtual Array get_connectivity() const override; + virtual Array get_points() const override; + virtual int elem_type() const override; + virtual int num_elems() const override; + virtual int np_elem() const override; + virtual int num_points() const override; + virtual void read_file(const std::string& file_name) override; + + void copy_points(Array& points) override; + void copy_point_data(const std::string& data_name, Array& mesh_data) override; + void copy_point_data(const std::string& data_name, Vector& mesh_data) override; void copy_point_data(const std::string& data_name, Vector& mesh_data); Array get_point_data(const std::string& data_name); std::vector get_point_data_names(); - bool has_point_data(const std::string& data_name); - virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0); + bool has_point_data(const std::string& data_name) override; + virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0) override; - virtual void set_element_data(const std::string& data_name, const Array& data); - virtual void set_element_data(const std::string& data_name, const Array& data); + virtual void set_element_data(const std::string& data_name, const Array& data) override; + virtual void set_element_data(const std::string& data_name, const Array& data) override; - virtual void set_point_data(const std::string& data_name, const Array& data); - virtual void set_point_data(const std::string& data_name, const Array& data); - virtual void set_point_data(const std::string& data_name, const Vector& data); + virtual void set_point_data(const std::string& data_name, const Array& data) override; + virtual void set_point_data(const std::string& data_name, const Array& data) override; + virtual void set_point_data(const std::string& data_name, const Vector& data) override; - virtual void set_points(const Array& points); - virtual void write(); + virtual void set_points(const Array& points) override; + virtual void write() override; private: class VtkVtpDataImpl; @@ -116,33 +121,33 @@ class VtkVtuData : public VtkData { VtkVtuData(const std::string& file_name, bool reader=true); ~VtkVtuData(); - virtual Array get_connectivity(); - virtual int elem_type(); - virtual int num_elems(); - virtual int np_elem(); - virtual int num_points(); - virtual void read_file(const std::string& file_name); + virtual Array get_connectivity() const override; + virtual int elem_type() const override; + virtual int num_elems() const override; + virtual int np_elem() const override; + virtual int num_points() const override; + virtual void read_file(const std::string& file_name) override; - void copy_points(Array& points); - void copy_point_data(const std::string& data_name, Array& mesh_data); - void copy_point_data(const std::string& data_name, Vector& mesh_data); + void copy_points(Array& points) override; + void copy_point_data(const std::string& data_name, Array& mesh_data) override; + void copy_point_data(const std::string& data_name, Vector& mesh_data) override; void copy_point_data(const std::string& data_name, Vector& mesh_data); Array get_point_data(const std::string& data_name); std::vector get_point_data_names(); - virtual Array get_points(); - bool has_point_data(const std::string& data_name); - virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0); + virtual Array get_points() const override; + bool has_point_data(const std::string& data_name) override; + virtual void set_connectivity(const int nsd, const Array& conn, const int pid = 0) override; - virtual void set_element_data(const std::string& data_name, const Array& data); - virtual void set_element_data(const std::string& data_name, const Array& data); + virtual void set_element_data(const std::string& data_name, const Array& data) override; + virtual void set_element_data(const std::string& data_name, const Array& data) override; - virtual void set_point_data(const std::string& data_name, const Array& data); - virtual void set_point_data(const std::string& data_name, const Array& data); - virtual void set_point_data(const std::string& data_name, const Vector& data); + virtual void set_point_data(const std::string& data_name, const Array& data) override; + virtual void set_point_data(const std::string& data_name, const Array& data) override; + virtual void set_point_data(const std::string& data_name, const Vector& data) override; - virtual void set_points(const Array& points); - virtual void write(); + virtual void set_points(const Array& points) override; + virtual void write() override; private: class VtkVtuDataImpl; diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index 55ce241ef..65e65018a 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -274,8 +274,10 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l return; } - if (btest(lBc.bType, enum_int(BoundaryConditionType::bType_Robin)) && not com_mod.dFlag) { - throw std::runtime_error("Robin BC can be set for a displacement-based eqn only"); + if (btest(lBc.bType, enum_int(BoundaryConditionType::bType_Robin))) { + if (!com_mod.dFlag) { + throw std::runtime_error("Robin BC can be set for a displacement-based eqn only"); + } } int iM = lFa.iM; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index a503f2b2d..3643f474e 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -33,6 +33,7 @@ #include "distribute.h" #include "all_fun.h" +#include "ComMod.h" #include "consts.h" #include "nn.h" #include "utils.h" @@ -654,13 +655,10 @@ void dist_bc(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bcType& lBc cm.bcast(cm_mod, &lBc.iM); cm.bcast(cm_mod, &lBc.r); cm.bcast(cm_mod, &lBc.g); - cm.bcast(cm_mod, &lBc.k); - cm.bcast(cm_mod, &lBc.c); cm.bcast(cm_mod, lBc.h); cm.bcast(cm_mod, &lBc.weakDir); cm.bcast(cm_mod, lBc.tauB); cm.bcast(cm_mod, &lBc.flwP); - cm.bcast(cm_mod, &lBc.rbnN); if (utils::btest(lBc.bType, static_cast(BoundaryConditionType::bType_RCR))) { @@ -802,6 +800,16 @@ void dist_bc(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bcType& lBc } } + // Communicating Robin BC + // + bool has_robin_bc = lBc.robin_bc.is_initialized(); + cm.bcast(cm_mod, &has_robin_bc); // Master process broadcasts the flag to all processes + + if (has_robin_bc) { + lBc.robin_bc.distribute(com_mod, cm_mod, cm, com_mod.msh[lBc.iM].fa[lBc.iFa]); + } + + // Communicating and reordering master node data for // undeforming Neumann BC faces. // @@ -1795,8 +1803,9 @@ void dist_solid_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmT void part_face(Simulation* simulation, mshType& lM, faceType& lFa, faceType& gFa, Vector& gmtl) { + #define n_debug_part_face #ifdef debug_part_face - DebugMsg dmsg(__func__, com_mod.cm.idcm()); + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); dmsg.banner(); dmsg << "lFa.name: " << lFa.name; #endif @@ -1880,6 +1889,7 @@ void part_face(Simulation* simulation, mshType& lM, faceType& lFa, faceType& gFa int i = gFa.nEl*(2+eNoNb) + gFa.nNo; Vector part(i); + // Broadcast data in gFa from master to all processes if (cm.mas(cm_mod)) { for (int e = 0; e < gFa.nEl; e++) { int Ec = gFa.gE[e]; diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index a87fb3789..8ec89f3b2 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -392,9 +392,18 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, // Stiffness and damping parameters for Robin BC if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Robin))) { - lBc.k = bc_params->stiffness.value(); - lBc.c = bc_params->damping.value(); - lBc.rbnN = bc_params->apply_along_normal_direction.value(); + + // 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(), + bc_params->apply_along_normal_direction.value(), + com_mod.msh[lBc.iM].fa[lBc.iFa]); + } 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]); + } } // To impose value or flux diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 266b99608..4613efb53 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1470,19 +1470,24 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const } } // Now treat Robin BC (stiffness and damping) here - // - if (utils::btest(lBc.bType,iBC_Robin)) { - set_bc_rbnl(com_mod, lFa, lBc.k, lBc.c, lBc.rbnN, Yg, Dg); + if (lBc.robin_bc.is_initialized()) { + set_bc_rbnl(com_mod, lFa, lBc.robin_bc, Yg, Dg); } } /// @brief Set Robin BC contribution to residual and tangent // -void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const double cs, const bool isN, +void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondition& robin_bc, const Array& Yg, const Array& Dg) { using namespace consts; + #define n_debug_set_bc_rbnl_l + #ifdef debug_set_bc_rbnl_l + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + const int cEq = com_mod.cEq; const auto& eq = com_mod.eq[cEq]; const double dt = com_mod.dt; @@ -1544,9 +1549,30 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const do auto nDn = mat_fun::mat_id(nsd); Vector h; - h = ks*u + cs*ud; + + #ifdef debug_set_bc_rbnl_l + dmsg << "Calculating weighted average of stiffness and damping for this integration point"; + #endif + // Calculate weighted average of stiffness and damping for this integration point + double ks_avg = 0.0; + double cs_avg = 0.0; + for (int a = 0; a < eNoN; a++) { + int Ac = lFa.IEN(a,e); + #ifdef debug_set_bc_rbnl_l + dmsg << "e: " << e << std::endl; + dmsg << "a: " << a << std::endl; + dmsg << "Local position: " << com_mod.x.col(Ac) << std::endl; + dmsg << "Stiffness: " << robin_bc.get_stiffness(robin_bc.get_local_index(Ac)) << std::endl; + #endif + ks_avg += N(a) * robin_bc.get_stiffness(robin_bc.get_local_index(Ac)); + cs_avg += N(a) * robin_bc.get_damping(robin_bc.get_local_index(Ac)); + } + + + + h = ks_avg*u + cs_avg*ud; - if (isN) { + if (robin_bc.normal_direction_only()) { h = utils::norm(h, nV) * nV; for (int a = 0; a < nsd; a++) { for (int b = 0; b < nsd; b++) { @@ -1568,8 +1594,11 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const do for (int a = 0; a < eNoN; a++) { for (int b = 0; b < eNoN; b++) { double T1 = wl*N(a)*N(b); - double T2 = (afm*ks + cs)*T1; - T1 = T1*ks; + int Bc = lFa.IEN(b,e); + double ks_b = robin_bc.get_stiffness(robin_bc.get_local_index(Bc)); + double cs_b = robin_bc.get_damping(robin_bc.get_local_index(Bc)); + double T2 = (afm*ks_b + cs_b)*T1; + T1 = T1*ks_b; // dM_1/dV_1 + af/am*dM_1/dU_1 lKd(0,a,b) = lKd(0,a,b) + T1*nDn(0,0); @@ -1612,7 +1641,8 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const do // cPhys != EquationType::phys_ustruct // } else { - double wl = w * (ks*afu + cs*afv); + // Use average stiffness and damping for this integration point + double wl = w * (ks_avg*afu + cs_avg*afv); for (int a = 0; a < eNoN; a++) { for (int b = 0; b < eNoN; b++) { @@ -1644,8 +1674,11 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const do for (int a = 0; a < eNoN; a++) { for (int b = 0; b < eNoN; b++) { double T1 = wl*N(a)*N(b); - double T2 = (afm*ks + cs)*T1; - T1 = T1*ks; + int Bc = lFa.IEN(b,e); + double ks_b = robin_bc.get_stiffness(robin_bc.get_local_index(Bc)); + double cs_b = robin_bc.get_damping(robin_bc.get_local_index(Bc)); + double T2 = (afm*ks_b + cs_b)*T1; + T1 = T1*ks_b; // dM_1/dV_1 + af/am*dM_1/dU_1 lKd(0,a,b) = lKd(0,a,b) + T1*nDn(0,0); @@ -1666,7 +1699,8 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const do } } else { - double wl = w * (ks*afu + cs*afv); + // Use average stiffness and damping for this integration point + double wl = w * (ks_avg*afu + cs_avg*afv); for (int a = 0; a < eNoN; a++) { for (int b = 0; b < eNoN; b++) { diff --git a/Code/Source/solver/set_bc.h b/Code/Source/solver/set_bc.h index d07d46413..85b11f530 100644 --- a/Code/Source/solver/set_bc.h +++ b/Code/Source/solver/set_bc.h @@ -33,6 +33,7 @@ #include "Simulation.h" #include "consts.h" +#include "RobinBoundaryCondition.h" #include @@ -61,7 +62,7 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array& Yg, const Array& Dg); void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, const Array& Yg, const Array& Dg); -void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const double ks, const double cs, const bool isN, +void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondition& robin_bc, const Array& Yg, const Array& Dg); void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa); diff --git a/tests/cases/struct/spatially_variable_robin/README.md b/tests/cases/struct/spatially_variable_robin/README.md new file mode 100644 index 000000000..3398e55b7 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/README.md @@ -0,0 +1,25 @@ +This test case simulates a spatially variable Robin boundary condition on a slab of material described by the Guccione material model. This case is identical to `ustruct/spatially_variable_robin`, except it uses `struct` physics. + +- Primary fibers run along the length of the slab (z-direction) and secondary fibers run across the width of the slab +(x-direction). + +- The slab is loaded on the +Y surface with a uniform pressure load. The load profile is a ramp to 10 dynes/cm^2 over 0.5 seconds, then held there until +2 seconds. The load is defined in `load.dat`, which can be generated with +`utilities/generate_boundary_condition_data/generate_load.py`. The load tends to push the slab downward. + +![Load Profile](load.png) + +- This is resisted by a spatially varying Robin boundary condition on the -Y surface. The stiffness is 0 at z = 0, and 50 at the far end. This is provided in `Y0_spatially_varying_robin.vtp`, which can be generated with `utilities/generate_boundary_condition_data/generate_spatially_variable_robin.py`. + +![Spatially varying Robin BC](Y0_spatially_varying_robin.png) + + +- The slab is also constrained by Dirichlet boundary conditions on the +-X and +-Z +surfaces, applied in the normal direction. These prevent the slab from moving +in the x and z directions. + +- The resulting deformation is shown in the video below: +![Deformation](animation.gif) +The black outline shows the initial configuration. As you can see, the displacement of the slab is greatest at z = 0, where the Robin BC stiffness is zero. At the far end of the slab, the displacement is very little, where the stiffness is greatest. The oscillations are due to the absence of any damping in the Robin BC. + +- Note, the deformation is fairly different from `ustruct/spatially_variable_robin`, like due to the different underlying physics formulation and the coarse mesh. Therefore, the `result_002.vtu` files are not identical. diff --git a/tests/cases/struct/spatially_variable_robin/Y0_spatially_varying_robin.png b/tests/cases/struct/spatially_variable_robin/Y0_spatially_varying_robin.png new file mode 100644 index 000000000..e802e1029 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/Y0_spatially_varying_robin.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8b8a4e6cc7a50194ea5d8385dcc96173133f8eb975f2d825ba20b234be18c519 +size 45762 diff --git a/tests/cases/struct/spatially_variable_robin/Y0_spatially_varying_robin.vtp b/tests/cases/struct/spatially_variable_robin/Y0_spatially_varying_robin.vtp new file mode 100644 index 000000000..dcc931105 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/Y0_spatially_varying_robin.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5b7fae7ab8ac1ba68df93cb69d6e427a9151148ab4e608bef4763fe09fd26642 +size 6131 diff --git a/tests/cases/struct/spatially_variable_robin/animation.avi b/tests/cases/struct/spatially_variable_robin/animation.avi new file mode 100644 index 000000000..8f196bf9e --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/animation.avi @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:52cdec41c9383bcfa36d65034a9b93686992f53dfb55f3f7aaaf28bdb54a52fc +size 6815946 diff --git a/tests/cases/struct/spatially_variable_robin/animation.gif b/tests/cases/struct/spatially_variable_robin/animation.gif new file mode 100644 index 000000000..2104af3e6 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/animation.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:80fe0515a87c25f3e1213733eddc72af2a27e939ee3b89ec17589623f73f3912 +size 3661841 diff --git a/tests/cases/struct/spatially_variable_robin/load.dat b/tests/cases/struct/spatially_variable_robin/load.dat new file mode 100644 index 000000000..cbf2bef2d --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/load.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1dbc63f619fec439c888c1a0e60c51061777e2a422043ad62155945ceb42f430 +size 1295 diff --git a/tests/cases/struct/spatially_variable_robin/load.png b/tests/cases/struct/spatially_variable_robin/load.png new file mode 100644 index 000000000..ef29d3980 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/load.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:a943274367734c7e77eaa120c300ffd89ca2ec4d743a5562e7360bf9dd122569 +size 17190 diff --git a/tests/cases/struct/spatially_variable_robin/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/spatially_variable_robin/mesh/mesh-complete.mesh.vtu new file mode 100644 index 000000000..f09242eb6 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d74418f571657908b5c4e73f7cc5ce040fd300ac444b1c893ab23001bb6e09ae +size 7499 diff --git a/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/X0.vtp b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/X0.vtp new file mode 100644 index 000000000..904401ad3 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f71480fecee6f177719ab68b14bb2b61bbbce8fc745c87576fd51f3fb0b3ef8 +size 3655 diff --git a/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/X1.vtp b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/X1.vtp new file mode 100644 index 000000000..20b2d0e58 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd46bcf90705c0e062aa658846aacb9091334a56fffac81715631628cd1f097f +size 3669 diff --git a/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Y0.vtp b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Y0.vtp new file mode 100644 index 000000000..19c7f0081 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ef44a9a441563692f58baccf576b3c10e6df4277046d1054db5ab821d498bf0 +size 4795 diff --git a/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Y1.vtp b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Y1.vtp new file mode 100644 index 000000000..7ac40cb0a --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4d2ad3c14f29583d9d2015b0c0ba681b66c0d2380e4b8851ca30a2cb81773183 +size 4629 diff --git a/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Z0.vtp b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Z0.vtp new file mode 100644 index 000000000..49112c41e --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6ee78716770c3ed486e5e81d32f90226e215bc095b5cda4f5c708218b84ee407 +size 3188 diff --git a/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Z1.vtp b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Z1.vtp new file mode 100644 index 000000000..22122b6fb --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/mesh/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c2bdc7813853e641c5a71c20f1cea215c9c0215523f51272c5650077614dcb29 +size 3245 diff --git a/tests/cases/struct/spatially_variable_robin/result_002.vtu b/tests/cases/struct/spatially_variable_robin/result_002.vtu new file mode 100644 index 000000000..3bc9c120a --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/result_002.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fac37b324afe68b2ad6039dce51d18713f9400d83da5ab1a8b610067d3771d75 +size 73861 diff --git a/tests/cases/struct/spatially_variable_robin/solver.xml b/tests/cases/struct/spatially_variable_robin/solver.xml new file mode 100644 index 000000000..546863d89 --- /dev/null +++ b/tests/cases/struct/spatially_variable_robin/solver.xml @@ -0,0 +1,149 @@ + + + + + 0 + 3 + 2 + 0.01 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/X0.vtp + + + + mesh/mesh-surfaces/X1.vtp + + + + mesh/mesh-surfaces/Y0.vtp + + + + mesh/mesh-surfaces/Y1.vtp + + + + mesh/mesh-surfaces/Z0.vtp + + + + mesh/mesh-surfaces/Z1.vtp + + + (0.0, 0.0, 1.0) + (1.0, 0.0, 0.0) + + + + + + + + + 1 + 4 + 1e-6 + + 1.0 + 1.0e6 + 0.48333 + + + 100.0 + 8.0 + 6.0 + 12.0 + + + ST91 + + + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-9 + 400 + + + + Dir + 0.0 + (0, 0, 1) + true + false + + + + Dir + 0.0 + (0, 0, 1) + true + false + + + + Dir + 0.0 + (1, 0, 0) + true + false + + + + Dir + 0.0 + (1, 0, 0) + true + false + + + + Robin + Y0_spatially_varying_robin.vtp + false + + + + Neu + Unsteady + load.dat + false + false + + + + + diff --git a/tests/cases/ustruct/spatially_variable_robin/README.md b/tests/cases/ustruct/spatially_variable_robin/README.md new file mode 100644 index 000000000..09aa6d5cc --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/README.md @@ -0,0 +1,26 @@ +This test case simulates a spatially variable Robin boundary condition on a slab of material described by the Guccione material model. This case is identical to `struct/spatially_variable_robin`, except it uses `ustruct` physics. + +- Primary fibers run along the length of the slab (z-direction) and secondary fibers run across the width of the slab +(x-direction). + +- The slab is loaded on the +Y surface with a uniform pressure load. The load profile is a ramp to 10 dynes/cm^2 over 0.5 seconds, then held there until +2 seconds. The load is defined in `load.dat`, which can be generated with +`utilities/generate_boundary_condition_data/generate_load.py`. The load tends to push the slab downward. + +![Load Profile](load.png) + +- This is resisted by a spatially varying Robin boundary condition on the -Y surface. The stiffness is 0 at z = 0, and 50 at the far end. This is provided in `Y0_spatially_varying_robin.vtp`, which can be generated with `utilities/generate_boundary_condition_data/generate_spatially_variable_robin.py`. + +![Spatially varying Robin BC](Y0_spatially_varying_robin.png) + + +- The slab is also constrained by Dirichlet boundary conditions on the +-X and +-Z +surfaces, applied in the normal direction. These prevent the slab from moving +in the x and z directions. + +- The resulting deformation is shown in the video below: +![Deformation](animation.gif) +The black outline shows the initial configuration. As you can see, the displacement of the slab is greatest at z = 0, where the Robin BC stiffness is zero. At the far end of the slab, the displacement is very little, where the stiffness is greatest. The oscillations are due to the absence of any damping in the Robin BC. + +- Note, the deformation is fairly different from `struct/spatially_variable_robin`, like due to the different underlying physics formulation and the coarse mesh. Therefore, the `result_002.vtu` files are not identical. + diff --git a/tests/cases/ustruct/spatially_variable_robin/Y0_spatially_varying_robin.png b/tests/cases/ustruct/spatially_variable_robin/Y0_spatially_varying_robin.png new file mode 100644 index 000000000..e802e1029 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/Y0_spatially_varying_robin.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8b8a4e6cc7a50194ea5d8385dcc96173133f8eb975f2d825ba20b234be18c519 +size 45762 diff --git a/tests/cases/ustruct/spatially_variable_robin/Y0_spatially_varying_robin.vtp b/tests/cases/ustruct/spatially_variable_robin/Y0_spatially_varying_robin.vtp new file mode 100644 index 000000000..dcc931105 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/Y0_spatially_varying_robin.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5b7fae7ab8ac1ba68df93cb69d6e427a9151148ab4e608bef4763fe09fd26642 +size 6131 diff --git a/tests/cases/ustruct/spatially_variable_robin/animation.avi b/tests/cases/ustruct/spatially_variable_robin/animation.avi new file mode 100644 index 000000000..762ff39b6 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/animation.avi @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1e25dd76f4e87c6e84d8ac4311709e8869351d2aba525f238627d9e56a217733 +size 6777356 diff --git a/tests/cases/ustruct/spatially_variable_robin/animation.gif b/tests/cases/ustruct/spatially_variable_robin/animation.gif new file mode 100644 index 000000000..f27e7e420 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/animation.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:60b561964c1381e8dafd1253e513618552ecfe8efc0348f787b84ca8a893ebbc +size 2662980 diff --git a/tests/cases/ustruct/spatially_variable_robin/load.dat b/tests/cases/ustruct/spatially_variable_robin/load.dat new file mode 100644 index 000000000..cbf2bef2d --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/load.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1dbc63f619fec439c888c1a0e60c51061777e2a422043ad62155945ceb42f430 +size 1295 diff --git a/tests/cases/ustruct/spatially_variable_robin/load.png b/tests/cases/ustruct/spatially_variable_robin/load.png new file mode 100644 index 000000000..ef29d3980 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/load.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:a943274367734c7e77eaa120c300ffd89ca2ec4d743a5562e7360bf9dd122569 +size 17190 diff --git a/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-complete.mesh.vtu b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-complete.mesh.vtu new file mode 100644 index 000000000..f09242eb6 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d74418f571657908b5c4e73f7cc5ce040fd300ac444b1c893ab23001bb6e09ae +size 7499 diff --git a/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/X0.vtp b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/X0.vtp new file mode 100644 index 000000000..904401ad3 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f71480fecee6f177719ab68b14bb2b61bbbce8fc745c87576fd51f3fb0b3ef8 +size 3655 diff --git a/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/X1.vtp b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/X1.vtp new file mode 100644 index 000000000..20b2d0e58 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd46bcf90705c0e062aa658846aacb9091334a56fffac81715631628cd1f097f +size 3669 diff --git a/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Y0.vtp b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Y0.vtp new file mode 100644 index 000000000..19c7f0081 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ef44a9a441563692f58baccf576b3c10e6df4277046d1054db5ab821d498bf0 +size 4795 diff --git a/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Y1.vtp b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Y1.vtp new file mode 100644 index 000000000..7ac40cb0a --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4d2ad3c14f29583d9d2015b0c0ba681b66c0d2380e4b8851ca30a2cb81773183 +size 4629 diff --git a/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Z0.vtp b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Z0.vtp new file mode 100644 index 000000000..49112c41e --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6ee78716770c3ed486e5e81d32f90226e215bc095b5cda4f5c708218b84ee407 +size 3188 diff --git a/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Z1.vtp b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Z1.vtp new file mode 100644 index 000000000..22122b6fb --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/mesh/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c2bdc7813853e641c5a71c20f1cea215c9c0215523f51272c5650077614dcb29 +size 3245 diff --git a/tests/cases/ustruct/spatially_variable_robin/result_002.vtu b/tests/cases/ustruct/spatially_variable_robin/result_002.vtu new file mode 100644 index 000000000..03a54e2da --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/result_002.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:944fbebcb555aed331ae0f659370027cda7576041961bd6fb39bb105e44680ea +size 78251 diff --git a/tests/cases/ustruct/spatially_variable_robin/solver.xml b/tests/cases/ustruct/spatially_variable_robin/solver.xml new file mode 100644 index 000000000..2a1a7bee5 --- /dev/null +++ b/tests/cases/ustruct/spatially_variable_robin/solver.xml @@ -0,0 +1,150 @@ + + + + + 0 + 3 + 2 + 0.01 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/X0.vtp + + + + mesh/mesh-surfaces/X1.vtp + + + + mesh/mesh-surfaces/Y0.vtp + + + + mesh/mesh-surfaces/Y1.vtp + + + + mesh/mesh-surfaces/Z0.vtp + + + + mesh/mesh-surfaces/Z1.vtp + + + (0.0, 0.0, 1.0) + (1.0, 0.0, 0.0) + + + + + + + + + 1 + 4 + 1e-6 + + 1.0 + 1.0e6 + 0.48333 + + + 100.0 + 8.0 + 6.0 + 12.0 + + + ST91 + + + true + true + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-9 + 400 + + + + Dir + 0.0 + (0, 0, 1) + true + false + + + + Dir + 0.0 + (0, 0, 1) + true + false + + + + Dir + 0.0 + (1, 0, 0) + true + false + + + + Dir + 0.0 + (1, 0, 0) + true + false + + + + Robin + Y0_spatially_varying_robin.vtp + + + + Neu + Unsteady + load.dat + false + false + + + + + diff --git a/tests/test_struct.py b/tests/test_struct.py index 6ba9db89f..35f1824a3 100644 --- a/tests/test_struct.py +++ b/tests/test_struct.py @@ -55,6 +55,10 @@ def test_robin(n_proc): test_folder = "robin" run_with_reference(base_folder, test_folder, fields, n_proc) +def test_spatially_variable_robin(n_proc): + test_folder = "spatially_variable_robin" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=2) + def test_LV_NeoHookean_passive(n_proc): test_folder = "LV_NeoHookean_passive" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=5) diff --git a/tests/test_ustruct.py b/tests/test_ustruct.py index e5d96993c..014a46d27 100644 --- a/tests/test_ustruct.py +++ b/tests/test_ustruct.py @@ -70,6 +70,10 @@ def test_LV_NeoHookean_passive_sv0D(n_proc): run_with_reference(base_folder, test_folder, fields, n_proc, t_max=3) +def test_spatially_variable_robin(n_proc): + test_folder = "spatially_variable_robin" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=2) + def test_tensile_adventitia_Newtonian_viscosity(n_proc): test_folder = "tensile_adventitia_Newtonian_viscosity" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) diff --git a/utilities/generate_boundary_condition_data/generate_load.py b/utilities/generate_boundary_condition_data/generate_load.py new file mode 100644 index 000000000..d08344dbc --- /dev/null +++ b/utilities/generate_boundary_condition_data/generate_load.py @@ -0,0 +1,32 @@ +import numpy as np +import os + +# Go to directory of this script +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Number of timesteps and number of Fourier modes +n_timesteps = 101 +n_modes = 64 + +# Generate time values from 0 to 2 +time = np.linspace(0, 2, n_timesteps) + +# Generate ramp from 0 to 10 in 0.5 seconds, then hold until 2 seconds. +load = np.zeros(n_timesteps) +load[time < 0.5] = 10 * time[time < 0.5] / 0.5 +load[(time >= 0.5)] = 10 + + +# Write the time and stress values to a text file +with open("load.dat", "w") as file: + file.write(f"{n_timesteps} {n_modes}\n") + for t, s in zip(time, load): + file.write(f"{t:.3f} {s:.3f}\n") + +# Plot the stress values +import matplotlib.pyplot as plt +plt.plot(time, load) +plt.xlabel("Time (s)") +plt.ylabel("Load (dynes/cm^2)") +plt.title("Load") +plt.savefig("load.png") \ No newline at end of file diff --git a/utilities/generate_boundary_condition_data/generate_spatially_variable_robin.py b/utilities/generate_boundary_condition_data/generate_spatially_variable_robin.py new file mode 100644 index 000000000..10bb76a06 --- /dev/null +++ b/utilities/generate_boundary_condition_data/generate_spatially_variable_robin.py @@ -0,0 +1,295 @@ +#!/usr/bin/env python3 +""" +Script to generate spatially varying Robin boundary condition VTP files. + +This script reads a VTP file and creates a new VTP file with spatially varying +stiffness and damping values based on user-defined functions of the coordinates (x, y, z). + +To use this script: +1. Edit the configuration section below to set your input/output files and functions +2. Run: python generate_spatially_variable_robin.py + +The script supports mathematical expressions using x, y, z coordinates and common +mathematical functions (sin, cos, exp, sqrt, etc.). + +Requirements: + - numpy + - pyvista +""" + +import numpy as np +import pyvista as pv +import os +import math + + +def safe_eval(expression: str, x: float, y: float, z: float) -> float: + """ + Safely evaluate a mathematical expression with x, y, z variables. + + Args: + expression: Mathematical expression as string + x, y, z: Coordinate values + + Returns: + Evaluated result + + Raises: + ValueError: If expression is invalid or contains unsafe operations + """ + # Define safe namespace for evaluation + safe_dict = { + 'x': x, 'y': y, 'z': z, + 'sin': math.sin, 'cos': math.cos, 'tan': math.tan, + 'exp': math.exp, 'log': math.log, 'log10': math.log10, + 'sqrt': math.sqrt, 'abs': abs, 'pow': pow, + 'pi': math.pi, 'e': math.e, + 'min': min, 'max': max, + '__builtins__': {} + } + + try: + return float(eval(expression, safe_dict)) + except Exception as e: + raise ValueError(f"Error evaluating expression '{expression}': {e}") + + +def read_vtp_file(filepath: str) -> pv.PolyData: + """ + Read a VTP file and return the PyVista polydata object. + + Args: + filepath: Path to the VTP file + + Returns: + PyVista polydata object + + Raises: + FileNotFoundError: If file doesn't exist + RuntimeError: If file cannot be read + """ + if not os.path.exists(filepath): + raise FileNotFoundError(f"VTP file not found: {filepath}") + + try: + return pv.read(filepath) + except Exception as e: + raise RuntimeError(f"Failed to read VTP file '{filepath}': {e}") + + +def write_vtp_file(polydata: pv.PolyData, filepath: str) -> None: + """ + Write a polydata object to a VTP file. + + Args: + polydata: PyVista polydata object + filepath: Output file path + + Raises: + RuntimeError: If file cannot be written + """ + try: + polydata.save(filepath, binary=True) + except Exception as e: + raise RuntimeError(f"Failed to write VTP file '{filepath}': {e}") + + +def get_coordinates(polydata: pv.PolyData) -> np.ndarray: + """ + Extract coordinates from polydata. + + Args: + polydata: PyVista polydata object + + Returns: + Array of shape (n_points, 3) with x, y, z coordinates + """ + if polydata.n_points == 0: + raise ValueError("No points found in VTP file") + + return polydata.points + + +def add_point_data(polydata: pv.PolyData, name: str, data: np.ndarray) -> None: + """ + Add point data array to polydata. + + Args: + polydata: PyVista polydata object + name: Name of the data array + data: 1D numpy array of data values + """ + polydata[name] = data + + +def generate_spatially_varying_robin_bc( + input_vtp: str, + output_vtp: str, + stiffness_func: str, + damping_func: str, + stiffness_scale: float = 1.0, + damping_scale: float = 1.0, + min_stiffness: float = 0.0, + min_damping: float = 0.0, + verbose: bool = False +) -> None: + """ + Generate spatially varying Robin boundary condition VTP file. + + Args: + input_vtp: Path to input VTP file + output_vtp: Path to output VTP file + stiffness_func: Mathematical expression for stiffness as function of x, y, z + damping_func: Mathematical expression for damping as function of x, y, z + stiffness_scale: Scaling factor for stiffness values + damping_scale: Scaling factor for damping values + min_stiffness: Minimum allowed stiffness value + min_damping: Minimum allowed damping value + verbose: Print detailed information + """ + if verbose: + print(f"Reading VTP file: {input_vtp}") + + # Read input VTP file + polydata = read_vtp_file(input_vtp) + coords = get_coordinates(polydata) + n_points = polydata.n_points + + if verbose: + print(f"Found {n_points} points in the mesh") + print(f"Coordinate ranges:") + print(f" X: [{coords[:, 0].min():.6f}, {coords[:, 0].max():.6f}]") + print(f" Y: [{coords[:, 1].min():.6f}, {coords[:, 1].max():.6f}]") + print(f" Z: [{coords[:, 2].min():.6f}, {coords[:, 2].max():.6f}]") + + # Initialize arrays for stiffness and damping + stiffness_values = np.zeros(n_points) + damping_values = np.zeros(n_points) + + # Evaluate functions at each point + if verbose: + print("Evaluating stiffness and damping functions...") + + for i in range(n_points): + x, y, z = coords[i, :] + + try: + # Evaluate stiffness function + stiffness_raw = safe_eval(stiffness_func, x, y, z) + stiffness_values[i] = max(min_stiffness, stiffness_scale * stiffness_raw) + + # Evaluate damping function + damping_raw = safe_eval(damping_func, x, y, z) + damping_values[i] = max(min_damping, damping_scale * damping_raw) + + except ValueError as e: + raise ValueError(f"Error at point {i} (x={x:.6f}, y={y:.6f}, z={z:.6f}): {e}") + + if verbose: + print(f"Stiffness range: [{stiffness_values.min():.6e}, {stiffness_values.max():.6e}]") + print(f"Damping range: [{damping_values.min():.6e}, {damping_values.max():.6e}]") + + # Create output polydata (copy of input) + output_polydata = polydata.copy() + + # Add stiffness and damping arrays + add_point_data(output_polydata, "Stiffness", stiffness_values) + add_point_data(output_polydata, "Damping", damping_values) + + # Show mesh with Stiffness + output_polydata.plot(scalars="Stiffness", show_edges=True, cmap="viridis") + + # Write output file + if verbose: + print(f"Writing output VTP file: {output_vtp}") + + write_vtp_file(output_polydata, output_vtp) + + if verbose: + print("Successfully generated spatially varying Robin BC VTP file!") + + +# ============================================================================= +# CONFIGURATION SECTION - EDIT THESE VALUES +# ============================================================================= + +# Input and output file paths +os.chdir(os.path.dirname(os.path.realpath(__file__))) +INPUT_VTP_FILE = "mesh/mesh-surfaces/Y0.vtp" # Path to input VTP file +OUTPUT_VTP_FILE = "Y0_spatially_varying_robin.vtp" # Path to output VTP file + +# Mathematical expressions for stiffness and damping as functions of x, y, z +# Available functions: sin, cos, tan, exp, log, log10, sqrt, abs, pow, min, max +# Available constants: pi, e +STIFFNESS_FUNCTION = "100 * z" # Example: linear variation in z-direction +DAMPING_FUNCTION = "0" # Example: no damping + +# Scaling factors +STIFFNESS_SCALE = 1.0 # Scaling factor for stiffness values +DAMPING_SCALE = 1.0 # Scaling factor for damping values + +# Minimum values (stiffness and damping cannot be negative) +MIN_STIFFNESS = 0.0 +MIN_DAMPING = 0.0 + +# Verbose output +VERBOSE = True + +# ============================================================================= +# EXAMPLE FUNCTIONS - UNCOMMENT AND MODIFY AS NEEDED +# ============================================================================= + +# Linear variation in x-direction +# STIFFNESS_FUNCTION = "1000 * x" +# DAMPING_FUNCTION = "10 * x" + +# Quadratic variation from center +# STIFFNESS_FUNCTION = "1000 * (x**2 + y**2 + z**2)" +# DAMPING_FUNCTION = "10" + +# Sinusoidal variation +# STIFFNESS_FUNCTION = "1000 * (1 + sin(2*pi*x))" +# DAMPING_FUNCTION = "10 * cos(pi*y)" + +# Exponential decay +# STIFFNESS_FUNCTION = "1000 * exp(-x)" +# DAMPING_FUNCTION = "10 * exp(-0.5*x)" + +# Step function +# STIFFNESS_FUNCTION = "1000 * (x > 2.5)" +# DAMPING_FUNCTION = "10 * (y > 2.5)" + +# Gaussian-like distribution +# STIFFNESS_FUNCTION = "1000 * exp(-(x**2 + y**2 + z**2))" +# DAMPING_FUNCTION = "10" + +# ============================================================================= +# MAIN EXECUTION +# ============================================================================= + +if __name__ == "__main__": + try: + print("Generating spatially varying Robin boundary condition VTP file...") + print(f"Input file: {INPUT_VTP_FILE}") + print(f"Output file: {OUTPUT_VTP_FILE}") + print(f"Stiffness function: {STIFFNESS_FUNCTION}") + print(f"Damping function: {DAMPING_FUNCTION}") + print() + + generate_spatially_varying_robin_bc( + input_vtp=INPUT_VTP_FILE, + output_vtp=OUTPUT_VTP_FILE, + stiffness_func=STIFFNESS_FUNCTION, + damping_func=DAMPING_FUNCTION, + stiffness_scale=STIFFNESS_SCALE, + damping_scale=DAMPING_SCALE, + min_stiffness=MIN_STIFFNESS, + min_damping=MIN_DAMPING, + verbose=VERBOSE + ) + + print("✓ Successfully generated spatially varying Robin BC VTP file!") + + except Exception as e: + print(f"✗ Error: {e}") + exit(1)