diff --git a/Code/Source/solver/BoundaryCondition.cpp b/Code/Source/solver/BoundaryCondition.cpp index 736f40382..ffd0bd4b5 100644 --- a/Code/Source/solver/BoundaryCondition.cpp +++ b/Code/Source/solver/BoundaryCondition.cpp @@ -40,7 +40,7 @@ #define n_debug_bc -BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std::vector& array_names, const StringBoolMap& flags, const faceType& face) +BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std::vector& array_names, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger) : face_(&face) , global_num_nodes_(face.nNo) , local_num_nodes_(0) @@ -48,6 +48,7 @@ BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std , spatially_variable(true) , vtp_file_path_(vtp_file_path) , flags_(flags) + , logger_(&logger) { try { global_data_ = read_data_from_vtp_file(vtp_file_path, array_names); @@ -73,13 +74,14 @@ BoundaryCondition::BoundaryCondition(const std::string& vtp_file_path, const std } } -BoundaryCondition::BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face) +BoundaryCondition::BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger) : face_(&face) , global_num_nodes_(face.nNo) , local_num_nodes_(0) , spatially_variable(false) , vtp_file_path_("") , flags_(flags) + , logger_(&logger) { try { // Store array names, validate and store values @@ -107,6 +109,7 @@ BoundaryCondition::BoundaryCondition(const BoundaryCondition& other) , vtp_file_path_(other.vtp_file_path_) , flags_(other.flags_) , global_node_map_(other.global_node_map_) + , logger_(other.logger_) { if (other.vtp_data_) { vtp_data_ = std::make_unique(*other.vtp_data_); @@ -126,6 +129,7 @@ void swap(BoundaryCondition& lhs, BoundaryCondition& rhs) noexcept { swap(lhs.flags_, rhs.flags_); swap(lhs.global_node_map_, rhs.global_node_map_); swap(lhs.vtp_data_, rhs.vtp_data_); + swap(lhs.logger_, rhs.logger_); } BoundaryCondition& BoundaryCondition::operator=(BoundaryCondition other) { @@ -145,6 +149,7 @@ BoundaryCondition::BoundaryCondition(BoundaryCondition&& other) noexcept , flags_(std::move(other.flags_)) , global_node_map_(std::move(other.global_node_map_)) , vtp_data_(std::move(other.vtp_data_)) + , logger_(other.logger_) { other.face_ = nullptr; other.global_num_nodes_ = 0; @@ -203,6 +208,11 @@ BoundaryCondition::StringArrayMap BoundaryCondition::read_data_from_vtp_file(con #endif } + logger_ -> log_message("[BoundaryCondition] Loaded from VTP file"); + logger_ -> log_message("\t File path:", vtp_file_path); + logger_ -> log_message("\t Arrays:", array_names); + logger_ -> log_message("\t Face:", face_->name); + return result; } @@ -375,11 +385,20 @@ void BoundaryCondition::distribute_spatially_variable(const ComMod& com_mod, con // Get VTP points for position matching Array vtp_points = vtp_data_->get_points(); - // Look up data for all nodes using point matching + // Get mesh scale factor from the face's mesh + double mesh_scale_factor = 1.0; // Default scale factor + if (face_ != nullptr) { + mesh_scale_factor = com_mod.msh[face_->iM].scF; + #ifdef debug_distribute_spatially_variable + dmsg << "Mesh scale factor: " << mesh_scale_factor << std::endl; + #endif + } + + // 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); + int vtp_idx = find_vtp_point_index(all_positions(0,i), all_positions(1,i), all_positions(2,i), vtp_points, mesh_scale_factor); all_values[array_name](i) = global_data_[array_name](vtp_idx, 0); } } @@ -474,10 +493,13 @@ void BoundaryCondition::distribute_flags(const CmMod& cm_mod, const cmType& cm, } int BoundaryCondition::find_vtp_point_index(double x, double y, double z, - const Array& vtp_points) const + const Array& vtp_points, double mesh_scale_factor) const { const int num_points = vtp_points.ncols(); - Vector target_point{x, y, z}; + + // Scale down the target coordinates to match the unscaled VTP coordinates + // The simulation coordinates are scaled by mesh_scale_factor, but VTP coordinates are not + Vector target_point{x / mesh_scale_factor, y / mesh_scale_factor, z / mesh_scale_factor}; // Simple linear search through all points in the VTP file for (int i = 0; i < num_points; i++) { @@ -490,6 +512,7 @@ int BoundaryCondition::find_vtp_point_index(double x, double y, double z, #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 << "Scaled target position (" << target_point(0) << ", " << target_point(1) << ", " << target_point(2) << ")" << std::endl; dmsg << "VTP point index: " << i << std::endl; #endif diff --git a/Code/Source/solver/BoundaryCondition.h b/Code/Source/solver/BoundaryCondition.h index 7105c8e0e..536496e43 100644 --- a/Code/Source/solver/BoundaryCondition.h +++ b/Code/Source/solver/BoundaryCondition.h @@ -36,6 +36,7 @@ #include "CmMod.h" #include "VtkData.h" #include +#include "SimulationLogger.h" #include #include #include @@ -79,7 +80,7 @@ class BoundaryCondition { using StringDoubleMap = std::map; /// @brief Data members for BC - const faceType* face_ = nullptr; ///< Face associated with the BC (can be null) + const faceType* face_ = nullptr; ///< Face associated with the BC (not owned by BoundaryCondition) 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 @@ -90,6 +91,7 @@ class BoundaryCondition { 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 + const SimulationLogger* logger_ = nullptr; ///< Logger for warnings/info (not owned by BoundaryCondition) public: /// @brief Tolerance for point matching in VTP files @@ -102,13 +104,15 @@ class BoundaryCondition { /// @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 + /// @param logger Simulation logger used to write warnings /// @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); + BoundaryCondition(const std::string& vtp_file_path, const std::vector& array_names, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger); /// @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); + /// @param logger Simulation logger used to write warnings + BoundaryCondition(const StringDoubleMap& uniform_values, const StringBoolMap& flags, const faceType& face, SimulationLogger& logger); /// @brief Copy constructor BoundaryCondition(const BoundaryCondition& other); @@ -198,9 +202,10 @@ class BoundaryCondition { /// @param y Y coordinate /// @param z Z coordinate /// @param vtp_points VTP points array + /// @param mesh_scale_factor Scale factor applied to mesh coordinates /// @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; + int find_vtp_point_index(double x, double y, double z, const Array& vtp_points, double mesh_scale_factor) const; /// @brief Hook for derived classes to validate array values /// @param array_name Name of the array being validated diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 02cffa0a4..a818ab9c6 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -443,7 +443,6 @@ 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 e9fd424d1..b70b72246 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -806,7 +806,6 @@ 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 index dd16a9f2b..fabb8e63f 100644 --- a/Code/Source/solver/RobinBoundaryCondition.cpp +++ b/Code/Source/solver/RobinBoundaryCondition.cpp @@ -29,6 +29,17 @@ */ #include "RobinBoundaryCondition.h" +#include #define n_debug_robin_bc +RobinBoundaryCondition::RobinBoundaryCondition(double uniform_stiffness, double uniform_damping, bool normal_only, const faceType& face, SimulationLogger& logger) + : BoundaryCondition({{"Stiffness", uniform_stiffness}, {"Damping", uniform_damping}}, StringBoolMap{{"normal_direction_only", normal_only}}, face, logger) +{ + // Warning if both stiffness and damping are zero + if (uniform_stiffness == 0.0 && uniform_damping == 0.0) { + logger_ -> log_message("WARNING [RobinBoundaryCondition] Both stiffness and damping values set to zero. " + "This will result in effectively no boundary condition being applied."); + } +} + diff --git a/Code/Source/solver/RobinBoundaryCondition.h b/Code/Source/solver/RobinBoundaryCondition.h index 13b53241c..702adbd26 100644 --- a/Code/Source/solver/RobinBoundaryCondition.h +++ b/Code/Source/solver/RobinBoundaryCondition.h @@ -70,11 +70,12 @@ class RobinBoundaryCondition : public BoundaryCondition { /// @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 + /// @param logger Simulation logger used to write warnings /// @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) {} + RobinBoundaryCondition(const std::string& vtp_file_path, bool normal_only, const faceType& face, SimulationLogger& logger) + : BoundaryCondition(vtp_file_path, std::vector{"Stiffness", "Damping"}, StringBoolMap{{"normal_direction_only", normal_only}}, face, logger) {} /// @brief Constructor for uniform values @@ -82,9 +83,9 @@ class RobinBoundaryCondition : public BoundaryCondition { /// @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 + /// @param logger Simulation logger used to write warnings /// @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) {} + RobinBoundaryCondition(double uniform_stiffness, double uniform_damping, bool normal_only, const faceType& face, SimulationLogger& logger); /// @brief Apply only along normal direction (getter) /// @return true if BC should be applied only along normal direction diff --git a/Code/Source/solver/SimulationLogger.h b/Code/Source/solver/SimulationLogger.h index 26369e87a..536d62fec 100755 --- a/Code/Source/solver/SimulationLogger.h +++ b/Code/Source/solver/SimulationLogger.h @@ -47,7 +47,7 @@ class SimulationLogger { this->initialize(file_name, cout_write); } - void initialize(const std::string& file_name, bool cout_write=false) + void initialize(const std::string& file_name, bool cout_write=false) const { log_file_.open(file_name); if (log_file_.fail()) { @@ -58,12 +58,14 @@ class SimulationLogger { file_name_ = file_name; } + bool is_initialized() const { return log_file_.is_open(); } + ~SimulationLogger() { log_file_.close(); } - template SimulationLogger& operator<< (const T& value) + template const SimulationLogger& operator<< (const T& value) const { if (file_name_ == "") { return *this; @@ -78,7 +80,28 @@ class SimulationLogger { return *this; } - SimulationLogger& operator<<(std::ostream&(*value)(std::ostream& o)) + // Special handling for vector + const SimulationLogger& operator<< (const std::vector& values) const + { + if (file_name_ == "") { + return *this; + } + + bool first = true; + for (const auto& value : values) { + if (!first) { + log_file_ << ", "; + if (cout_write_) std::cout << ", "; + } + log_file_ << value; + if (cout_write_) std::cout << value; + first = false; + } + + return *this; + } + + const SimulationLogger& operator<<(std::ostream&(*value)(std::ostream& o)) const { if (file_name_ == "") { return *this; @@ -93,10 +116,27 @@ class SimulationLogger { return *this; }; + /// @brief Log a message with automatic space separation + /// @param args Arguments to log + template + const SimulationLogger& log_message(const Args&... args) const { + if (file_name_ == "") return *this; + + bool first = true; + ((first ? (first = false, (*this << args)) : (*this << ' ' << args)), ...); + *this << std::endl; + return *this; + } + private: - bool cout_write_ = false; - std::string file_name_; - std::ofstream log_file_; + // These members are marked mutable because they can be modified in const member functions. + // While logging writes to a file (modifying these members), it doesn't change the logical + // state of the logger - this is exactly what mutable is designed for: implementation + // details that need to change even when the object's public state remains constant. + mutable bool cout_write_ = false; + mutable bool initialized_ = false; + mutable std::string file_name_; + mutable std::ofstream log_file_; }; #endif diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index e871c4bda..0e29f5dec 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -371,22 +371,6 @@ void initialize(Simulation* simulation, Vector& timeP) dmsg.banner(); #endif - // Setup logging to history file. - if (!simulation->com_mod.cm.slv(simulation->cm_mod)) { - std::string hist_file_name; - - if (chnl_mod.appPath != "") { - auto mkdir_arg = std::string("mkdir -p ") + chnl_mod.appPath; - std::system(mkdir_arg.c_str()); - hist_file_name = chnl_mod.appPath + "/" + simulation->history_file_name; - } else { - hist_file_name = simulation->history_file_name; - } - - bool output_to_cout = true; - simulation->logger.initialize(hist_file_name, output_to_cout); - } - #ifdef debug_initialize dmsg << "Set time " << std::endl; dmsg << "com_mod.timeP[0]: " << com_mod.timeP[0] << std::endl; diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 8ec89f3b2..f7c2ccb14 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -394,15 +394,17 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Robin))) { // 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(), + if (bc_params->spatial_values_file_path.defined()) { + lBc.robin_bc = RobinBoundaryCondition(bc_params->spatial_values_file_path.value(), bc_params->apply_along_normal_direction.value(), - com_mod.msh[lBc.iM].fa[lBc.iFa]); + com_mod.msh[lBc.iM].fa[lBc.iFa], + simulation->logger); } 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]); + com_mod.msh[lBc.iM].fa[lBc.iFa], + simulation->logger); } } @@ -1703,7 +1705,23 @@ void read_files(Simulation* simulation, const std::string& file_name) com_mod.stFileFlag = gen_params.continue_previous_simulation.value(); } - // Set simulatioin and module member data from XML parameters. + // Setup logging to history file. + if (!simulation->com_mod.cm.slv(simulation->cm_mod)) { + std::string hist_file_name; + + if (chnl_mod.appPath != "") { + auto mkdir_arg = std::string("mkdir -p ") + chnl_mod.appPath; + std::system(mkdir_arg.c_str()); + hist_file_name = chnl_mod.appPath + "/" + simulation->history_file_name; + } else { + hist_file_name = simulation->history_file_name; + } + + bool output_to_cout = true; + simulation->logger.initialize(hist_file_name, output_to_cout); + } + + // Set simulation and module member data from XML parameters. simulation->set_module_parameters(); // Read mesh and BCs data. diff --git a/tests/cases/struct/spatially_variable_robin/solver.xml b/tests/cases/struct/spatially_variable_robin/solver.xml index 546863d89..6b26a3309 100644 --- a/tests/cases/struct/spatially_variable_robin/solver.xml +++ b/tests/cases/struct/spatially_variable_robin/solver.xml @@ -132,7 +132,7 @@ Robin - Y0_spatially_varying_robin.vtp + Y0_spatially_varying_robin.vtp false diff --git a/tests/cases/ustruct/spatially_variable_robin/solver.xml b/tests/cases/ustruct/spatially_variable_robin/solver.xml index 2a1a7bee5..57ac72dea 100644 --- a/tests/cases/ustruct/spatially_variable_robin/solver.xml +++ b/tests/cases/ustruct/spatially_variable_robin/solver.xml @@ -134,7 +134,7 @@ Robin - Y0_spatially_varying_robin.vtp + Y0_spatially_varying_robin.vtp