From adeb2b4a3aa44c8735ba6ea62faa539897821909 Mon Sep 17 00:00:00 2001 From: Han Zhao Date: Mon, 10 Feb 2025 17:41:00 -0800 Subject: [PATCH 01/19] Fixed spatial dimension bug in bw_fluid_3d and index error in set_bc_dir_wl --- .../pipe_RCR_weak_dir_3d/lumen_inlet.fcs | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100755 tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs diff --git a/tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs b/tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs new file mode 100755 index 000000000..844483140 --- /dev/null +++ b/tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs @@ -0,0 +1,21 @@ + 0.0000 + 1.0000 + 0.000000E+00 + 0.000000E+00 + 16 + -6.283185E+01 0.000000E+00 + 6.263025E+01 -1.868984E-07 + -1.851028E-07 1.979808E-15 + -3.293147E-07 1.821408E-07 + 2.967629E-07 4.499564E-17 + -7.134611E-07 -1.729158E-07 + -1.667927E-07 5.199497E-16 + -1.920733E-06 1.597801E-07 + 5.066059E-08 -2.587250E-16 + 9.364629E-07 -1.435117E-07 + -1.344913E-07 8.999129E-17 + -9.410540E-07 1.250477E-07 + 1.921845E-07 2.874722E-17 + 9.540060E-07 -1.054107E-07 + -9.547580E-08 -1.106526E-16 + 1.016783E-06 8.563008E-08 From 306c9b4202a93c7f2555ed6a5e9fd97f4e16b846 Mon Sep 17 00:00:00 2001 From: Han Zhao Date: Tue, 25 Feb 2025 14:19:19 -0800 Subject: [PATCH 02/19] Adding explicit RIS code. Resolving conflits --- Code/Source/solver/CmMod.h | 1 + Code/Source/solver/ComMod.h | 3 +++ Code/Source/solver/Parameters.cpp | 42 +++++++++++++++++++++++++++++++ Code/Source/solver/Parameters.h | 26 +++++++++++++++++++ Code/Source/solver/initialize.cpp | 18 +++++++++++++ 5 files changed, 90 insertions(+) diff --git a/Code/Source/solver/CmMod.h b/Code/Source/solver/CmMod.h index 36ac60e18..4e355d83d 100644 --- a/Code/Source/solver/CmMod.h +++ b/Code/Source/solver/CmMod.h @@ -79,6 +79,7 @@ class cmType { void bcast(const CmMod& cm_mod, double* data) const; void bcast(const CmMod& cm_mod, Vector& data, const std::string& name="") const; void bcast(const CmMod& cm_mod, Array& data, const std::string& name="") const; + void bcast(const CmMod& cm_mod, Array& data, const std::string& name="") const; void bcast(const CmMod& cm_mod, int* data) const; void bcast(const CmMod& cm_mod, Vector& data) const; diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 5a9de882d..a688145b6 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -135,6 +135,9 @@ class bcType // Strong/Weak application of Dirichlet BC int clsFlgRis = 0; + // Strong/Weak application of Dirichlet BC + int clsFlgRis = 0; + // Pre/Res/Flat/Para... boundary types // // This stores differnt BCs as bitwise values. diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index a70042a72..cd8f20694 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2606,6 +2606,48 @@ void PrecomputedSolutionParameters::set_values(tinyxml2::XMLElement* xml_elem) xml_util_set_parameters(ftpr, xml_elem, error_msg); } +////////////////////////////////////////////////////////// +// RISProjectionParameters // +////////////////////////////////////////////////////////// + +/// @brief Define the XML element name for mesh parameters. +const std::string RISProjectionParameters::xml_element_name_ = "Add_RIS_projection"; + +RISProjectionParameters::RISProjectionParameters() +{ + // A parameter that must be defined. + bool required = true; + + name = Parameter("name", "", required); + + set_parameter("Project_from_face", "", required, project_from_face); + set_parameter("Resistance", 1.e6, !required, resistance); + set_parameter("Projection_tolerance", 0.0, !required, projection_tolerance); +} + +void RISProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + + // Get the 'type' from the element. + const char* sname; + auto result = xml_elem->QueryStringAttribute("name", &sname); + if (sname == nullptr) { + throw std::runtime_error("No TYPE given in the XML element."); + } + name.set(std::string(sname)); + + using std::placeholders::_1; + using std::placeholders::_2; + + std::function ftpr = + std::bind( &RISProjectionParameters::set_parameter_value, *this, _1, _2); + + xml_util_set_parameters(ftpr, xml_elem, error_msg); +} + + ////////////////////////////////////////////////////////// // P r o j e c t i o n P a r a m e t e r s // ////////////////////////////////////////////////////////// diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 23bb44c41..9aa79bae1 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -864,6 +864,30 @@ class PrecomputedSolutionParameters: public ParameterLists Parameter use_precomputed_solution; }; +/// @brief The RISProjectionParameters class stores parameters for the +/// 'Add_RIS_projection' XML element used for RIS valve simulations. +/// \code {.xml} +/// +/// right_ris +/// 1.e6 +/// +/// \endcode +class RISProjectionParameters : public ParameterLists +{ + public: + RISProjectionParameters(); + + void set_values(tinyxml2::XMLElement* xml_elem); + + static const std::string xml_element_name_; + + Parameter name; + + Parameter project_from_face; + Parameter resistance; + Parameter projection_tolerance; +}; + /// @brief The ProjectionParameters class stores parameters for the /// 'Add_projection' XML element used for fluid-structure interaction /// simulations. @@ -1636,6 +1660,7 @@ class Parameters { void set_equation_values(tinyxml2::XMLElement* root_element); void set_mesh_values(tinyxml2::XMLElement* root_element); void set_precomputed_solution_values(tinyxml2::XMLElement* root_element); + void set_RIS_projection_values(tinyxml2::XMLElement* root_element); void set_projection_values(tinyxml2::XMLElement* root_element); void set_svzerodsolver_interface_values(tinyxml2::XMLElement* root_element); @@ -1647,6 +1672,7 @@ class Parameters { GeneralSimulationParameters general_simulation_parameters; std::vector mesh_parameters; std::vector equation_parameters; + std::vector RIS_projection_parameters; std::vector projection_parameters; PrecomputedSolutionParameters precomputed_solution_parameters; diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7f1e1ae84..b99b7efd4 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -123,6 +123,7 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); } else if (risFlag) { bin_file.read((char*)Ad.data(), Ad.msize()); +<<<<<<< HEAD std::vector clsFlagChar(com_mod.ris.clsFlg.size()); bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { @@ -136,6 +137,10 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< for (int i = 0; i < com_mod.nUris; i++) { com_mod.uris[i].cnt = urisCnt(i); com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} +======= + // [HZ] [TODO] How to read std::vector from file? + // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); +>>>>>>> 1da8805 (Adding explicit RIS code.) } else { bin_file.read((char*)Ad.data(), Ad.msize()); } @@ -147,6 +152,7 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bin_file.read((char*)Xion.data(), Xion.msize()); bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); } else if (risFlag) { +<<<<<<< HEAD std::vector clsFlagChar(com_mod.ris.clsFlg.size()); bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { @@ -159,6 +165,10 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< for (int i = 0; i < com_mod.nUris; i++) { com_mod.uris[i].cnt = urisCnt(i); com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} +======= + // [HZ] [TODO] How to read std::vector from file? + // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); +>>>>>>> 1da8805 (Adding explicit RIS code.) } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao, Do } @@ -170,6 +180,7 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< if (cepEq) { bin_file.read((char*)Xion.data(), Xion.msize()); } else if (risFlag) { +<<<<<<< HEAD std::vector clsFlagChar(com_mod.ris.clsFlg.size()); bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { @@ -182,6 +193,10 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< for (int i = 0; i < com_mod.nUris; i++) { com_mod.uris[i].cnt = urisCnt(i); com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} +======= + // [HZ] [TODO] How to read std::vector from file? + // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); +>>>>>>> 1da8805 (Adding explicit RIS code.) } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao } @@ -512,9 +527,12 @@ void initialize(Simulation* simulation, Vector& timeP) if (com_mod.risFlag) { i = i + com_mod.ris.nbrRIS; } +<<<<<<< HEAD if (com_mod.urisFlag) { i = i + com_mod.nUris * 2; } +======= +>>>>>>> 1da8805 (Adding explicit RIS code.) i = sizeof(int)*(1+com_mod.stamp.size()) + sizeof(double)*(2 + com_mod.nEq + com_mod.cplBC.nX + i*com_mod.tnNo); From 0d04e4751678d35ba276bf148f4e769e02aec67b Mon Sep 17 00:00:00 2001 From: Han Zhao Date: Mon, 3 Mar 2025 10:43:16 -0800 Subject: [PATCH 03/19] Updating output and restrat functions for explicit RIS Resolving conflicts --- Code/Source/solver/initialize.cpp | 49 ------------------------------- 1 file changed, 49 deletions(-) diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index b99b7efd4..bc326df43 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -123,24 +123,8 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); } else if (risFlag) { bin_file.read((char*)Ad.data(), Ad.msize()); -<<<<<<< HEAD - std::vector clsFlagChar(com_mod.ris.clsFlg.size()); - bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); - for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { - com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} - } else if (urisFlag) { - bin_file.read((char*)Ad.data(), Ad.msize()); - Vector urisCnt(com_mod.nUris); - bin_file.read((char*)urisCnt.data(), urisCnt.msize()); - std::vector urisClsFlagChar(com_mod.nUris); - bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); - for (int i = 0; i < com_mod.nUris; i++) { - com_mod.uris[i].cnt = urisCnt(i); - com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} -======= // [HZ] [TODO] How to read std::vector from file? // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); ->>>>>>> 1da8805 (Adding explicit RIS code.) } else { bin_file.read((char*)Ad.data(), Ad.msize()); } @@ -152,23 +136,8 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bin_file.read((char*)Xion.data(), Xion.msize()); bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); } else if (risFlag) { -<<<<<<< HEAD - std::vector clsFlagChar(com_mod.ris.clsFlg.size()); - bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); - for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { - com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} - } else if (urisFlag) { - Vector urisCnt(com_mod.nUris); - bin_file.read((char*)urisCnt.data(), urisCnt.msize()); - std::vector urisClsFlagChar(com_mod.nUris); - bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); - for (int i = 0; i < com_mod.nUris; i++) { - com_mod.uris[i].cnt = urisCnt(i); - com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} -======= // [HZ] [TODO] How to read std::vector from file? // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); ->>>>>>> 1da8805 (Adding explicit RIS code.) } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao, Do } @@ -180,23 +149,8 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< if (cepEq) { bin_file.read((char*)Xion.data(), Xion.msize()); } else if (risFlag) { -<<<<<<< HEAD - std::vector clsFlagChar(com_mod.ris.clsFlg.size()); - bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); - for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { - com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} - } else if (urisFlag) { - Vector urisCnt(com_mod.nUris); - bin_file.read((char*)urisCnt.data(), urisCnt.msize()); - std::vector urisClsFlagChar(com_mod.nUris); - bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); - for (int i = 0; i < com_mod.nUris; i++) { - com_mod.uris[i].cnt = urisCnt(i); - com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} -======= // [HZ] [TODO] How to read std::vector from file? // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); ->>>>>>> 1da8805 (Adding explicit RIS code.) } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao } @@ -527,12 +481,9 @@ void initialize(Simulation* simulation, Vector& timeP) if (com_mod.risFlag) { i = i + com_mod.ris.nbrRIS; } -<<<<<<< HEAD if (com_mod.urisFlag) { i = i + com_mod.nUris * 2; } -======= ->>>>>>> 1da8805 (Adding explicit RIS code.) i = sizeof(int)*(1+com_mod.stamp.size()) + sizeof(double)*(2 + com_mod.nEq + com_mod.cplBC.nX + i*com_mod.tnNo); From 8beff2c865fb1e31e4de24e75d87df6a477f7ba0 Mon Sep 17 00:00:00 2001 From: Han Zhao Date: Thu, 13 Mar 2025 09:41:28 -0700 Subject: [PATCH 04/19] Adding implicit RIS functions for serial simulation. Resovling conflicts --- Code/Source/solver/Parameters.cpp | 172 +++++++++++++++++++++++++++--- Code/Source/solver/Parameters.h | 29 +---- Code/Source/solver/fluid.cpp | 6 ++ Code/Source/solver/initialize.cpp | 43 ++++++-- Code/Source/solver/vtk_xml.cpp | 5 + 5 files changed, 206 insertions(+), 49 deletions(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index cd8f20694..be2bc883b 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2607,13 +2607,16 @@ void PrecomputedSolutionParameters::set_values(tinyxml2::XMLElement* xml_elem) } ////////////////////////////////////////////////////////// -// RISProjectionParameters // +// P r o j e c t i o n P a r a m e t e r s // ////////////////////////////////////////////////////////// +// The ProjectionParameters class stores parameters for the +// 'Add_projection' XML element used for fluid-structure interaction simulations. + /// @brief Define the XML element name for mesh parameters. -const std::string RISProjectionParameters::xml_element_name_ = "Add_RIS_projection"; +const std::string ProjectionParameters::xml_element_name_ = "Add_projection"; -RISProjectionParameters::RISProjectionParameters() +ProjectionParameters::ProjectionParameters() { // A parameter that must be defined. bool required = true; @@ -2621,16 +2624,15 @@ RISProjectionParameters::RISProjectionParameters() name = Parameter("name", "", required); set_parameter("Project_from_face", "", required, project_from_face); - set_parameter("Resistance", 1.e6, !required, resistance); set_parameter("Projection_tolerance", 0.0, !required, projection_tolerance); } -void RISProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) +void ProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) { using namespace tinyxml2; std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; - // Get the 'type' from the element. + // Get the 'type' from the element. const char* sname; auto result = xml_elem->QueryStringAttribute("name", &sname); if (sname == nullptr) { @@ -2642,23 +2644,19 @@ void RISProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) using std::placeholders::_2; std::function ftpr = - std::bind( &RISProjectionParameters::set_parameter_value, *this, _1, _2); + std::bind( &ProjectionParameters::set_parameter_value, *this, _1, _2); xml_util_set_parameters(ftpr, xml_elem, error_msg); } - ////////////////////////////////////////////////////////// -// P r o j e c t i o n P a r a m e t e r s // +// RISProjectionParameters // ////////////////////////////////////////////////////////// -// The ProjectionParameters class stores parameters for the -// 'Add_projection' XML element used for fluid-structure interaction simulations. - /// @brief Define the XML element name for mesh parameters. -const std::string ProjectionParameters::xml_element_name_ = "Add_projection"; +const std::string RISProjectionParameters::xml_element_name_ = "Add_RIS_projection"; -ProjectionParameters::ProjectionParameters() +RISProjectionParameters::RISProjectionParameters() { // A parameter that must be defined. bool required = true; @@ -2666,15 +2664,16 @@ ProjectionParameters::ProjectionParameters() name = Parameter("name", "", required); set_parameter("Project_from_face", "", required, project_from_face); + set_parameter("Resistance", 1.e6, !required, resistance); set_parameter("Projection_tolerance", 0.0, !required, projection_tolerance); } -void ProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) +void RISProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) { using namespace tinyxml2; std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; - // Get the 'type' from the element. + // Get the 'type' from the element. const char* sname; auto result = xml_elem->QueryStringAttribute("name", &sname); if (sname == nullptr) { @@ -2686,11 +2685,150 @@ void ProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) using std::placeholders::_2; std::function ftpr = - std::bind( &ProjectionParameters::set_parameter_value, *this, _1, _2); + std::bind( &RISProjectionParameters::set_parameter_value, *this, _1, _2); xml_util_set_parameters(ftpr, xml_elem, error_msg); } + +////////////////////////////////////////////////////////// +// URIS Mesh Parameters // +////////////////////////////////////////////////////////// +// [HZ] implemente URIS parameters here + +// Process parameters for the 'Add_URIS_mesh' XML element used for defining URIS mesh elements. + +/// @brief Define the XML element name for mesh parameters. +const std::string URISMeshParameters::xml_element_name_ = "Add_URIS_mesh"; + +URISMeshParameters::URISMeshParameters() +{ + bool required = true; + + // Mesh name from Add_mesh element. + name = Parameter("name", "", required); + + // Parameters under Add_mesh element. + // + set_parameter("Mesh_scale_factor", 1.0, !required, mesh_scale_factor); + set_parameter("Thickness", 1.0, !required, thickness); + set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); + set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); +} + +void URISMeshParameters::print_parameters() +{ + std::cout << std::endl; + std::cout << "---------------" << std::endl; + std::cout << "URIS Mesh Parameters" << std::endl; + std::cout << "---------------" << std::endl; + std::cout << name.name() << ": " << name.value() << std::endl; + + auto params_name_value = get_parameter_list(); + for (auto& [ key, value ] : params_name_value) { + std::cout << key << ": " << value << std::endl; + } + + for (auto& face : URIS_face_parameters) { + face->print_parameters(); + } +} + +void URISMeshParameters::set_values(tinyxml2::XMLElement* mesh_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + auto item = mesh_elem->FirstChildElement(); + + while (item != nullptr) { + auto name = std::string(item->Value()); + + // Add_face sub-element. + if (name == URISFaceParameters::xml_element_name_) { + auto URIS_face_params = new URISFaceParameters(); + URIS_face_params->set_values(item); + URIS_face_parameters.push_back(URIS_face_params); + } else if (item->GetText() != nullptr) { + auto value = item->GetText(); + try { + set_parameter_value(name, value); + } catch (const std::bad_function_call& exception) { + throw std::runtime_error(error_msg + name + "'."); + } + } else { + throw std::runtime_error(error_msg + name + "'."); + } + + item = item->NextSiblingElement(); + } +} + + +////////////////////////////////////////////////////////// +// URIS Face Parameters // +////////////////////////////////////////////////////////// + +/// @brief Process parameters for the 'Add_URIS_face' XML element. +/// +/// Define the XML element name for face parameters. +const std::string URISFaceParameters::xml_element_name_ = "Add_URIS_face"; + +URISFaceParameters::URISFaceParameters() +{ + set_xml_element_name(xml_element_name_); + + // A parameter that must be defined. + bool required = true; + + name = Parameter("name", "", required); + + set_parameter("Face_file_path", "", !required, face_file_path); + set_parameter("Open_motion_file_path", "", !required, open_motion_file_path); + set_parameter("Close_motion_file_path", "", !required, close_motion_file_path); + + // set_parameter("End_nodes_face_file_path", "", !required, end_nodes_face_file_path); + // set_parameter("Quadrature_modifier_TRI3", (2.0/3.0), !required, quadrature_modifier_TRI3); +} + +void URISFaceParameters::print_parameters() +{ + std::cout << std::endl; + std::cout << "---------------" << std::endl; + std::cout << "URIS Face Parameters" << std::endl; + std::cout << "---------------" << std::endl; + std::cout << name.name() << ": " << name.value() << std::endl; + std::cout << face_file_path.name() << ": " << face_file_path.value() << std::endl; +} + +void URISFaceParameters::set_values(tinyxml2::XMLElement* face_elem) +{ + using namespace tinyxml2; + + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + const char* face_name; + auto result = face_elem->QueryStringAttribute("name", &face_name); + name.set(std::string(face_name)); + auto item = face_elem->FirstChildElement(); + + while (item != nullptr) { + auto name = std::string(item->Value()); + auto value = item->GetText(); + + if (value == nullptr) { + throw std::runtime_error(error_msg + name + "'."); + } + + try { + set_parameter_value(name, value); + } catch (const std::bad_function_call& exception) { + throw std::runtime_error(error_msg + name + "'."); + } + + item = item->NextSiblingElement(); + } +} + + ////////////////////////////////////////////////////////// // RISProjectionParameters // ////////////////////////////////////////////////////////// diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 9aa79bae1..fd176330a 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -864,30 +864,6 @@ class PrecomputedSolutionParameters: public ParameterLists Parameter use_precomputed_solution; }; -/// @brief The RISProjectionParameters class stores parameters for the -/// 'Add_RIS_projection' XML element used for RIS valve simulations. -/// \code {.xml} -/// -/// right_ris -/// 1.e6 -/// -/// \endcode -class RISProjectionParameters : public ParameterLists -{ - public: - RISProjectionParameters(); - - void set_values(tinyxml2::XMLElement* xml_elem); - - static const std::string xml_element_name_; - - Parameter name; - - Parameter project_from_face; - Parameter resistance; - Parameter projection_tolerance; -}; - /// @brief The ProjectionParameters class stores parameters for the /// 'Add_projection' XML element used for fluid-structure interaction /// simulations. @@ -1660,19 +1636,20 @@ class Parameters { void set_equation_values(tinyxml2::XMLElement* root_element); void set_mesh_values(tinyxml2::XMLElement* root_element); void set_precomputed_solution_values(tinyxml2::XMLElement* root_element); - void set_RIS_projection_values(tinyxml2::XMLElement* root_element); void set_projection_values(tinyxml2::XMLElement* root_element); void set_svzerodsolver_interface_values(tinyxml2::XMLElement* root_element); void set_RIS_projection_values(tinyxml2::XMLElement* root_element); void set_URIS_mesh_values(tinyxml2::XMLElement* root_element); + void set_RIS_projection_values(tinyxml2::XMLElement* root_element); + void set_URIS_mesh_values(tinyxml2::XMLElement* root_element); + // Objects representing each parameter section of XML file. ContactParameters contact_parameters; GeneralSimulationParameters general_simulation_parameters; std::vector mesh_parameters; std::vector equation_parameters; - std::vector RIS_projection_parameters; std::vector projection_parameters; PrecomputedSolutionParameters precomputed_solution_parameters; diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index b064ea762..0fb9e43dc 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -1472,7 +1472,10 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; +<<<<<<< HEAD if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} +======= +>>>>>>> 3257cf7 (Adding implicit RIS functions for serial simulation.) } double rho = dmn.prop[PhysicalProperyType::fluid_density]; @@ -1802,7 +1805,10 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; +<<<<<<< HEAD if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} +======= +>>>>>>> 3257cf7 (Adding implicit RIS functions for serial simulation.) } double ctM = 1.0; diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index bc326df43..7f1e1ae84 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -123,8 +123,19 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); } else if (risFlag) { bin_file.read((char*)Ad.data(), Ad.msize()); - // [HZ] [TODO] How to read std::vector from file? - // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} + } else if (urisFlag) { + bin_file.read((char*)Ad.data(), Ad.msize()); + Vector urisCnt(com_mod.nUris); + bin_file.read((char*)urisCnt.data(), urisCnt.msize()); + std::vector urisClsFlagChar(com_mod.nUris); + bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.nUris; i++) { + com_mod.uris[i].cnt = urisCnt(i); + com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} } else { bin_file.read((char*)Ad.data(), Ad.msize()); } @@ -136,8 +147,18 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bin_file.read((char*)Xion.data(), Xion.msize()); bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); } else if (risFlag) { - // [HZ] [TODO] How to read std::vector from file? - // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} + } else if (urisFlag) { + Vector urisCnt(com_mod.nUris); + bin_file.read((char*)urisCnt.data(), urisCnt.msize()); + std::vector urisClsFlagChar(com_mod.nUris); + bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.nUris; i++) { + com_mod.uris[i].cnt = urisCnt(i); + com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao, Do } @@ -149,8 +170,18 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< if (cepEq) { bin_file.read((char*)Xion.data(), Xion.msize()); } else if (risFlag) { - // [HZ] [TODO] How to read std::vector from file? - // bin_file.read((char*)com_mod.RIS.clsFlg.data(), com_mod.RIS.clsFlg.msize()); + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} + } else if (urisFlag) { + Vector urisCnt(com_mod.nUris); + bin_file.read((char*)urisCnt.data(), urisCnt.msize()); + std::vector urisClsFlagChar(com_mod.nUris); + bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.nUris; i++) { + com_mod.uris[i].cnt = urisCnt(i); + com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao } diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 8cb79d9b0..bd96c9be0 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1379,6 +1379,8 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, const Arrayset_point_data(outNames[iOut], tmpV); } From 1a5b2adff8262a867659087505eec907530579a0 Mon Sep 17 00:00:00 2001 From: Han Zhao Date: Fri, 28 Mar 2025 09:07:51 -0700 Subject: [PATCH 05/19] Parallelization for implicit RIS code. --- Code/Source/solver/distribute.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 932d301ea..4cc9aa2f9 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -340,11 +340,17 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &simulation->cep_mod.cepEq); cm.bcast(cm_mod, &com_mod.risFlag); +<<<<<<< HEAD cm.bcast(cm_mod, &com_mod.ris0DFlag); cm.bcast(cm_mod, &com_mod.urisFlag); cm.bcast(cm_mod, &com_mod.urisActFlag); cm.bcast(cm_mod, &com_mod.urisRes); cm.bcast(cm_mod, &com_mod.urisResClose); +======= + cm.bcast(cm_mod, &com_mod.urisFlag); + cm.bcast(cm_mod, &com_mod.urisActFlag); + cm.bcast(cm_mod, &com_mod.urisRes); +>>>>>>> 582b719 (Parallelization for implicit RIS code.) cm.bcast(cm_mod, &com_mod.usePrecomp); if (com_mod.rmsh.isReqd) { auto& rmsh = com_mod.rmsh; @@ -1168,7 +1174,10 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].nFa); cm.bcast(cm_mod, &uris[iUris].sdf_default); cm.bcast(cm_mod, &uris[iUris].sdf_deps); +<<<<<<< HEAD cm.bcast(cm_mod, &uris[iUris].sdf_deps_close); +======= +>>>>>>> 582b719 (Parallelization for implicit RIS code.) cm.bcast(cm_mod, &uris[iUris].clsFlg); cm.bcast(cm_mod, &uris[iUris].cnt); cm.bcast(cm_mod, &uris[iUris].scF); From ac46ddac49ea236c8d5252561202aadfc851b79f Mon Sep 17 00:00:00 2001 From: Han Zhao Date: Wed, 2 Apr 2025 11:05:45 -0700 Subject: [PATCH 06/19] Adding RIS0D functions. Resolving conflicts --- Code/Source/solver/ComMod.h | 3 +++ Code/Source/solver/distribute.cpp | 4 ---- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index a688145b6..1fdf5c530 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -171,6 +171,9 @@ class bcType // RIS0D: resistance double resistance = 0.0; + // RIS0D: resistance + double resistance = 0.0; + // Penalty parameters for weakly applied Dir BC Vector tauB{0.0, 0.0}; //double tauB[2]; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 4cc9aa2f9..bd65088c2 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -340,7 +340,6 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &simulation->cep_mod.cepEq); cm.bcast(cm_mod, &com_mod.risFlag); -<<<<<<< HEAD cm.bcast(cm_mod, &com_mod.ris0DFlag); cm.bcast(cm_mod, &com_mod.urisFlag); cm.bcast(cm_mod, &com_mod.urisActFlag); @@ -1174,10 +1173,7 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].nFa); cm.bcast(cm_mod, &uris[iUris].sdf_default); cm.bcast(cm_mod, &uris[iUris].sdf_deps); -<<<<<<< HEAD cm.bcast(cm_mod, &uris[iUris].sdf_deps_close); -======= ->>>>>>> 582b719 (Parallelization for implicit RIS code.) cm.bcast(cm_mod, &uris[iUris].clsFlg); cm.bcast(cm_mod, &uris[iUris].cnt); cm.bcast(cm_mod, &uris[iUris].scF); From 1056e82ea5f314ef5a9f7feb244226624e8e8675 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Mon, 9 Jun 2025 10:14:32 -0700 Subject: [PATCH 07/19] Improve the SDF function and add valve thickness and ressitance specification. Resolving conflicts --- Code/Source/solver/Parameters.cpp | 4 +++- Code/Source/solver/fluid.cpp | 8 ++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index be2bc883b..c69446482 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2711,7 +2711,9 @@ URISMeshParameters::URISMeshParameters() // Parameters under Add_mesh element. // set_parameter("Mesh_scale_factor", 1.0, !required, mesh_scale_factor); - set_parameter("Thickness", 1.0, !required, thickness); + set_parameter("Thickness", 0.04, !required, thickness); + set_parameter("Close_thickness", 0.25, !required, close_thickness); + set_parameter("Resistance", 1.0e5, !required, resistance); set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); } diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index 0fb9e43dc..d57bb15a2 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -1472,10 +1472,14 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; +<<<<<<< HEAD <<<<<<< HEAD if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} ======= >>>>>>> 3257cf7 (Adding implicit RIS functions for serial simulation.) +======= + if (com_mod.uris[0].clsFlg) {Res = 1.0e5;} +>>>>>>> f6358a3 (Improve the SDF function and add valve thickness and ressitance specification.) } double rho = dmn.prop[PhysicalProperyType::fluid_density]; @@ -1805,10 +1809,14 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; +<<<<<<< HEAD <<<<<<< HEAD if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} ======= >>>>>>> 3257cf7 (Adding implicit RIS functions for serial simulation.) +======= + if (com_mod.uris[0].clsFlg) {Res = 1.0e5;} +>>>>>>> f6358a3 (Improve the SDF function and add valve thickness and ressitance specification.) } double ctM = 1.0; From b1e7a095bdbaf9285372eced5523979eb021bbcc Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Tue, 12 Aug 2025 22:24:37 +0800 Subject: [PATCH 08/19] Adding option for using customized resistance value whenthe RIS valve is closed. --- Code/Source/solver/Parameters.cpp | 3 ++- Code/Source/solver/distribute.cpp | 5 +---- Code/Source/solver/fluid.cpp | 14 -------------- 3 files changed, 3 insertions(+), 19 deletions(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index c69446482..55211ca37 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2712,8 +2712,9 @@ URISMeshParameters::URISMeshParameters() // set_parameter("Mesh_scale_factor", 1.0, !required, mesh_scale_factor); set_parameter("Thickness", 0.04, !required, thickness); - set_parameter("Close_thickness", 0.25, !required, close_thickness); + set_parameter("Closed_thickness", 0.25, !required, close_thickness); set_parameter("Resistance", 1.0e5, !required, resistance); + set_parameter("Closed_resistance", 1.0e5, !required, resistance_close); set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); } diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index bd65088c2..50b99416b 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -346,10 +346,7 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.urisRes); cm.bcast(cm_mod, &com_mod.urisResClose); ======= - cm.bcast(cm_mod, &com_mod.urisFlag); - cm.bcast(cm_mod, &com_mod.urisActFlag); - cm.bcast(cm_mod, &com_mod.urisRes); ->>>>>>> 582b719 (Parallelization for implicit RIS code.) +>>>>>>> edd5c69 (Adding option for using customized resistance value whenthe RIS valve is closed.) cm.bcast(cm_mod, &com_mod.usePrecomp); if (com_mod.rmsh.isReqd) { auto& rmsh = com_mod.rmsh; diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index d57bb15a2..298e1e8c8 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -1472,14 +1472,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; -<<<<<<< HEAD -<<<<<<< HEAD - if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} -======= ->>>>>>> 3257cf7 (Adding implicit RIS functions for serial simulation.) -======= if (com_mod.uris[0].clsFlg) {Res = 1.0e5;} ->>>>>>> f6358a3 (Improve the SDF function and add valve thickness and ressitance specification.) } double rho = dmn.prop[PhysicalProperyType::fluid_density]; @@ -1809,14 +1802,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; -<<<<<<< HEAD -<<<<<<< HEAD - if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} -======= ->>>>>>> 3257cf7 (Adding implicit RIS functions for serial simulation.) -======= if (com_mod.uris[0].clsFlg) {Res = 1.0e5;} ->>>>>>> f6358a3 (Improve the SDF function and add valve thickness and ressitance specification.) } double ctM = 1.0; From f76aa865623aa85d00edabfe09535f1701bbbbb8 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Wed, 13 Aug 2025 09:18:25 +0800 Subject: [PATCH 09/19] Resolving a double definition of bcast for int array --- Code/Source/solver/CmMod.h | 1 - 1 file changed, 1 deletion(-) diff --git a/Code/Source/solver/CmMod.h b/Code/Source/solver/CmMod.h index 4e355d83d..36ac60e18 100644 --- a/Code/Source/solver/CmMod.h +++ b/Code/Source/solver/CmMod.h @@ -79,7 +79,6 @@ class cmType { void bcast(const CmMod& cm_mod, double* data) const; void bcast(const CmMod& cm_mod, Vector& data, const std::string& name="") const; void bcast(const CmMod& cm_mod, Array& data, const std::string& name="") const; - void bcast(const CmMod& cm_mod, Array& data, const std::string& name="") const; void bcast(const CmMod& cm_mod, int* data) const; void bcast(const CmMod& cm_mod, Vector& data) const; From f9bfaa26eee638ce5082aa3803610c2bf6e66c44 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Sat, 20 Sep 2025 18:54:56 -0700 Subject: [PATCH 10/19] Adding functions to support fixed surface scaffold in unfitted RIS valve simulation. --- Code/Source/solver/ComMod.h | 8 +- Code/Source/solver/Parameters.cpp | 3 +- Code/Source/solver/Parameters.h | 1 + Code/Source/solver/distribute.cpp | 26 ++++++ Code/Source/solver/fluid.cpp | 17 +++- Code/Source/solver/fsi.cpp | 18 +++- Code/Source/solver/uris.cpp | 149 ++++++++++++++++++++++++++++-- Code/Source/solver/vtk_xml.cpp | 21 +++++ 8 files changed, 230 insertions(+), 13 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 1fdf5c530..fa144e8e3 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1485,7 +1485,7 @@ class urisType // Iteration count. int cnt = 1000000; - // URIS: signed distance function of each node to the uris (1D array). + // URIS: signed distance function of each node to the uris valves (1D array). Vector sdf; // Mesh scale factor. @@ -1510,6 +1510,12 @@ class urisType // IB meshes std::vector msh; + // Scaffold mesh + bool scaffold_flag = false; + mshType scaffold_mesh; + + Vector sdf_scaffold; + }; /// @brief The ComMod class duplicates the data structures in the Fortran COMMOD module diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 55211ca37..75d229d2f 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2711,12 +2711,13 @@ URISMeshParameters::URISMeshParameters() // Parameters under Add_mesh element. // set_parameter("Mesh_scale_factor", 1.0, !required, mesh_scale_factor); - set_parameter("Thickness", 0.04, !required, thickness); + set_parameter("Thickness", 0.25, !required, thickness); set_parameter("Closed_thickness", 0.25, !required, close_thickness); set_parameter("Resistance", 1.0e5, !required, resistance); set_parameter("Closed_resistance", 1.0e5, !required, resistance_close); set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); + set_parameter("Scaffold_file_path", "", !required, scaffold_file_path); } void URISMeshParameters::print_parameters() diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index fd176330a..6d26437e3 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1613,6 +1613,7 @@ class URISMeshParameters : public ParameterLists Parameter resistance_close; // Resistance of the valve when it is closed Parameter valve_starts_as_closed; // Whether the valve starts as closed Parameter positive_flow_normal_file_path; // File path for the positive flow normal + Parameter scaffold_file_path; // File path for the scaffold }; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 50b99416b..393bc6022 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1175,6 +1175,32 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].cnt); cm.bcast(cm_mod, &uris[iUris].scF); cm.bcast(cm_mod, uris[iUris].nrm); + + cm.bcast(cm_mod, &uris[iUris].scaffold_flag); + if (uris[iUris].scaffold_flag) { + cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.lShl); + cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.nEl); + cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.gnEl); + cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.eNoN); + cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.nNo); + cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.gnNo); + } + } + + if (cm.slv(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (uris[iUris].scaffold_flag) { + uris[iUris].scaffold_mesh.x.resize(com_mod.nsd, uris[iUris].scaffold_mesh.gnNo); + uris[iUris].scaffold_mesh.IEN.resize(uris[iUris].scaffold_mesh.eNoN, uris[iUris].scaffold_mesh.gnEl); + } + } + } + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (uris[iUris].scaffold_flag) { + cm.bcast(cm_mod, uris[iUris].scaffold_mesh.x); + cm.bcast(cm_mod, uris[iUris].scaffold_mesh.IEN); + } } std::vector> lM_gN_flat(com_mod.nUris); diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index 298e1e8c8..646c8d0c0 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -651,17 +651,22 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag // Plot the coordinates of the quad point in the current configuration if (com_mod.urisFlag) { Vector distSrf(com_mod.nUris); + Vector distSrf_scaffold(com_mod.nUris); distSrf = 0.0; + distSrf_scaffold = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < com_mod.nUris; iUris++) { distSrf(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + if (com_mod.uris[iUris].scaffold_flag) { + distSrf_scaffold(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf_scaffold(Ac)); + } } } DDir = 0.0; + double sdf_deps_temp = 0; double DDirTmp = 0.0; - double sdf_deps_temp = 0.0; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { if (com_mod.uris[iUris].clsFlg) { sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; @@ -673,8 +678,16 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag (2*sdf_deps_temp*sdf_deps_temp); if (DDirTmp > DDir) {DDir = DDirTmp;} } - } + if (com_mod.uris[iUris].scaffold_flag) { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + } + if (distSrf_scaffold(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + if (DDirTmp > DDir) {DDir = DDirTmp;} + } + } if (!com_mod.urisActFlag) {DDir = 0.0;} } diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 18a207d7e..990eb6bcc 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -189,11 +189,16 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // Plot the coordinates of the quad point in the current configuration if (com_mod.urisFlag) { Vector distSrf(com_mod.nUris); + Vector distSrf_scaffold(com_mod.nUris); distSrf = 0.0; + distSrf_scaffold = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < com_mod.nUris; iUris++) { distSrf(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + if (com_mod.uris[iUris].scaffold_flag) { + distSrf_scaffold(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf_scaffold(Ac)); + } } } @@ -217,14 +222,21 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar (2*sdf_deps_temp*sdf_deps_temp); if (DDirTmp > DDir) {DDir = DDirTmp;} } + + if (com_mod.uris[iUris].scaffold_flag) { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + } + if (distSrf_scaffold(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + if (DDirTmp > DDir) {DDir = DDirTmp;} + } + } if (!com_mod.urisActFlag) {DDir = 0.0;} - - // std::cout << "===== DDir: " << DDir << std::endl; } - if (nsd == 3) { switch (cPhys) { case Equation_fluid: { diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index cadbb34f6..66b9e796e 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -518,9 +518,46 @@ void uris_read_msh(Simulation* simulation) { } file_stream.close(); - uris_obj.sdf_default = uris_obj.sdf_default * uris_obj.scF; - uris_obj.sdf_deps = param->thickness() * uris_obj.scF; - uris_obj.sdf_deps_close = param->close_thickness() * uris_obj.scF; + std::string scaffold_file_path = param->scaffold_file_path(); + if (scaffold_file_path != "") { + uris_obj.scaffold_flag = true; + uris_obj.scaffold_mesh.lShl = true; + vtk_xml::read_vtu(scaffold_file_path, uris_obj.scaffold_mesh); + nn::select_ele(simulation->com_mod, uris_obj.scaffold_mesh); + read_msh_ns::check_ien(simulation, uris_obj.scaffold_mesh); + + int b = 0; + auto& scaffold_mesh = uris_obj.scaffold_mesh; + scaffold_mesh.nNo = scaffold_mesh.gnNo; + scaffold_mesh.gN.resize(scaffold_mesh.nNo); + scaffold_mesh.gN = 0; + scaffold_mesh.lN.resize(scaffold_mesh.nNo); + scaffold_mesh.lN = 0; + for (int a = 0; a < scaffold_mesh.nNo; a++) { + scaffold_mesh.gN(a) = b; + scaffold_mesh.lN(b) = a; + b++; + } + + // Remap msh%gIEN array + scaffold_mesh.nEl = scaffold_mesh.gnEl; + scaffold_mesh.IEN.resize(scaffold_mesh.eNoN, scaffold_mesh.nEl); + for (int e = 0; e < scaffold_mesh.nEl; e++) { + for (int a = 0; a < scaffold_mesh.eNoN; a++) { + int Ac = scaffold_mesh.gIEN(a,e); + Ac = scaffold_mesh.gN(Ac); + scaffold_mesh.IEN(a,e) = Ac; + } + } + scaffold_mesh.gIEN.clear(); + + std::cout << "Scaffold mesh is included for: " << uris_obj.name << std::endl; + std::cout << "Scaffold mesh nodes: " << uris_obj.scaffold_mesh.gnNo << std::endl; + std::cout << "Scaffold mesh elements: " << uris_obj.scaffold_mesh.gnEl << std::endl; + } + + uris_obj.sdf_deps = param->thickness(); + uris_obj.sdf_deps_close = param->close_thickness(); uris_obj.clsFlg = param->valve_starts_as_closed(); // uris_obj.tnNo = 0; @@ -714,6 +751,10 @@ void uris_read_msh(Simulation* simulation) { } } mesh.gIEN.clear(); + + std::cout << "Number of uris elements nel: " << mesh.nEl << std::endl; + std::cout << "URIS mesh IEN size: " << mesh.IEN.nrows() << ", " << mesh.IEN.ncols() << std::endl; + } if (uris_obj.nFa > 0) { @@ -896,6 +937,8 @@ void uris_calc_sdf(ComMod& com_mod) { for (int iUris = 0; iUris < nUris; iUris++) { // We need to check if the valve needs to move auto& uris_obj = uris[iUris]; + + // Compute the SDF for the uris valves int cnt = 0; if (!uris_obj.clsFlg) { cnt = std::min(uris_obj.cnt, uris_obj.DxOpen.nrows()); @@ -923,8 +966,8 @@ void uris_calc_sdf(ComMod& com_mod) { max_eNoN = mesh.eNoN; } } - Array lX(nsd, max_eNoN); + // if (!uris_obj.sdf.allocated()) { if (uris_obj.sdf.size() <= 0) { uris_obj.sdf.resize(com_mod.tnNo); @@ -987,7 +1030,9 @@ void uris_calc_sdf(ComMod& com_mod) { int jM = -1; Vector xb(nsd); for (int iM = 0; iM < uris_obj.nFa; iM++) { + auto& mesh = uris_obj.msh[iM]; + for (int e = 0; e < mesh.nEl; e++) { xb = 0.0; for (int a = 0; a < mesh.eNoN; a++) { @@ -1015,6 +1060,7 @@ void uris_calc_sdf(ComMod& com_mod) { // We also need to compute the sign (above or below the valve). // Compute the element normal + auto& mesh = uris_obj.msh[jM]; xXi = 0.0; lX = 0.0; @@ -1069,6 +1115,97 @@ void uris_calc_sdf(ComMod& com_mod) { } } } + + + // Compute the SDF for the scaffold mesh + for (int iUris = 0; iUris < nUris; iUris++) { + auto& uris_obj = uris[iUris]; + + if (!uris_obj.scaffold_flag || uris_obj.sdf_scaffold.size() > 0) { + continue; + } + + auto& scaffold_mesh = uris_obj.scaffold_mesh; + + if (uris_obj.sdf_scaffold.size() <= 0) { + uris_obj.sdf_scaffold.resize(com_mod.tnNo); + uris_obj.sdf_scaffold = 0.0; + } + uris_obj.sdf_scaffold = uris_obj.sdf_default; + + Vector minb_scaf(nsd); + Vector maxb_scaf(nsd); + Vector extra_scaf(nsd); + for (int i = 0; i < nsd; i++) { + minb_scaf(i) = std::numeric_limits::max(); + maxb_scaf(i) = std::numeric_limits::lowest(); + } + + // For each coordinate dimension, find the minimum and maximum in uris_obj.x. + double extra_scaf_val = 0.1; // [HZ] The BBox is 10% larger than the actual valve, default is 0.1 + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < scaffold_mesh.x.ncols(); j++) { + double val = scaffold_mesh.x(i,j); + if (val < minb_scaf(i)) + minb_scaf(i) = val; + if (val > maxb_scaf(i)) + maxb_scaf(i) = val; + } + extra_scaf(i) = (maxb_scaf(i) - minb_scaf(i)) * extra_scaf_val; + } + + // The SDF is computed on the reference configuration, which + // means that the valves will be morphed based on the fluid mesh + // motion. If the fluid mesh stretches near the valve, the valve + // leaflets will also be streched. Note that + // this is a simplifying assumption. + Vector xp_scaf(nsd); + for (int ca = 0; ca < com_mod.tnNo; ca++) { + // std::cout << "URIS scaffold SDF index: " << ca << std::endl; + double minS_scaf = std::numeric_limits::max(); + for (int i = 0; i < nsd; i++) { + xp_scaf(i) = com_mod.x(i,ca); + } + // Is the node inside the BBox? + bool inside = true; + for (int i = 0; i < nsd; i++) { + if (xp_scaf(i) < (minb_scaf(i) - extra_scaf(i)) || xp_scaf(i) > (maxb_scaf(i) + extra_scaf(i))) { + inside = false; + break; + } + } + + if (inside) { + // This point is inside the BBox + // Find the closest URIS face centroid + Vector xb_scaf(nsd); + + for (int e = 0; e < scaffold_mesh.nEl; e++) { + xb_scaf = 0.0; + for (int a = 0; a < scaffold_mesh.eNoN; a++) { + int Ac = scaffold_mesh.IEN(a,e); + for (int i = 0; i < nsd; i++) { + xb_scaf(i) += scaffold_mesh.x(i,Ac); + } + } + for (int i = 0; i < nsd; i++) { + xb_scaf(i) /= static_cast(scaffold_mesh.eNoN); + } + double dS_scaf = 0.0; + for (int i = 0; i < nsd; i++) { + dS_scaf += (xp_scaf[i] - xb_scaf[i]) * (xp_scaf[i] - xb_scaf[i]); + } + dS_scaf = std::sqrt(dS_scaf); + + if (dS_scaf < minS_scaf) { + minS_scaf = dS_scaf; + } + } + uris_obj.sdf_scaffold[ca] = minS_scaf; + } + } + } + } @@ -1177,8 +1314,8 @@ int in_poly(Vector& P, Array& P1, bool ext) { /// @brief Chech if a point is on the same side of anotehr point wrt a triangle in 3D int same_side(Vector& v1, Vector& v2, Vector& v3, Vector& v4, Vector& p, bool ext) { - #define n_dbg_read_sv - #ifdef dbg_read_sv + #define n_dbg_same_side + #ifdef dbg_same_side DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); dmsg.banner(); dmsg << "checking same side"; diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index bd96c9be0..aac1ff8e9 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -969,6 +969,12 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array outNames(nOut); @@ -1291,11 +1297,26 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array(com_mod.uris[iUris].sdf(Ac)); + // d[iM].x(is,a) = static_cast(com_mod.uris[iUris].sdf_scaffold(Ac)); } + + if (com_mod.uris[iUris].scaffold_flag) { + cOut = cOut + 1; + int is = outS[cOut]; + int ie = is; + outS[cOut+1] = ie + 1; + outNames[cOut] = "URIS_SDF_SCAF_" + com_mod.uris[iUris].name; + for (int a = 0; a < msh.nNo; a++) { + int Ac = msh.gN(a); + d[iM].x(is,a) = static_cast(com_mod.uris[iUris].sdf_scaffold(Ac)); + } + } + } } From 4f93fa05d64cb0315c4e0fea75f9d77339aa9900 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Sun, 26 Oct 2025 21:07:11 -0700 Subject: [PATCH 11/19] Adding quadratic transition between RIS valve open and close thickness --- Code/Source/solver/fluid.cpp | 36 +++++++++++++++++++++++++------ Code/Source/solver/fsi.cpp | 42 +++++++++++++++++++++++++----------- 2 files changed, 60 insertions(+), 18 deletions(-) diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index 646c8d0c0..e71672e05 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -669,9 +669,32 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag double DDirTmp = 0.0; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { if (com_mod.uris[iUris].clsFlg) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + // Gradual transition using quadratic interpolation + if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { + // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps + double progress = static_cast(com_mod.uris[iUris].cnt) / + static_cast(com_mod.uris[iUris].DxClose.nrows()); + // Quadratic interpolation: smooth transition + double quad_progress = progress * progress; + sdf_deps_temp = com_mod.uris[iUris].sdf_deps + + quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); + } else { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + // // Gradual transition using quadratic interpolation + // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxOpen.nrows()) { + // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxOpen.nrows() steps + // double progress = static_cast(com_mod.uris[iUris].DxOpen.nrows() - com_mod.uris[iUris].cnt) / + // static_cast(com_mod.uris[iUris].DxOpen.nrows()); + // // Quadratic interpolation: smooth transition + // double quad_progress = progress * progress; + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + + // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); + // } else { + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + // } } if (distSrf(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ @@ -681,11 +704,12 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag if (com_mod.uris[iUris].scaffold_flag) { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; - } - if (distSrf_scaffold(iUris) <= sdf_deps_temp) { - DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ - (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} + + if (distSrf_scaffold(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + if (DDirTmp > DDir) {DDir = DDirTmp;} + } } } if (!com_mod.urisActFlag) {DDir = 0.0;} diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 990eb6bcc..ce149a842 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -206,16 +206,33 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar double sdf_deps_temp = 0; double DDirTmp = 0.0; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - // if (distSrf(iUris) <= com_mod.uris[iUris].sdf_deps) { - // DDirTmp = (1 + cos(pi*distSrf(iUris)/com_mod.uris[iUris].sdf_deps))/ - // (2*com_mod.uris[iUris].sdf_deps*com_mod.uris[iUris].sdf_deps); - // if (DDirTmp > DDir) {DDir = DDirTmp;} - // } - if (com_mod.uris[iUris].clsFlg) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + // Gradual transition using quadratic interpolation + if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { + // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps + double progress = static_cast(com_mod.uris[iUris].cnt) / + static_cast(com_mod.uris[iUris].DxClose.nrows()); + // Quadratic interpolation: smooth transition + double quad_progress = progress * progress; + sdf_deps_temp = com_mod.uris[iUris].sdf_deps + + quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); + } else { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + // // Gradual transition using quadratic interpolation + // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxOpen.nrows()) { + // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxOpen.nrows() steps + // double progress = static_cast(com_mod.uris[iUris].DxOpen.nrows() - com_mod.uris[iUris].cnt) / + // static_cast(com_mod.uris[iUris].DxOpen.nrows()); + // // Quadratic interpolation: smooth transition + // double quad_progress = progress * progress; + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + + // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); + // } else { + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + // } } if (distSrf(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ @@ -225,11 +242,12 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar if (com_mod.uris[iUris].scaffold_flag) { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; - } - if (distSrf_scaffold(iUris) <= sdf_deps_temp) { - DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ - (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} + + if (distSrf_scaffold(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + if (DDirTmp > DDir) {DDir = DDirTmp;} + } } } From 37f8e5800f2ccb16ece02e5ff4c50d19388805a7 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Thu, 20 Nov 2025 17:14:38 -0800 Subject: [PATCH 12/19] Adding valve velocity in the RIS penalty term. --- Code/Source/solver/ComMod.h | 28 +++-- Code/Source/solver/Parameters.cpp | 184 --------------------------- Code/Source/solver/Parameters.h | 3 - Code/Source/solver/distribute.cpp | 10 +- Code/Source/solver/fluid.cpp | 93 ++++++++++---- Code/Source/solver/fluid.h | 4 +- Code/Source/solver/fsi.cpp | 41 +++--- Code/Source/solver/uris.cpp | 203 +++++++++++++++++++++++------- Code/Source/solver/vtk_xml.cpp | 31 +++-- 9 files changed, 301 insertions(+), 296 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index fa144e8e3..e38a7f4cd 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -135,9 +135,6 @@ class bcType // Strong/Weak application of Dirichlet BC int clsFlgRis = 0; - // Strong/Weak application of Dirichlet BC - int clsFlgRis = 0; - // Pre/Res/Flat/Para... boundary types // // This stores differnt BCs as bitwise values. @@ -171,9 +168,6 @@ class bcType // RIS0D: resistance double resistance = 0.0; - // RIS0D: resistance - double resistance = 0.0; - // Penalty parameters for weakly applied Dir BC Vector tauB{0.0, 0.0}; //double tauB[2]; @@ -1458,6 +1452,15 @@ class urisType // Position coordinates (2D array: rows x columns). Array x; + // Position coordinates of last time step. + Array x_prev; + + // Position coordinates of next time step. + Array x_next; + + // Time derivative of the position coordinates (2D array). + Array v; + // Displacement (new) (2D array). Array Yd; @@ -1467,8 +1470,11 @@ class urisType // Default distance value of the valve boundary (valve thickness). double sdf_deps = 0.25; - // Default distance value of the valve boundary when the valve is closed. - double sdf_deps_close = 0.25; + // Default half valve thickness when valve is closed. + double sdf_deps_close = 0.2; + + // Default half scaffold thickness. + double sdf_deps_scaffold = 0.1; // Displacements of the valve when it opens (3D array). Array3 DxOpen; @@ -1485,9 +1491,15 @@ class urisType // Iteration count. int cnt = 1000000; + // URIS: signed distance function of each node to the uris valves (1D array). + Vector dirac_delta_func; + // URIS: signed distance function of each node to the uris valves (1D array). Vector sdf; + // URIS: time derivative of the signed distance field + Array sdf_t; + // Mesh scale factor. double scF = 1.0; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 75d229d2f..756abab2d 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2832,190 +2832,6 @@ void URISFaceParameters::set_values(tinyxml2::XMLElement* face_elem) } } - -////////////////////////////////////////////////////////// -// RISProjectionParameters // -////////////////////////////////////////////////////////// - -/// @brief Define the XML element name for mesh parameters. -const std::string RISProjectionParameters::xml_element_name_ = "Add_RIS_projection"; - -RISProjectionParameters::RISProjectionParameters() -{ - // A parameter that must be defined. - bool required = true; - - name = Parameter("name", "", required); - - set_parameter("Project_from_face", "", required, project_from_face); - set_parameter("Resistance", 1.e6, !required, resistance); - set_parameter("Projection_tolerance", 0.0, !required, projection_tolerance); -} - -void RISProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) -{ - using namespace tinyxml2; - std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; - - // Get the 'type' from the element. - const char* sname; - auto result = xml_elem->QueryStringAttribute("name", &sname); - if (sname == nullptr) { - throw std::runtime_error("No TYPE given in the XML element."); - } - name.set(std::string(sname)); - - using std::placeholders::_1; - using std::placeholders::_2; - - std::function ftpr = - std::bind( &RISProjectionParameters::set_parameter_value, *this, _1, _2); - - xml_util_set_parameters(ftpr, xml_elem, error_msg); -} - - -////////////////////////////////////////////////////////// -// URIS Mesh Parameters // -////////////////////////////////////////////////////////// -// [HZ] implemente URIS parameters here - -// Process parameters for the 'Add_URIS_mesh' XML element used for defining URIS mesh elements. - -/// @brief Define the XML element name for mesh parameters. -const std::string URISMeshParameters::xml_element_name_ = "Add_URIS_mesh"; - -URISMeshParameters::URISMeshParameters() -{ - bool required = true; - - // Mesh name from Add_mesh element. - name = Parameter("name", "", required); - - // Parameters under Add_mesh element. - // - set_parameter("Mesh_scale_factor", 1.0, !required, mesh_scale_factor); - set_parameter("Thickness", 0.04, !required, thickness); - set_parameter("Closed_thickness", 0.25, !required, close_thickness); - set_parameter("Resistance", 1.0e5, !required, resistance); - set_parameter("Closed_resistance", 1.0e5, !required, resistance_close); - set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); - set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); -} - -void URISMeshParameters::print_parameters() -{ - std::cout << std::endl; - std::cout << "---------------" << std::endl; - std::cout << "URIS Mesh Parameters" << std::endl; - std::cout << "---------------" << std::endl; - std::cout << name.name() << ": " << name.value() << std::endl; - - auto params_name_value = get_parameter_list(); - for (auto& [ key, value ] : params_name_value) { - std::cout << key << ": " << value << std::endl; - } - - for (auto& face : URIS_face_parameters) { - face->print_parameters(); - } -} - -void URISMeshParameters::set_values(tinyxml2::XMLElement* mesh_elem) -{ - using namespace tinyxml2; - std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; - auto item = mesh_elem->FirstChildElement(); - - while (item != nullptr) { - auto name = std::string(item->Value()); - - // Add_face sub-element. - if (name == URISFaceParameters::xml_element_name_) { - auto URIS_face_params = new URISFaceParameters(); - URIS_face_params->set_values(item); - URIS_face_parameters.push_back(URIS_face_params); - } else if (item->GetText() != nullptr) { - auto value = item->GetText(); - try { - set_parameter_value(name, value); - } catch (const std::bad_function_call& exception) { - throw std::runtime_error(error_msg + name + "'."); - } - } else { - throw std::runtime_error(error_msg + name + "'."); - } - - item = item->NextSiblingElement(); - } -} - - -////////////////////////////////////////////////////////// -// URIS Face Parameters // -////////////////////////////////////////////////////////// - -/// @brief Process parameters for the 'Add_URIS_face' XML element. -/// -/// Define the XML element name for face parameters. -const std::string URISFaceParameters::xml_element_name_ = "Add_URIS_face"; - -URISFaceParameters::URISFaceParameters() -{ - set_xml_element_name(xml_element_name_); - - // A parameter that must be defined. - bool required = true; - - name = Parameter("name", "", required); - - set_parameter("Face_file_path", "", !required, face_file_path); - set_parameter("Open_motion_file_path", "", !required, open_motion_file_path); - set_parameter("Close_motion_file_path", "", !required, close_motion_file_path); - - // set_parameter("End_nodes_face_file_path", "", !required, end_nodes_face_file_path); - // set_parameter("Quadrature_modifier_TRI3", (2.0/3.0), !required, quadrature_modifier_TRI3); -} - -void URISFaceParameters::print_parameters() -{ - std::cout << std::endl; - std::cout << "---------------" << std::endl; - std::cout << "URIS Face Parameters" << std::endl; - std::cout << "---------------" << std::endl; - std::cout << name.name() << ": " << name.value() << std::endl; - std::cout << face_file_path.name() << ": " << face_file_path.value() << std::endl; -} - -void URISFaceParameters::set_values(tinyxml2::XMLElement* face_elem) -{ - using namespace tinyxml2; - - std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; - const char* face_name; - auto result = face_elem->QueryStringAttribute("name", &face_name); - name.set(std::string(face_name)); - auto item = face_elem->FirstChildElement(); - - while (item != nullptr) { - auto name = std::string(item->Value()); - auto value = item->GetText(); - - if (value == nullptr) { - throw std::runtime_error(error_msg + name + "'."); - } - - try { - set_parameter_value(name, value); - } catch (const std::bad_function_call& exception) { - throw std::runtime_error(error_msg + name + "'."); - } - - item = item->NextSiblingElement(); - } -} - - ////////////////////////////////////////////////////////// // LinearAlgebraParameters // ////////////////////////////////////////////////////////// diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 6d26437e3..ec1853c0c 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1643,9 +1643,6 @@ class Parameters { void set_RIS_projection_values(tinyxml2::XMLElement* root_element); void set_URIS_mesh_values(tinyxml2::XMLElement* root_element); - void set_RIS_projection_values(tinyxml2::XMLElement* root_element); - void set_URIS_mesh_values(tinyxml2::XMLElement* root_element); - // Objects representing each parameter section of XML file. ContactParameters contact_parameters; GeneralSimulationParameters general_simulation_parameters; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 393bc6022..39defb1c6 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -345,8 +345,7 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.urisActFlag); cm.bcast(cm_mod, &com_mod.urisRes); cm.bcast(cm_mod, &com_mod.urisResClose); -======= ->>>>>>> edd5c69 (Adding option for using customized resistance value whenthe RIS valve is closed.) + cm.bcast(cm_mod, &com_mod.usePrecomp); if (com_mod.rmsh.isReqd) { auto& rmsh = com_mod.rmsh; @@ -1171,6 +1170,7 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].sdf_default); cm.bcast(cm_mod, &uris[iUris].sdf_deps); cm.bcast(cm_mod, &uris[iUris].sdf_deps_close); + cm.bcast(cm_mod, &uris[iUris].sdf_deps_scaffold); cm.bcast(cm_mod, &uris[iUris].clsFlg); cm.bcast(cm_mod, &uris[iUris].cnt); cm.bcast(cm_mod, &uris[iUris].scF); @@ -1356,6 +1356,9 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { uris[iUris].DxOpen.resize(OpenN(iUris), com_mod.nsd, uris[iUris].tnNo); uris[iUris].DxClose.resize(CloseN(iUris), com_mod.nsd, uris[iUris].tnNo); uris[iUris].x.resize(com_mod.nsd, uris[iUris].tnNo); + uris[iUris].x_prev.resize(com_mod.nsd, uris[iUris].tnNo); + uris[iUris].x_next.resize(com_mod.nsd, uris[iUris].tnNo); + uris[iUris].v.resize(com_mod.nsd, uris[iUris].tnNo); uris[iUris].Yd.resize(com_mod.nsd, uris[iUris].tnNo); DxOpenFlat[iUris].resize(DxOpenFlatSize(iUris)); DxCloseFlat[iUris].resize(DxCloseFlatSize(iUris)); @@ -1364,6 +1367,9 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { for (int iUris = 0; iUris < com_mod.nUris; iUris++) { cm.bcast(cm_mod, uris[iUris].x); + cm.bcast(cm_mod, uris[iUris].x_prev); + cm.bcast(cm_mod, uris[iUris].x_next); + cm.bcast(cm_mod, uris[iUris].v); cm.bcast(cm_mod, uris[iUris].Yd); cm.bcast(cm_mod, DxOpenFlat[iUris]); cm.bcast(cm_mod, DxCloseFlat[iUris]); diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index e71672e05..a68be0cdb 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -298,7 +298,6 @@ void bw_fluid_3d(ComMod& com_mod, const int eNoNw, const int eNoNq, const double u(2) = u(2) + Nw(a)*yl(2,a); ux(0,0) = ux(0,0) + Nwx(0,a)*yl(0,a); - // ux(1,0) = ux(1,1) + Nwx(1,a)*yl(0,a); ux(1,0) = ux(1,0) + Nwx(1,a)*yl(0,a); ux(2,0) = ux(2,0) + Nwx(2,a)*yl(0,a); ux(0,1) = ux(0,1) + Nwx(0,a)*yl(1,a); @@ -535,6 +534,15 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag // local tangent matrix (for a single element) Array3 lK(dof*dof,eNoN,eNoN); + // local velocity vector (for a single element) + Array vValve; + if (com_mod.urisFlag) { + vValve.resize(com_mod.nUris, nsd); + } else { + vValve.resize(0, 0); + } + vValve = 0.0; + double DDir = 0.0; // Loop over all elements of mesh @@ -653,11 +661,15 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag Vector distSrf(com_mod.nUris); Vector distSrf_scaffold(com_mod.nUris); distSrf = 0.0; + vValve = 0.0; distSrf_scaffold = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < com_mod.nUris; iUris++) { distSrf(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + for (int i = 0; i < nsd; i++) { + vValve(iUris,i) += fs[0].N(a,g) * com_mod.uris[iUris].sdf_t(i, Ac); + } if (com_mod.uris[iUris].scaffold_flag) { distSrf_scaffold(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf_scaffold(Ac)); } @@ -669,18 +681,19 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag double DDirTmp = 0.0; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { if (com_mod.uris[iUris].clsFlg) { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; // Gradual transition using quadratic interpolation - if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { - // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps - double progress = static_cast(com_mod.uris[iUris].cnt) / - static_cast(com_mod.uris[iUris].DxClose.nrows()); - // Quadratic interpolation: smooth transition - double quad_progress = progress * progress; - sdf_deps_temp = com_mod.uris[iUris].sdf_deps + - quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); - } else { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; - } + // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { + // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps + // double progress = static_cast(com_mod.uris[iUris].cnt) / + // static_cast(com_mod.uris[iUris].DxClose.nrows()); + // // Quadratic interpolation: smooth transition + // double quad_progress = progress * progress; + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + + // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); + // } else { + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + // } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; // // Gradual transition using quadratic interpolation @@ -703,7 +716,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag } if (com_mod.uris[iUris].scaffold_flag) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_scaffold; if (distSrf_scaffold(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ @@ -721,7 +734,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, - Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir); + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir, vValve); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); @@ -771,7 +784,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag if (nsd == 3) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); - fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir); + fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir, vValve); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); @@ -1478,7 +1491,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve) { #define n_debug_fluid3d_c #ifdef debug_fluid3d_c @@ -1509,7 +1522,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; - if (com_mod.uris[0].clsFlg) {Res = 1.0e5;} + if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} } double rho = dmn.prop[PhysicalProperyType::fluid_density]; @@ -1735,6 +1748,14 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2] + (Res*DDir)*u[2]); + if (com_mod.urisFlag) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + up[0] += -tauM * (Res * DDir)*(-vValve(iUris, 0)); + up[1] += -tauM * (Res * DDir)*(-vValve(iUris, 1)); + up[2] += -tauM * (Res * DDir)*(-vValve(iUris, 2)); + } + } + for (int a = 0; a < eNoNw; a++) { double uNx = u[0]*Nwx(0,a) + u[1]*Nwx(1,a) + u[2]*Nwx(2,a); // T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); @@ -1810,7 +1831,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve) { #define n_debug_fluid_3d_m #ifdef debug_fluid_3d_m @@ -1839,7 +1860,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e Res = 0.0; } else { Res = com_mod.urisRes; - if (com_mod.uris[0].clsFlg) {Res = 1.0e5;} + if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} } double ctM = 1.0; @@ -2095,6 +2116,14 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2] + (Res*DDir)*u[2]); + if (com_mod.urisFlag) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + up[0] += -tauM * (Res * DDir)*(-vValve(iUris, 0)); + up[1] += -tauM * (Res * DDir)*(-vValve(iUris, 1)); + up[2] += -tauM * (Res * DDir)*(-vValve(iUris, 2)); + } + } + double tauC, tauB, pa; double eps = std::numeric_limits::epsilon(); double ua[3] = {}; @@ -2174,8 +2203,8 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - - mu*K_inverse_darcy_permeability*Nw(a) - - (Res*DDir)*Nw(a); + - mu*K_inverse_darcy_permeability*Nw(a); + - (Res*DDir)*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[1]*mu_g*esNx[0][a]; @@ -2274,13 +2303,23 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Residual contribution Birkman term // Local residue + for (int a = 0; a < eNoNw; a++) { - lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[0]+up[0]) - + Res*DDir*w*Nw(a)*(u[0]+up[0]); - lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[1]+up[1]) - + Res*DDir*w*Nw(a)*(u[1]+up[1]); - lR(2,a) = lR(2,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[2]+up[2]) - + Res*DDir*w*Nw(a)*(u[2]+up[2]); + lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[0]+up[0]) + + Res*DDir*w*Nw(a)*(u[0]); //+up[0]); + lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[1]+up[1]) + + Res*DDir*w*Nw(a)*(u[1]); //+up[1]); + lR(2,a) = lR(2,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[2]+up[2]) + + Res*DDir*w*Nw(a)*(u[2]); //+up[2]); + + if (com_mod.urisFlag) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + lR(0,a) += (Res * DDir)*w*Nw(a)*(-vValve(iUris, 0)); + lR(1,a) += (Res * DDir)*w*Nw(a)*(-vValve(iUris, 1)); + lR(2,a) += (Res * DDir)*w*Nw(a)*(-vValve(iUris, 2)); + } + } + } } diff --git a/Code/Source/solver/fluid.h b/Code/Source/solver/fluid.h index 077e14373..198a23c92 100644 --- a/Code/Source/solver/fluid.h +++ b/Code/Source/solver/fluid.h @@ -36,12 +36,12 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir=0.0); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve); void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir=0.0); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve); void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, double& mu, double& mu_s, double& mu_x); diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index ce149a842..ec5d9c778 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -79,6 +79,13 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar double struct_3d_time = 0.0; double fluid_3d_time = 0.0; double DDir = 0.0; + Array vValve; + if (com_mod.urisFlag) { + vValve.resize(com_mod.nUris, nsd); + } else { + vValve.resize(0, 0); + } + vValve = 0.0; for (int e = 0; e < lM.nEl; e++) { // setting globals @@ -196,6 +203,9 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < com_mod.nUris; iUris++) { distSrf(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + for (int i = 0; i < nsd; i++) { + vValve(iUris,i) += fs_1[0].N(a,g) * com_mod.uris[iUris].sdf_t(i, Ac); + } if (com_mod.uris[iUris].scaffold_flag) { distSrf_scaffold(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf_scaffold(Ac)); } @@ -207,18 +217,19 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar double DDirTmp = 0.0; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { if (com_mod.uris[iUris].clsFlg) { - // Gradual transition using quadratic interpolation - if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { - // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps - double progress = static_cast(com_mod.uris[iUris].cnt) / - static_cast(com_mod.uris[iUris].DxClose.nrows()); - // Quadratic interpolation: smooth transition - double quad_progress = progress * progress; - sdf_deps_temp = com_mod.uris[iUris].sdf_deps + - quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); - } else { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; - } + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + // // Gradual transition using quadratic interpolation + // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { + // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps + // double progress = static_cast(com_mod.uris[iUris].cnt) / + // static_cast(com_mod.uris[iUris].DxClose.nrows()); + // // Quadratic interpolation: smooth transition + // double quad_progress = progress * progress; + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + + // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); + // } else { + // sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + // } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; // // Gradual transition using quadratic interpolation @@ -241,7 +252,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar } if (com_mod.uris[iUris].scaffold_flag) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_scaffold; if (distSrf_scaffold(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ @@ -263,7 +274,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman // fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); - fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir); + fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir, vValve); } break; @@ -341,7 +352,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman //fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); - fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir); + fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir, vValve); } break; case Equation_ustruct: diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 66b9e796e..64dacb715 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -54,8 +54,6 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { // Let's compute left side Array sUPS(1,com_mod.tnNo); - // std::cout << "com_mod.tnNo: " << com_mod.tnNo << std::endl; - // std::cout << "uris_obj.sdf size: " << uris_obj.sdf.size() << std::endl; sUPS = 0.0; for (size_t j = 0; j < sUPS.size(); j++) { @@ -481,7 +479,7 @@ void uris_read_msh(Simulation* simulation) { com_mod.urisResClose = param->resistance_close(); std::cout << "URIS resistance: " << com_mod.urisRes << std::endl; - std::cout << "URIS resistance when the valve is closed: " << com_mod.urisResClose << std::endl; + std::cout << "URIS resistance (closed): " << com_mod.urisResClose << std::endl; int nUris = simulation->parameters.URIS_mesh_parameters.size(); com_mod.nUris = nUris; @@ -502,11 +500,7 @@ void uris_read_msh(Simulation* simulation) { Array gX(0,0); std::string positive_flow_normal_file_path = param->positive_flow_normal_file_path(); - // [HZ] Need to read flow normal file (*.dat) into uris_obj.nrm - // lPtr => lPM%get(fTmp, "Positive flow normal file") - // fid = fTmp%open() - // READ (fid,*) uris(iUris)%nrm(:) - // CLOSE (fid) + std::ifstream file_stream; file_stream.open(positive_flow_normal_file_path); if (!file_stream.is_open()) { @@ -550,19 +544,13 @@ void uris_read_msh(Simulation* simulation) { } } scaffold_mesh.gIEN.clear(); - - std::cout << "Scaffold mesh is included for: " << uris_obj.name << std::endl; - std::cout << "Scaffold mesh nodes: " << uris_obj.scaffold_mesh.gnNo << std::endl; - std::cout << "Scaffold mesh elements: " << uris_obj.scaffold_mesh.gnEl << std::endl; } uris_obj.sdf_deps = param->thickness(); uris_obj.sdf_deps_close = param->close_thickness(); uris_obj.clsFlg = param->valve_starts_as_closed(); - // uris_obj.tnNo = 0; for (int iM = 0; iM < uris_obj.nFa; iM++) { - // Set as shell auto mesh_param = param->URIS_face_parameters[iM]; auto& mesh = uris_obj.msh[iM]; mesh.lShl = true; @@ -600,8 +588,6 @@ void uris_read_msh(Simulation* simulation) { } int dispNtOpen, dispNnOpen; file_stream >> dispNtOpen >> dispNnOpen; - // std::cout << "dispNtOpen: " << dispNtOpen << std::endl; - // std::cout << "dispNnOpen: " << dispNnOpen << std::endl; if (dispNnOpen != mesh.gnNo) { throw std::runtime_error("Mismatch in node numbers between URIS mesh and displacements."); @@ -612,7 +598,6 @@ void uris_read_msh(Simulation* simulation) { for (int a = 0; a < dispNnOpen; a++) { for (int i = 0; i < nsd; i++) { file_stream >> dispOpen(t,i,a); - // std::cout << "dispOpen: " << dispOpen(t,n,a) << std::endl; } } } @@ -627,8 +612,6 @@ void uris_read_msh(Simulation* simulation) { } int dispNtClose, dispNnClose; file_stream >> dispNtClose >> dispNnClose; - // std::cout << "dispNtClose: " << dispNtClose << std::endl; - // std::cout << "dispNnClose: " << dispNnClose << std::endl; if (dispNnClose != mesh.gnNo) { throw std::runtime_error("Mismatch in node numbers between URIS mesh and displacements."); @@ -639,7 +622,6 @@ void uris_read_msh(Simulation* simulation) { for (int a = 0; a < dispNnClose; a++) { for (int i = 0; i < nsd; i++) { file_stream >> dispClose(t,i,a); - // std::cout << "dispClose: " << dispClose(t,n,a) << std::endl; } } } @@ -647,9 +629,6 @@ void uris_read_msh(Simulation* simulation) { // To scale the mesh, while attaching x to gX int a = uris_obj.tnNo + mesh.gnNo; - // std::cout << "uris obj tnNo: " << uris_obj.tnNo << std::endl; - // std::cout << "mesh gnNo: " << mesh.gnNo << std::endl; - // std::cout << "mesh x size: " << mesh.x.nrows() << ", " << mesh.x.ncols() << std::endl; if (iM == 0) { gX.resize(nsd, a); @@ -714,7 +693,29 @@ void uris_read_msh(Simulation* simulation) { // dispClose.clear(); } uris_obj.x.resize(nsd, uris_obj.tnNo); - uris_obj.x = gX; + uris_obj.x_prev.resize(nsd, uris_obj.tnNo); + uris_obj.x_next.resize(nsd, uris_obj.tnNo); + if (!uris_obj.clsFlg) { + int nrowsDxOpen = uris_obj.DxOpen.nrows(); + for (int i = 0; i < uris_obj.x.nrows(); i++) { + for (int j = 0; j < uris_obj.x.ncols(); j++) { + uris_obj.x(i,j) = uris_obj.DxOpen(nrowsDxOpen-1,i,j); + uris_obj.x_prev(i,j) = uris_obj.DxOpen(nrowsDxOpen-1,i,j); + uris_obj.x_next(i,j) = uris_obj.DxOpen(nrowsDxOpen-1,i,j); + } + } + } else { + int nrowsDxClose = uris_obj.DxClose.nrows(); + for (int i = 0; i < uris_obj.x.nrows(); i++) { + for (int j = 0; j < uris_obj.x.ncols(); j++) { + uris_obj.x(i,j) = uris_obj.DxClose(nrowsDxClose-1,i,j); + uris_obj.x_prev(i,j) = uris_obj.DxClose(nrowsDxClose-1,i,j); + uris_obj.x_next(i,j) = uris_obj.DxClose(nrowsDxClose-1,i,j); + } + } + } + uris_obj.v.resize(nsd, uris_obj.tnNo); + uris_obj.v = 0.0; uris_obj.Yd.resize(nsd, uris_obj.tnNo); uris_obj.Yd = 0.0; // gX.clear(); @@ -786,8 +787,8 @@ void uris_write_vtus(ComMod& com_mod) { const int nsd = com_mod.nsd; const int nUris = com_mod.nUris; - // we plot coord + displacement - int nOut = 2; + // we plot coord + displacement + velocity + int nOut = 3; int outDof = nOut * nsd; std::vector outNames(nOut); Vector outS(nOut+1); @@ -802,11 +803,8 @@ void uris_write_vtus(ComMod& com_mod) { for (int iM = 0; iM < uris_obj.nFa; iM++) { auto& mesh = uris_obj.msh[iM]; int cOut = 0; - outS(cOut) = 0; // [HZ] Need to check this if it's 1 or 0 + outS(cOut) = 0; outS(cOut+1) = nsd; - // outNames[cOut] = ""; - - // outS = [0, 3] if (mesh.eType == ElementType::NRB) { throw std::runtime_error("Outputs for NURBS data is under development."); @@ -839,8 +837,6 @@ void uris_write_vtus(ComMod& com_mod) { int is = outS(cOut); // is = outS[1] = 3 int ie = is + l; // ie = 3 + 3 = 6 outS(cOut+1) = ie; // outS[2] = 7, outS = [0, 3, 6] - outNames[0] = "coordinates"; - outNames[1] = "URIS_displacement"; for (int a = 0; a < mesh.nNo; a++) { int Ac = mesh.gN(a); @@ -849,6 +845,24 @@ void uris_write_vtus(ComMod& com_mod) { } } + // For velocity + cOut += 1; // cOut = 2 + is = outS(cOut); // is = outS[2] = 6 + ie = is + l; // ie = 6 + 3 = 9 + outS(cOut+1) = ie; // outS[3] = 7, outS = [0, 3, 6, 9] + + for (int a = 0; a < mesh.nNo; a++) { + int Ac = mesh.gN(a); + for (int i = 0; i < nsd; i++) { + d[iM].x(is+i,a) = uris_obj.v(s+i,Ac); + } + } + + // [HZ] Need to check this if it's 1 or 0 + outNames[0] = "coordinates"; + outNames[1] = "URIS_displacement"; + outNames[2] = "URIS_velocity"; + nNo += mesh.nNo; nEl += mesh.nEl; } @@ -927,6 +941,11 @@ void uris_calc_sdf(ComMod& com_mod) { dmsg.banner(); #endif + using namespace consts; + + int cEq = com_mod.cEq; + auto& eq = com_mod.eq[cEq]; + auto& cm = com_mod.cm; auto& uris = com_mod.uris; const int nsd = com_mod.nsd; @@ -945,6 +964,17 @@ void uris_calc_sdf(ComMod& com_mod) { for (int i = 0; i < uris_obj.x.nrows(); i++) { for (int j = 0; j < uris_obj.x.ncols(); j++) { uris_obj.x(i,j) = uris_obj.DxOpen(cnt-1,i,j); + // Obtain x_n-1 and x_n+1 from the displacement history + if (cnt == uris_obj.DxOpen.nrows()) { + uris_obj.x_next(i,j) = uris_obj.DxOpen(cnt-1,i,j); + } else { + uris_obj.x_next(i,j) = uris_obj.DxOpen(cnt,i,j); + } + if (cnt == 1) { + uris_obj.x_prev(i,j) = uris_obj.DxOpen(cnt-1,i,j); + } else { + uris_obj.x_prev(i,j) = uris_obj.DxOpen(cnt-2,i,j); + } } } } else { @@ -952,12 +982,38 @@ void uris_calc_sdf(ComMod& com_mod) { for (int i = 0; i < uris_obj.x.nrows(); i++) { for (int j = 0; j < uris_obj.x.ncols(); j++) { uris_obj.x(i,j) = uris_obj.DxClose(cnt-1,i,j); + // Obtain x_n-1 and x_n+1 from the displacement history + if (cnt == uris_obj.DxClose.nrows()) { + uris_obj.x_next(i,j) = uris_obj.DxClose(cnt-1,i,j); + } else { + uris_obj.x_next(i,j) = uris_obj.DxClose(cnt,i,j); + } + if (cnt == 1) { + uris_obj.x_prev(i,j) = uris_obj.DxClose(cnt-1,i,j); + } else { + uris_obj.x_prev(i,j) = uris_obj.DxClose(cnt-2,i,j); + } } } } + + // uris(iUris)%v = (uris(iUris)%x - uris(iUris)%x_prev)/dt + // uris(iUris)%x_prev = uris(iUris)%x + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < uris_obj.tnNo; j++) { + uris_obj.v(i,j) = (uris_obj.x_next(i,j) - uris_obj.x_prev(i,j)) / (2*com_mod.dt); + // uris_obj.v(i,j) = (-eq.af * (1 - eq.af) * uris_obj.x_prev(i,j) + // + eq.af * (1 - eq.af) * uris_obj.x_next(i,j) + // + (1 - 2 * eq.af) * uris_obj.x(i,j)) / com_mod.dt; + // uris_obj.v(i,j) = (uris_obj.x(i,j) - uris_obj.x_prev(i,j)) / com_mod.dt; + // uris_obj.x_prev(i,j) = uris_obj.x(i,j); + } + } - // if (uris_obj.sdf.allocated() && cnt < uris_obj.cnt) {continue;} - if (uris_obj.sdf.size() > 0 && cnt < uris_obj.cnt) {continue;} + if (uris_obj.sdf.size() > 0 && cnt < uris_obj.cnt) { + uris_obj.sdf_t = 0.0; + continue; + } int max_eNoN = 0; for (int iM = 0; iM < uris_obj.nFa; iM++) { @@ -972,12 +1028,18 @@ void uris_calc_sdf(ComMod& com_mod) { if (uris_obj.sdf.size() <= 0) { uris_obj.sdf.resize(com_mod.tnNo); uris_obj.sdf = 0.0; + uris_obj.sdf_t.resize(nsd, com_mod.tnNo); + uris_obj.sdf_t = 0.0; + uris_obj.dirac_delta_func.resize(com_mod.tnNo); + uris_obj.dirac_delta_func = 0.0; } if (cm.idcm() == 0) { std::cout << "Recomputing SDF for " << uris_obj.name << std::endl; } uris_obj.sdf = uris_obj.sdf_default; + uris_obj.sdf_t = 0.0; + uris_obj.dirac_delta_func = 0.0; // Each time when the URIS moves (open/close), we need to // recompute the signed distance function. @@ -1010,10 +1072,14 @@ void uris_calc_sdf(ComMod& com_mod) { // leaflets will also be streched. Note that // this is a simplifying assumption. Vector xp(nsd); + for (int ca = 0; ca < com_mod.tnNo; ca++) { + // int mesh_ind = 0; + // for (int ca = 0; ca < com_mod.msh[mesh_ind].nNo; ca++) { double minS = std::numeric_limits::max(); for (int i = 0; i < nsd; i++) { xp(i) = com_mod.x(i,ca); + // xp(i) = com_mod.x(i,com_mod.msh[mesh_ind].gN(ca)); } // Is the node inside the BBox? bool inside = true; @@ -1051,16 +1117,15 @@ void uris_calc_sdf(ComMod& com_mod) { dS = std::sqrt(dS); if (dS < minS) { - minS = dS; - Ec = e; - jM = iM; + minS = dS; // closest distance + Ec = e; // closest URIS element index + jM = iM; // closest URIS mesh index } } } // We also need to compute the sign (above or below the valve). // Compute the element normal - auto& mesh = uris_obj.msh[jM]; xXi = 0.0; lX = 0.0; @@ -1094,29 +1159,74 @@ void uris_calc_sdf(ComMod& com_mod) { // } else { // dotP = 1.0; // } + // uris_obj.sdf[ca] = dotP * minS; // [HZ] Improved implementation for SDF sign + double sdf_sign = 1.0; if (uris_obj.clsFlg) { auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); if (dot_nrm < 0.0 && dotP < 0.0) { - dotP = -1.0; + sdf_sign = -1.0; } else { - dotP = 1.0; + sdf_sign = 1.0; } } else { - if (dotP < 0.0) { - dotP = -1.0; + auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); + // if (dotP < 0.0) { + if (dot_nrm < 0.0 && dotP < 0.0) { + sdf_sign = -1.0; } else { - dotP = 1.0; + sdf_sign = 1.0; } } - uris_obj.sdf[ca] = dotP * minS; + // uris_obj.sdf[com_mod.msh[mesh_ind].gN(ca)] = sdf_sign * minS; + uris_obj.sdf[ca] = sdf_sign * minS; + + Vector xp_plane(nsd), E1(nsd), E2(nsd), v(nsd), u_interp(nsd); + + // We need to interpolate valve velocity + // project xp onto the triangle plane + for (int i = 0; i < nsd; i++) { + xp_plane(i) = xp(i) - dotP * nV(i); + } + // compute barycentric (\xi, \eta) via the planar solve + for (int i = 0; i < nsd; i++) { + E1(i) = lX(i,1) - lX(i,0); + E2(i) = lX(i,2) - lX(i,0); + v(i) = xp_plane(i) - lX(i,0); + } + auto g11 = utils::norm(E1, E1); + auto g12 = utils::norm(E1, E2); + auto g22 = utils::norm(E2, E2); + auto b1 = utils::norm(v, E1); + auto b2 = utils::norm(v, E2); + double det = g11 * g22 - g12 * g12; + + double xi = (g22 * b1 - g12 * b2) / det; + double eta = (g11 * b2 - g12 * b1) / det; + + // shape functions: + double N1 = 1.0 - xi - eta; + double N2 = xi; + double N3 = eta; + + // interpolate the valve velocity: + for (int i = 0; i < nsd; i++) { + u_interp(i) = N1 * uris_obj.v(i, mesh.IEN(0,Ec)) + + N2 * uris_obj.v(i, mesh.IEN(1,Ec)) + + N3 * uris_obj.v(i, mesh.IEN(2,Ec)); + } + + // Store interpolated velocity + for (int i = 0; i < nsd; i++) { + uris_obj.sdf_t(i,ca) = u_interp(i); + } + } } } - // Compute the SDF for the scaffold mesh for (int iUris = 0; iUris < nUris; iUris++) { auto& uris_obj = uris[iUris]; @@ -1161,6 +1271,8 @@ void uris_calc_sdf(ComMod& com_mod) { // this is a simplifying assumption. Vector xp_scaf(nsd); for (int ca = 0; ca < com_mod.tnNo; ca++) { + // int mesh_ind = 0; + // for (int ca = 0; ca < com_mod.msh[mesh_ind].nNo; ca++) { // std::cout << "URIS scaffold SDF index: " << ca << std::endl; double minS_scaf = std::numeric_limits::max(); for (int i = 0; i < nsd; i++) { @@ -1202,6 +1314,7 @@ void uris_calc_sdf(ComMod& com_mod) { } } uris_obj.sdf_scaffold[ca] = minS_scaf; + // uris_obj.sdf_scaffold[com_mod.msh[mesh_ind].gN(ca)] = minS_scaf; } } } diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index aac1ff8e9..7911d1204 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -967,8 +967,8 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, const Array(com_mod.uris[iUris].sdf(Ac)); - // d[iM].x(is,a) = static_cast(com_mod.uris[iUris].sdf_scaffold(Ac)); } + } + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + cOut = cOut + 1; + int is = outS[cOut]; + int ie = is + nsd - 1; + outS[cOut+1] = ie + 1; + outNames[cOut] = "URIS_SDF_VEL_" + com_mod.uris[iUris].name; + for (int a = 0; a < msh.nNo; a++) { + int Ac = msh.gN(a); + for (int b = is; b <= ie; b++) { + d[iM].x(b,a) = static_cast(com_mod.uris[iUris].sdf_t(b-is, Ac)); + } + } + } + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { if (com_mod.uris[iUris].scaffold_flag) { cOut = cOut + 1; int is = outS[cOut]; @@ -1316,9 +1328,8 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array(com_mod.uris[iUris].sdf_scaffold(Ac)); } } - - } - } + } + } } // iM for loop From 5c92b35372c9dbb86d2d4d21885f645204e59fbb Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Tue, 13 Jan 2026 14:24:54 -0800 Subject: [PATCH 13/19] Adding options considering valve velocity and different resistances values for RIS valves --- Code/Source/solver/ComMod.h | 19 ++- Code/Source/solver/Parameters.cpp | 1 + Code/Source/solver/Parameters.h | 1 + Code/Source/solver/distribute.cpp | 5 +- Code/Source/solver/fluid.cpp | 124 ++++++++-------- Code/Source/solver/fluid.h | 4 +- Code/Source/solver/fsi.cpp | 63 +++++--- Code/Source/solver/initialize.cpp | 11 ++ Code/Source/solver/uris.cpp | 237 +++++++++++------------------- 9 files changed, 213 insertions(+), 252 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index e38a7f4cd..adb2daf30 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1467,8 +1467,8 @@ class urisType // Default signed distance value away from the valve. double sdf_default = 1000.0; - // Default distance value of the valve boundary (valve thickness). - double sdf_deps = 0.25; + // Default distance value of the valve boundary (half valve thickness). + double sdf_deps = 0.2; // Default half valve thickness when valve is closed. double sdf_deps_close = 0.2; @@ -1476,6 +1476,15 @@ class urisType // Default half scaffold thickness. double sdf_deps_scaffold = 0.1; + /// @brief URIS resistance + double resistance; + + /// @brief URIS resistance when the valve is closed + double resistance_close; + + // Whether to use the valve velocity. + bool use_valve_velocity = false; + // Displacements of the valve when it opens (3D array). Array3 DxOpen; @@ -1615,12 +1624,6 @@ class ComMod { /// @brief Number of URIS surfaces (uninitialized, to be set later) int nUris; - /// @brief URIS resistance - double urisRes; - - /// @brief URIS resistance when the valve is closed - double urisResClose; - /// @brief Whether to use precomputed state-variable solutions bool usePrecomp = false; //----- int members -----// diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 756abab2d..9ec827ec7 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2715,6 +2715,7 @@ URISMeshParameters::URISMeshParameters() set_parameter("Closed_thickness", 0.25, !required, close_thickness); set_parameter("Resistance", 1.0e5, !required, resistance); set_parameter("Closed_resistance", 1.0e5, !required, resistance_close); + set_parameter("Use_valve_velocity", false, !required, use_valve_velocity); set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); set_parameter("Scaffold_file_path", "", !required, scaffold_file_path); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index ec1853c0c..1de79c085 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1611,6 +1611,7 @@ class URISMeshParameters : public ParameterLists Parameter close_thickness; // Thickness of the valve when it is closed Parameter resistance; // Resistance of the valve Parameter resistance_close; // Resistance of the valve when it is closed + Parameter use_valve_velocity; // Whether to use the valve velocity Parameter valve_starts_as_closed; // Whether the valve starts as closed Parameter positive_flow_normal_file_path; // File path for the positive flow normal Parameter scaffold_file_path; // File path for the scaffold diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 39defb1c6..e3a941276 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -343,8 +343,6 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.ris0DFlag); cm.bcast(cm_mod, &com_mod.urisFlag); cm.bcast(cm_mod, &com_mod.urisActFlag); - cm.bcast(cm_mod, &com_mod.urisRes); - cm.bcast(cm_mod, &com_mod.urisResClose); cm.bcast(cm_mod, &com_mod.usePrecomp); if (com_mod.rmsh.isReqd) { @@ -1171,12 +1169,15 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].sdf_deps); cm.bcast(cm_mod, &uris[iUris].sdf_deps_close); cm.bcast(cm_mod, &uris[iUris].sdf_deps_scaffold); + cm.bcast(cm_mod, &uris[iUris].resistance); + cm.bcast(cm_mod, &uris[iUris].resistance_close); cm.bcast(cm_mod, &uris[iUris].clsFlg); cm.bcast(cm_mod, &uris[iUris].cnt); cm.bcast(cm_mod, &uris[iUris].scF); cm.bcast(cm_mod, uris[iUris].nrm); cm.bcast(cm_mod, &uris[iUris].scaffold_flag); + cm.bcast(cm_mod, &uris[iUris].use_valve_velocity); if (uris[iUris].scaffold_flag) { cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.lShl); cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.nEl); diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index a68be0cdb..a4f40127b 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -535,15 +535,13 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag Array3 lK(dof*dof,eNoN,eNoN); // local velocity vector (for a single element) - Array vValve; + Array vValve(0, 0); if (com_mod.urisFlag) { vValve.resize(com_mod.nUris, nsd); - } else { - vValve.resize(0, 0); } vValve = 0.0; - double DDir = 0.0; + double ris_factor = 0.0; // Loop over all elements of mesh // @@ -661,14 +659,16 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag Vector distSrf(com_mod.nUris); Vector distSrf_scaffold(com_mod.nUris); distSrf = 0.0; - vValve = 0.0; distSrf_scaffold = 0.0; + vValve = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < com_mod.nUris; iUris++) { distSrf(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); - for (int i = 0; i < nsd; i++) { - vValve(iUris,i) += fs[0].N(a,g) * com_mod.uris[iUris].sdf_t(i, Ac); + if (com_mod.uris[iUris].use_valve_velocity) { + for (int i = 0; i < nsd; i++) { + vValve(iUris,i) += fs[0].N(a,g) * com_mod.uris[iUris].sdf_t(i, Ac); + } } if (com_mod.uris[iUris].scaffold_flag) { distSrf_scaffold(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf_scaffold(Ac)); @@ -676,12 +676,24 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag } } - DDir = 0.0; + double DDir = 0.0; + double ris_resistance = 0.0; double sdf_deps_temp = 0; double DDirTmp = 0.0; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (com_mod.uris[iUris].scaffold_flag) { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_scaffold; + ris_resistance = com_mod.uris[iUris].resistance_close; + if (distSrf_scaffold(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + // if (DDirTmp > DDir) {DDir = DDirTmp;} + } + } + if (com_mod.uris[iUris].clsFlg) { sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + ris_resistance = com_mod.uris[iUris].resistance_close; // Gradual transition using quadratic interpolation // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps @@ -696,6 +708,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag // } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + ris_resistance = com_mod.uris[iUris].resistance; // // Gradual transition using quadratic interpolation // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxOpen.nrows()) { // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxOpen.nrows() steps @@ -712,20 +725,15 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag if (distSrf(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} - } - - if (com_mod.uris[iUris].scaffold_flag) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_scaffold; - - if (distSrf_scaffold(iUris) <= sdf_deps_temp) { - DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ - (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} - } + // if (DDirTmp > DDir) {DDir = DDirTmp;} } + DDir = DDirTmp; + } + if (com_mod.urisActFlag) { + ris_factor = ris_resistance * DDir; + } else { + ris_factor = 0.0; } - if (!com_mod.urisActFlag) {DDir = 0.0;} } // Compute momentum residual and tangent matrix. @@ -734,7 +742,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, - Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir, vValve); + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, ris_factor, vValve); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); @@ -784,7 +792,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag if (nsd == 3) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); - fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir, vValve); + fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, ris_factor, vValve); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); @@ -1491,7 +1499,8 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, + double ris_factor, const Array& vValve) { #define n_debug_fluid3d_c #ifdef debug_fluid3d_c @@ -1517,14 +1526,6 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e const double ctM = 1.0; const double ctC = 36.0; - double Res; - if (!com_mod.urisFlag) { - Res = 0.0; - } else { - Res = com_mod.urisRes; - if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} - } - double rho = dmn.prop[PhysicalProperyType::fluid_density]; double f[3]; f[0] = dmn.prop[PhysicalProperyType::f_x]; @@ -1714,7 +1715,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // In case of unfitted RIS, compute the delta function at the quad point, // add the additional value to the stabilization param - kT = kT + pow(Res*DDir, 2.0); + kT = kT + pow(ris_factor, 2.0); double kU = u[0]*u[0]*Kxi(0,0) + u[1]*u[0]*Kxi(1,0) + u[2]*u[0]*Kxi(2,0) + u[0]*u[1]*Kxi(0,1) + u[1]*u[1]*Kxi(1,1) + u[2]*u[1]*Kxi(2,1) @@ -1742,17 +1743,17 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2]); up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability*u[0] - + (Res*DDir)*u[0]); + + ris_factor*u[0]); up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability*u[1] - + (Res*DDir)*u[1]); + + ris_factor*u[1]); up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2] - + (Res*DDir)*u[2]); + + ris_factor*u[2]); if (com_mod.urisFlag) { for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - up[0] += -tauM * (Res * DDir)*(-vValve(iUris, 0)); - up[1] += -tauM * (Res * DDir)*(-vValve(iUris, 1)); - up[2] += -tauM * (Res * DDir)*(-vValve(iUris, 2)); + up[0] += -tauM * ris_factor*(-vValve(iUris, 0)); + up[1] += -tauM * ris_factor*(-vValve(iUris, 1)); + up[2] += -tauM * ris_factor*(-vValve(iUris, 2)); } } @@ -1763,7 +1764,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a) - - (Res*DDir)*Nw(a); + - ris_factor*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[1]*mu_g*esNx[0][a]; @@ -1831,7 +1832,8 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, + double ris_factor, const Array& vValve) { #define n_debug_fluid_3d_m #ifdef debug_fluid_3d_m @@ -1855,14 +1857,6 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double Res; - if (!com_mod.urisFlag) { - Res = 0.0; - } else { - Res = com_mod.urisRes; - if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} - } - double ctM = 1.0; double ctC = 36.0; @@ -2075,7 +2069,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // In case of unfitted RIS, compute the delta function at the quad point, // add the additional value to the stabilization param - kT = kT + pow(Res*DDir, 2.0); + kT = kT + pow(ris_factor, 2.0); double kU = u[0]*u[0]*Kxi(0,0) + u[1]*u[0]*Kxi(1,0) + u[2]*u[0]*Kxi(2,0) + u[0]*u[1]*Kxi(0,1) + u[1]*u[1]*Kxi(1,1) + u[2]*u[1]*Kxi(2,1) @@ -2110,17 +2104,17 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2]); up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability * u[0] - + (Res*DDir)*u[0]); + + ris_factor*u[0]); up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability * u[1] - + (Res*DDir)*u[1]); + + ris_factor*u[1]); up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2] - + (Res*DDir)*u[2]); + + ris_factor*u[2]); if (com_mod.urisFlag) { for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - up[0] += -tauM * (Res * DDir)*(-vValve(iUris, 0)); - up[1] += -tauM * (Res * DDir)*(-vValve(iUris, 1)); - up[2] += -tauM * (Res * DDir)*(-vValve(iUris, 2)); + up[0] += -tauM * ris_factor*(-vValve(iUris, 0)); + up[1] += -tauM * ris_factor*(-vValve(iUris, 1)); + up[2] += -tauM * ris_factor*(-vValve(iUris, 2)); } } @@ -2204,7 +2198,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); - - (Res*DDir)*Nw(a); + - ris_factor*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[1]*mu_g*esNx[0][a]; @@ -2241,7 +2235,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e lK(0,a,b) = lK(0,a,b) + wl*(T2 + T1); // lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) - + (Res*DDir)*wl*Nw(b)*Nw(a); + + ris_factor*wl*Nw(b)*Nw(a); // dRm_a1/du_b2 T2 = mu*rM[1][0] + tauC*rM[0][1] + esNx[0][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][0][b]; @@ -2260,7 +2254,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e lK(5,a,b) = lK(5,a,b) + wl*(T2 + T1); // lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) - + (Res*DDir)*wl*Nw(b)*Nw(a); + + ris_factor*wl*Nw(b)*Nw(a); // dRm_a2/du_b3 T2 = mu*rM[2][1] + tauC*rM[1][2] + esNx[1][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][1][b]; @@ -2279,7 +2273,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e lK(10,a,b) = lK(10,a,b) + wl*(T2 + T1); // lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) - + (Res*DDir)*wl*Nw(b)*Nw(a); + + ris_factor*wl*Nw(b)*Nw(a); //dmsg << "lK(10,a,b): " << lK(10,a,b); } } @@ -2306,17 +2300,17 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e for (int a = 0; a < eNoNw; a++) { lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[0]+up[0]) - + Res*DDir*w*Nw(a)*(u[0]); //+up[0]); + + ris_factor*w*Nw(a)*(u[0]); //+up[0]); lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[1]+up[1]) - + Res*DDir*w*Nw(a)*(u[1]); //+up[1]); + + ris_factor*w*Nw(a)*(u[1]); //+up[1]); lR(2,a) = lR(2,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[2]+up[2]) - + Res*DDir*w*Nw(a)*(u[2]); //+up[2]); + + ris_factor*w*Nw(a)*(u[2]); //+up[2]); if (com_mod.urisFlag) { for (int iUris = 0; iUris < com_mod.nUris; iUris++) { - lR(0,a) += (Res * DDir)*w*Nw(a)*(-vValve(iUris, 0)); - lR(1,a) += (Res * DDir)*w*Nw(a)*(-vValve(iUris, 1)); - lR(2,a) += (Res * DDir)*w*Nw(a)*(-vValve(iUris, 2)); + lR(0,a) += ris_factor*w*Nw(a)*(-vValve(iUris, 0)); + lR(1,a) += ris_factor*w*Nw(a)*(-vValve(iUris, 1)); + lR(2,a) += ris_factor*w*Nw(a)*(-vValve(iUris, 2)); } } diff --git a/Code/Source/solver/fluid.h b/Code/Source/solver/fluid.h index 198a23c92..48092af94 100644 --- a/Code/Source/solver/fluid.h +++ b/Code/Source/solver/fluid.h @@ -36,12 +36,12 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double ris_factor, const Array& vValve); void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir, const Array& vValve); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double ris_factor, const Array& vValve); void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, double& mu, double& mu_s, double& mu_x); diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index ec5d9c778..61cf3ea8b 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -23,7 +23,7 @@ namespace fsi { void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Yg, const Array& Dg) { - #define n_debug_construct_fsi + #define n_debug_construct_fsi #ifdef debug_construct_fsi DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); @@ -62,6 +62,15 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar dmsg << "pS0.size(): " << pS0.size(); #endif + // if (com_mod.cm.idcm() == 0) { + // std::cout << "tnNo: " << com_mod.tnNo << std::endl; + // std::cout << "tDof: " << com_mod.tDof << std::endl; + // std::cout << "dof: " << com_mod.dof << std::endl; + // std::cout << "mesh name: " << lM.name << std::endl; + // std::cout << "mesh gnno: " << lM.gnNo << std::endl; + // std::cout << "mesh nEl: " << lM.nEl << std::endl; + // } + Vector ptr(eNoN); Array3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN); Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), @@ -78,7 +87,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // double struct_3d_time = 0.0; double fluid_3d_time = 0.0; - double DDir = 0.0; + Array vValve; if (com_mod.urisFlag) { vValve.resize(com_mod.nUris, nsd); @@ -86,6 +95,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar vValve.resize(0, 0); } vValve = 0.0; + double ris_factor = 0.0; for (int e = 0; e < lM.nEl; e++) { // setting globals @@ -199,12 +209,15 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar Vector distSrf_scaffold(com_mod.nUris); distSrf = 0.0; distSrf_scaffold = 0.0; + vValve = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < com_mod.nUris; iUris++) { distSrf(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); - for (int i = 0; i < nsd; i++) { - vValve(iUris,i) += fs_1[0].N(a,g) * com_mod.uris[iUris].sdf_t(i, Ac); + if (com_mod.uris[iUris].use_valve_velocity) { + for (int i = 0; i < nsd; i++) { + vValve(iUris,i) += fs_1[0].N(a,g) * com_mod.uris[iUris].sdf_t(i, Ac); + } } if (com_mod.uris[iUris].scaffold_flag) { distSrf_scaffold(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf_scaffold(Ac)); @@ -212,13 +225,25 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar } } - DDir = 0.0; + double DDir = 0.0; + double ris_resistance = 0.0; double sdf_deps_temp = 0; double DDirTmp = 0.0; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (com_mod.uris[iUris].scaffold_flag) { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_scaffold; + ris_resistance = com_mod.uris[iUris].resistance_close; + if (distSrf_scaffold(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + // if (DDirTmp > DDir) {DDir = DDirTmp;} + } + } + if (com_mod.uris[iUris].clsFlg) { sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; - // // Gradual transition using quadratic interpolation + ris_resistance = com_mod.uris[iUris].resistance_close; + // Gradual transition using quadratic interpolation // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps // double progress = static_cast(com_mod.uris[iUris].cnt) / @@ -232,6 +257,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + ris_resistance = com_mod.uris[iUris].resistance; // // Gradual transition using quadratic interpolation // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxOpen.nrows()) { // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxOpen.nrows() steps @@ -248,22 +274,15 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar if (distSrf(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} + // if (DDirTmp > DDir) {DDir = DDirTmp;} } - - if (com_mod.uris[iUris].scaffold_flag) { - sdf_deps_temp = com_mod.uris[iUris].sdf_deps_scaffold; - - if (distSrf_scaffold(iUris) <= sdf_deps_temp) { - DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ - (2*sdf_deps_temp*sdf_deps_temp); - if (DDirTmp > DDir) {DDir = DDirTmp;} - } - } - + DDir = DDirTmp; + } + if (com_mod.urisActFlag) { + ris_factor = ris_resistance * DDir; + } else { + ris_factor = 0.0; } - - if (!com_mod.urisActFlag) {DDir = 0.0;} } if (nsd == 3) { @@ -274,7 +293,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman // fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); - fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir, vValve); + fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, ris_factor, vValve); } break; @@ -352,7 +371,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman //fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); - fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir, vValve); + fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, ris_factor, vValve); } break; case Equation_ustruct: diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7f1e1ae84..29e4b4909 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -847,6 +847,17 @@ void initialize(Simulation* simulation, Vector& timeP) std::fill(com_mod.rmsh.flag.begin(), com_mod.rmsh.flag.end(), false); com_mod.resetSim = false; + + if (com_mod.urisFlag) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + auto& uris_obj = com_mod.uris[iUris]; + uris_obj.sdf.resize(com_mod.tnNo); + uris_obj.sdf_t.resize(nsd, com_mod.tnNo); + if (uris_obj.scaffold_flag) { + uris_obj.sdf_scaffold.resize(com_mod.tnNo); + } + } + } } //----------- diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 64dacb715..3f3fcb50e 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -15,6 +15,7 @@ #include "read_msh.h" #include "VtkData.h" + namespace uris { /// @brief This subroutine computes the mean pressure and flux on the @@ -330,7 +331,6 @@ void uris_update_disp(ComMod& com_mod, CmMod& cm_mod) { uris_obj.Yd(i,nd) /= divisor; } } - } } @@ -474,13 +474,6 @@ void uris_read_msh(Simulation* simulation) { com_mod.urisFlag = true; com_mod.urisActFlag = true; - auto param = simulation->parameters.URIS_mesh_parameters[0]; - com_mod.urisRes = param->resistance(); - com_mod.urisResClose = param->resistance_close(); - - std::cout << "URIS resistance: " << com_mod.urisRes << std::endl; - std::cout << "URIS resistance (closed): " << com_mod.urisResClose << std::endl; - int nUris = simulation->parameters.URIS_mesh_parameters.size(); com_mod.nUris = nUris; @@ -497,6 +490,7 @@ void uris_read_msh(Simulation* simulation) { uris_obj.nFa = param->URIS_face_parameters.size(); uris_obj.msh.resize(uris_obj.nFa); uris_obj.nrm.resize(nsd); + Array gX(0,0); std::string positive_flow_normal_file_path = param->positive_flow_normal_file_path(); @@ -548,7 +542,10 @@ void uris_read_msh(Simulation* simulation) { uris_obj.sdf_deps = param->thickness(); uris_obj.sdf_deps_close = param->close_thickness(); + uris_obj.resistance = param->resistance(); + uris_obj.resistance_close = param->resistance_close(); uris_obj.clsFlg = param->valve_starts_as_closed(); + uris_obj.use_valve_velocity = param->use_valve_velocity(); for (int iM = 0; iM < uris_obj.nFa; iM++) { auto mesh_param = param->URIS_face_parameters[iM]; @@ -752,10 +749,6 @@ void uris_read_msh(Simulation* simulation) { } } mesh.gIEN.clear(); - - std::cout << "Number of uris elements nel: " << mesh.nEl << std::endl; - std::cout << "URIS mesh IEN size: " << mesh.IEN.nrows() << ", " << mesh.IEN.ncols() << std::endl; - } if (uris_obj.nFa > 0) { @@ -945,74 +938,54 @@ void uris_calc_sdf(ComMod& com_mod) { int cEq = com_mod.cEq; auto& eq = com_mod.eq[cEq]; - auto& cm = com_mod.cm; auto& uris = com_mod.uris; const int nsd = com_mod.nsd; const int nUris = com_mod.nUris; Array xXi(nsd, nsd-1); + Vector minb(nsd), maxb(nsd), extra(nsd), xp(nsd); + Vector minb_scaf(nsd), maxb_scaf(nsd), extra_scaf(nsd), xp_scaf(nsd); + Vector xp_plane(nsd), E1(nsd), E2(nsd), v(nsd), u_interp(nsd); + double extra_val = 0.1, extra_scaf_val = 0.1; // BBox size adjustment for (int iUris = 0; iUris < nUris; iUris++) { // We need to check if the valve needs to move auto& uris_obj = uris[iUris]; - // Compute the SDF for the uris valves - int cnt = 0; - if (!uris_obj.clsFlg) { - cnt = std::min(uris_obj.cnt, uris_obj.DxOpen.nrows()); - for (int i = 0; i < uris_obj.x.nrows(); i++) { - for (int j = 0; j < uris_obj.x.ncols(); j++) { - uris_obj.x(i,j) = uris_obj.DxOpen(cnt-1,i,j); - // Obtain x_n-1 and x_n+1 from the displacement history - if (cnt == uris_obj.DxOpen.nrows()) { - uris_obj.x_next(i,j) = uris_obj.DxOpen(cnt-1,i,j); - } else { - uris_obj.x_next(i,j) = uris_obj.DxOpen(cnt,i,j); - } - if (cnt == 1) { - uris_obj.x_prev(i,j) = uris_obj.DxOpen(cnt-1,i,j); - } else { - uris_obj.x_prev(i,j) = uris_obj.DxOpen(cnt-2,i,j); - } - } - } - } else { - cnt = std::min(uris_obj.cnt, uris_obj.DxClose.nrows()); - for (int i = 0; i < uris_obj.x.nrows(); i++) { - for (int j = 0; j < uris_obj.x.ncols(); j++) { - uris_obj.x(i,j) = uris_obj.DxClose(cnt-1,i,j); - // Obtain x_n-1 and x_n+1 from the displacement history - if (cnt == uris_obj.DxClose.nrows()) { - uris_obj.x_next(i,j) = uris_obj.DxClose(cnt-1,i,j); - } else { - uris_obj.x_next(i,j) = uris_obj.DxClose(cnt,i,j); - } - if (cnt == 1) { - uris_obj.x_prev(i,j) = uris_obj.DxClose(cnt-1,i,j); - } else { - uris_obj.x_prev(i,j) = uris_obj.DxClose(cnt-2,i,j); - } - } + int cnt = uris_obj.clsFlg ? std::min(uris_obj.cnt, uris_obj.DxClose.nrows()) : + std::min(uris_obj.cnt, uris_obj.DxOpen.nrows()); + Array3& Dx = uris_obj.clsFlg ? uris_obj.DxClose : uris_obj.DxOpen; + + for (int i = 0; i < uris_obj.x.nrows(); i++) { + for (int j = 0; j < uris_obj.x.ncols(); j++) { + uris_obj.x(i,j) = Dx(cnt-1,i,j); } } - // uris(iUris)%v = (uris(iUris)%x - uris(iUris)%x_prev)/dt - // uris(iUris)%x_prev = uris(iUris)%x for (int i = 0; i < nsd; i++) { for (int j = 0; j < uris_obj.tnNo; j++) { - uris_obj.v(i,j) = (uris_obj.x_next(i,j) - uris_obj.x_prev(i,j)) / (2*com_mod.dt); - // uris_obj.v(i,j) = (-eq.af * (1 - eq.af) * uris_obj.x_prev(i,j) - // + eq.af * (1 - eq.af) * uris_obj.x_next(i,j) - // + (1 - 2 * eq.af) * uris_obj.x(i,j)) / com_mod.dt; - // uris_obj.v(i,j) = (uris_obj.x(i,j) - uris_obj.x_prev(i,j)) / com_mod.dt; - // uris_obj.x_prev(i,j) = uris_obj.x(i,j); + if (uris_obj.use_valve_velocity) { + // uris_obj.v(i,j) = (uris_obj.x_next(i,j) - uris_obj.x_prev(i,j)) / (2*com_mod.dt); + + // uris_obj.v(i,j) = (-eq.af * (1 - eq.af) * uris_obj.x_prev(i,j) + // + eq.af * (1 - eq.af) * uris_obj.x_next(i,j) + // + (1 - 2 * eq.af) * uris_obj.x(i,j)) / com_mod.dt; + uris_obj.v(i,j) = (uris_obj.x(i,j) - uris_obj.x_prev(i,j)) / com_mod.dt; + } else { + uris_obj.v(i,j) = 0.0; + } + uris_obj.x_prev(i,j) = uris_obj.x(i,j); } } - if (uris_obj.sdf.size() > 0 && cnt < uris_obj.cnt) { + if (com_mod.cTS > 1 && cnt < uris_obj.cnt) { + // std::cout << "cts:" << com_mod.cTS << std::endl; uris_obj.sdf_t = 0.0; continue; + } else { + uris_obj.sdf = uris_obj.sdf_default; + uris_obj.sdf_t = 0.0; } int max_eNoN = 0; @@ -1022,24 +995,12 @@ void uris_calc_sdf(ComMod& com_mod) { max_eNoN = mesh.eNoN; } } - Array lX(nsd, max_eNoN); - // if (!uris_obj.sdf.allocated()) { - if (uris_obj.sdf.size() <= 0) { - uris_obj.sdf.resize(com_mod.tnNo); - uris_obj.sdf = 0.0; - uris_obj.sdf_t.resize(nsd, com_mod.tnNo); - uris_obj.sdf_t = 0.0; - uris_obj.dirac_delta_func.resize(com_mod.tnNo); - uris_obj.dirac_delta_func = 0.0; - } + Array lX(nsd, max_eNoN); if (cm.idcm() == 0) { std::cout << "Recomputing SDF for " << uris_obj.name << std::endl; } - uris_obj.sdf = uris_obj.sdf_default; - uris_obj.sdf_t = 0.0; - uris_obj.dirac_delta_func = 0.0; // Each time when the URIS moves (open/close), we need to // recompute the signed distance function. @@ -1054,14 +1015,11 @@ void uris_calc_sdf(ComMod& com_mod) { } // For each coordinate dimension, find the minimum and maximum in uris_obj.x. - double extra_val = 0.1; // [HZ] The BBox is 10% larger than the actual valve, default is 0.1 for (int i = 0; i < nsd; i++) { for (int j = 0; j < uris_obj.x.ncols(); j++) { double val = uris_obj.x(i,j); - if (val < minb(i)) - minb(i) = val; - if (val > maxb(i)) - maxb(i) = val; + minb(i) = std::min(minb(i), val); + maxb(i) = std::max(maxb(i), val); } extra(i) = (maxb(i) - minb(i)) * extra_val; } @@ -1070,9 +1028,7 @@ void uris_calc_sdf(ComMod& com_mod) { // means that the valves will be morphed based on the fluid mesh // motion. If the fluid mesh stretches near the valve, the valve // leaflets will also be streched. Note that - // this is a simplifying assumption. - Vector xp(nsd); - + // this is a simplifying assumption. for (int ca = 0; ca < com_mod.tnNo; ca++) { // int mesh_ind = 0; // for (int ca = 0; ca < com_mod.msh[mesh_ind].nNo; ca++) { @@ -1096,9 +1052,7 @@ void uris_calc_sdf(ComMod& com_mod) { int jM = -1; Vector xb(nsd); for (int iM = 0; iM < uris_obj.nFa; iM++) { - auto& mesh = uris_obj.msh[iM]; - for (int e = 0; e < mesh.nEl; e++) { xb = 0.0; for (int a = 0; a < mesh.eNoN; a++) { @@ -1143,9 +1097,9 @@ void uris_calc_sdf(ComMod& com_mod) { for (int a = 0; a < mesh.eNoN; a++) { for (int i = 0; i < nsd - 1; i++) { - double factor = mesh.Nx(i,a,0); - for (int j = 0; j < nsd; j++) - xXi(j,i) += lX(j,a) * factor; + for (int j = 0; j < nsd; j++) { + xXi(j,i) += lX(j,a) * mesh.Nx(i,a,0); + } } } @@ -1163,64 +1117,52 @@ void uris_calc_sdf(ComMod& com_mod) { // [HZ] Improved implementation for SDF sign double sdf_sign = 1.0; - if (uris_obj.clsFlg) { - auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); - if (dot_nrm < 0.0 && dotP < 0.0) { - sdf_sign = -1.0; - } else { - sdf_sign = 1.0; - } - } else { - auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); - // if (dotP < 0.0) { - if (dot_nrm < 0.0 && dotP < 0.0) { - sdf_sign = -1.0; - } else { - sdf_sign = 1.0; - } - } + auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); + if (dot_nrm < 0.0 && dotP < 0.0) { + sdf_sign = -1.0; + } // uris_obj.sdf[com_mod.msh[mesh_ind].gN(ca)] = sdf_sign * minS; uris_obj.sdf[ca] = sdf_sign * minS; - Vector xp_plane(nsd), E1(nsd), E2(nsd), v(nsd), u_interp(nsd); - - // We need to interpolate valve velocity - // project xp onto the triangle plane - for (int i = 0; i < nsd; i++) { - xp_plane(i) = xp(i) - dotP * nV(i); - } - // compute barycentric (\xi, \eta) via the planar solve - for (int i = 0; i < nsd; i++) { - E1(i) = lX(i,1) - lX(i,0); - E2(i) = lX(i,2) - lX(i,0); - v(i) = xp_plane(i) - lX(i,0); - } - auto g11 = utils::norm(E1, E1); - auto g12 = utils::norm(E1, E2); - auto g22 = utils::norm(E2, E2); - auto b1 = utils::norm(v, E1); - auto b2 = utils::norm(v, E2); - double det = g11 * g22 - g12 * g12; - - double xi = (g22 * b1 - g12 * b2) / det; - double eta = (g11 * b2 - g12 * b1) / det; - - // shape functions: - double N1 = 1.0 - xi - eta; - double N2 = xi; - double N3 = eta; - - // interpolate the valve velocity: - for (int i = 0; i < nsd; i++) { - u_interp(i) = N1 * uris_obj.v(i, mesh.IEN(0,Ec)) + - N2 * uris_obj.v(i, mesh.IEN(1,Ec)) + - N3 * uris_obj.v(i, mesh.IEN(2,Ec)); - } - - // Store interpolated velocity - for (int i = 0; i < nsd; i++) { - uris_obj.sdf_t(i,ca) = u_interp(i); + if (uris_obj.use_valve_velocity) { + // We need to interpolate valve velocity + // project xp onto the triangle plane + for (int i = 0; i < nsd; i++) { + xp_plane(i) = xp(i) - dotP * nV(i); + } + // compute barycentric (\xi, \eta) via the planar solve + for (int i = 0; i < nsd; i++) { + E1(i) = lX(i,1) - lX(i,0); + E2(i) = lX(i,2) - lX(i,0); + v(i) = xp_plane(i) - lX(i,0); + } + auto g11 = utils::norm(E1, E1); + auto g12 = utils::norm(E1, E2); + auto g22 = utils::norm(E2, E2); + auto b1 = utils::norm(v, E1); + auto b2 = utils::norm(v, E2); + double det = g11 * g22 - g12 * g12; + + double xi = (g22 * b1 - g12 * b2) / det; + double eta = (g11 * b2 - g12 * b1) / det; + + // shape functions: + double N1 = 1.0 - xi - eta; + double N2 = xi; + double N3 = eta; + + // interpolate the valve velocity: + for (int i = 0; i < nsd; i++) { + u_interp(i) = N1 * uris_obj.v(i, mesh.IEN(0,Ec)) + + N2 * uris_obj.v(i, mesh.IEN(1,Ec)) + + N3 * uris_obj.v(i, mesh.IEN(2,Ec)); + } + + // Store interpolated velocity + for (int i = 0; i < nsd; i++) { + uris_obj.sdf_t(i,ca) = u_interp(i); + } } } @@ -1231,35 +1173,24 @@ void uris_calc_sdf(ComMod& com_mod) { for (int iUris = 0; iUris < nUris; iUris++) { auto& uris_obj = uris[iUris]; - if (!uris_obj.scaffold_flag || uris_obj.sdf_scaffold.size() > 0) { + if (!uris_obj.scaffold_flag || com_mod.cTS > 1) { continue; } auto& scaffold_mesh = uris_obj.scaffold_mesh; - - if (uris_obj.sdf_scaffold.size() <= 0) { - uris_obj.sdf_scaffold.resize(com_mod.tnNo); - uris_obj.sdf_scaffold = 0.0; - } uris_obj.sdf_scaffold = uris_obj.sdf_default; - Vector minb_scaf(nsd); - Vector maxb_scaf(nsd); - Vector extra_scaf(nsd); for (int i = 0; i < nsd; i++) { minb_scaf(i) = std::numeric_limits::max(); maxb_scaf(i) = std::numeric_limits::lowest(); } // For each coordinate dimension, find the minimum and maximum in uris_obj.x. - double extra_scaf_val = 0.1; // [HZ] The BBox is 10% larger than the actual valve, default is 0.1 for (int i = 0; i < nsd; i++) { for (int j = 0; j < scaffold_mesh.x.ncols(); j++) { double val = scaffold_mesh.x(i,j); - if (val < minb_scaf(i)) - minb_scaf(i) = val; - if (val > maxb_scaf(i)) - maxb_scaf(i) = val; + minb_scaf(i) = std::min(minb_scaf(i), val); + maxb_scaf(i) = std::max(maxb_scaf(i), val); } extra_scaf(i) = (maxb_scaf(i) - minb_scaf(i)) * extra_scaf_val; } From 320eb44ee59a425676336bfecf179d4e4657d574 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Tue, 20 Jan 2026 18:19:12 -0800 Subject: [PATCH 14/19] Adding explicit geometric coupling for FSI problems --- Code/Source/solver/ComMod.h | 3 ++ Code/Source/solver/Parameters.cpp | 1 + Code/Source/solver/Parameters.h | 1 + Code/Source/solver/distribute.cpp | 1 + Code/Source/solver/pic.cpp | 63 +++++++++++++++++++++---------- Code/Source/solver/read_files.cpp | 5 ++- 6 files changed, 54 insertions(+), 20 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index adb2daf30..1c7e652ad 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1202,6 +1202,9 @@ class eqType /// @brief URIS: Outputs std::vector outURIS; + /// @brief Explicit geometry coupling + bool expl_geom_cpl = false; + /// @brief Body force associated with this equation std::vector bf; }; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 9ec827ec7..41790adb5 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2037,6 +2037,7 @@ EquationParameters::EquationParameters() set_parameter("Tolerance", 0.5, !required, tolerance); set_parameter("Use_taylor_hood_type_basis", false, !required, use_taylor_hood_type_basis); + set_parameter("Explicit_geometric_coupling", false, !required, explicit_geometric_coupling); } diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 1de79c085..556946a60 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1340,6 +1340,7 @@ class EquationParameters : public ParameterLists Parameter type; Parameter use_taylor_hood_type_basis; + Parameter explicit_geometric_coupling; // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman Parameter inverse_darcy_permeability; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index e3a941276..c7f0186e1 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1505,6 +1505,7 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: cm.bcast(cm_mod, &lEq.tol); cm.bcast(cm_mod, &lEq.useTLS); cm.bcast(cm_mod, &lEq.assmTLS); + cm.bcast(cm_mod, &lEq.expl_geom_cpl); #ifdef dist_eq dmsg << "lEq.nOutput: " << lEq.nOutput; diff --git a/Code/Source/solver/pic.cpp b/Code/Source/solver/pic.cpp index 34a274def..0c62e05fd 100644 --- a/Code/Source/solver/pic.cpp +++ b/Code/Source/solver/pic.cpp @@ -80,6 +80,7 @@ void picc(Simulation* simulation) int s = eq.s; int e = eq.e; + std::array coef; coef[0] = eq.gam * dt; coef[1] = eq.beta*dt*dt; @@ -143,7 +144,7 @@ void picc(Simulation* simulation) pic_eth(simulation); } - if (eq.phys == Equation_FSI) { + if (eq.phys == Equation_FSI && !eq.expl_geom_cpl) { int s = com_mod.eq[1].s; int e = com_mod.eq[1].e; #ifdef debug_picc @@ -263,31 +264,55 @@ void picc(Simulation* simulation) //if (ALL(eq.ok)) RETURN if (eq.coupled) { - cEq = cEq + 1; - #ifdef debug_picc - dmsg << "eq " << " coupled "; - dmsg << "1st update cEq: " << cEq; - #endif - - auto& eqs = com_mod.eq; - if (std::count_if(eqs.begin(),eqs.end(),[](eqType& eq){return !eq.coupled || eq.ok;}) == eqs.size()) { - while (cEq < com_mod.nEq) { - if (!eqs[cEq].coupled) { - break; - } + if (eq.expl_geom_cpl) { + if (eq.ok) { cEq = cEq + 1; + if (cEq >= com_mod.nEq) { + cEq = 0; + } } - - } else { - if (cEq >= com_mod.nEq) { - cEq = 0; + if (eq.ok && eq.phys == Equation_FSI) { + int s = com_mod.eq[1].s; + int e = com_mod.eq[1].e; + for (int Ac = 0; Ac < tnNo; Ac++) { + if (all_fun::is_domain(com_mod, eq, Ac, Equation_struct) || + all_fun::is_domain(com_mod, eq, Ac, Equation_ustruct) || + all_fun::is_domain(com_mod, eq, Ac, Equation_lElas)) { + for (int i = 0; i < e-s+1; i++) { + An(i+s,Ac) = An(i,Ac); + Yn(i+s,Ac) = Yn(i,Ac); + Dn(i+s,Ac) = Dn(i,Ac); + } + } + } } + } else { + cEq = cEq + 1; + #ifdef debug_picc + dmsg << "eq " << " coupled "; + dmsg << "1st update cEq: " << cEq; + #endif + + auto& eqs = com_mod.eq; + if (std::count_if(eqs.begin(),eqs.end(),[](eqType& eq){return !eq.coupled || eq.ok;}) == eqs.size()) { + while (cEq < com_mod.nEq) { + if (!eqs[cEq].coupled) { + break; + } + cEq = cEq + 1; + } - while (!eqs[cEq].coupled) { - cEq = cEq + 1; + } else { if (cEq >= com_mod.nEq) { cEq = 0; } + + while (!eqs[cEq].coupled) { + cEq = cEq + 1; + if (cEq >= com_mod.nEq) { + cEq = 0; + } + } } } diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 9a644e3db..07112afc3 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1381,6 +1381,7 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) lEq.minItr = eq_params->min_iterations.value(); lEq.maxItr = eq_params->max_iterations.value(); lEq.tol = eq_params->tolerance.value(); + lEq.expl_geom_cpl = eq_params->explicit_geometric_coupling.value(); // Initialize coupled BC. // @@ -1796,7 +1797,9 @@ void read_files(Simulation* simulation, const std::string& file_name) if (eq.phys == EquationType::phys_mesh) { if (!com_mod.mvMsh) { throw std::runtime_error("mesh equation can only be specified after FSI equation"); - } + } + // Use the explicit geometry coupling flag of the FSI equation. + eq.expl_geom_cpl = com_mod.eq[0].expl_geom_cpl; } } #ifdef debug_read_files From 79fb10d69aac4a905a329ce2861f7d069660aab0 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Thu, 22 Jan 2026 12:14:41 -0800 Subject: [PATCH 15/19] Fixed the wrong sign issue in sign distance function --- Code/Source/solver/uris.cpp | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 3f3fcb50e..d55a99df1 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -130,7 +130,8 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } // If the uris has passed the closing state - if (uris_obj.cnt > uris_obj.DxClose.nrows()) { + int close_cnt = 1; // The uris is considered to be closed if within this number of states + if (uris_obj.cnt > close_cnt*uris_obj.DxClose.nrows()) { if (uris_obj.meanPD > uris_obj.meanPU) { uris_obj.cnt = 1; uris_obj.clsFlg = false; @@ -226,7 +227,8 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } // If the uris has passed the open state - if (uris_obj.cnt > uris_obj.DxOpen.nrows()) { + int open_cnt = 1; // The uris is considered to be open if within this number of states + if (uris_obj.cnt > open_cnt*uris_obj.DxOpen.nrows()) { if (meanV < 0.0) { uris_obj.cnt = 1; uris_obj.clsFlg = true; @@ -1117,10 +1119,20 @@ void uris_calc_sdf(ComMod& com_mod) { // [HZ] Improved implementation for SDF sign double sdf_sign = 1.0; - auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); - if (dot_nrm < 0.0 && dotP < 0.0) { - sdf_sign = -1.0; - } + if (uris_obj.clsFlg) { + auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); + if (dot_nrm < 0.0 && dotP < 0.0) { + sdf_sign = -1.0; + } else { + sdf_sign = 1.0; + } + } else { + if (dotP < 0.0) { + sdf_sign = -1.0; + } else { + sdf_sign = 1.0; + } + } // uris_obj.sdf[com_mod.msh[mesh_ind].gN(ca)] = sdf_sign * minS; uris_obj.sdf[ca] = sdf_sign * minS; From 59928c296d5f2b57ea30f17db2245b5c53b64710 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Fri, 23 Jan 2026 23:50:14 -0800 Subject: [PATCH 16/19] Adding a parameter in the uris valve to allow reversed surface normal --- Code/Source/solver/ComMod.h | 12 +++++++++--- Code/Source/solver/Parameters.cpp | 2 ++ Code/Source/solver/Parameters.h | 3 ++- Code/Source/solver/distribute.cpp | 2 ++ Code/Source/solver/uris.cpp | 20 ++++++++++++++------ 5 files changed, 29 insertions(+), 10 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 1c7e652ad..e22f319f7 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1503,6 +1503,9 @@ class urisType // Iteration count. int cnt = 1000000; + // pressurization time + double pressurization_time = 0.0; + // URIS: signed distance function of each node to the uris valves (1D array). Vector dirac_delta_func; @@ -1524,16 +1527,19 @@ class urisType // Relaxation factor to compute weighted averages of pressure values. double relax_factor = 0.5; - // Array to store the fluid mesh elements that the uris node is in (2D array). - Array elemId; - // Array to count how many times a uris node is found in the fluid mesh of a processor (1D array). Vector elemCounter; + // Array to store the fluid mesh elements that the uris node is in (2D array). + Array elemId; + // Derived type variables // IB meshes std::vector msh; + // Reserved leaflet surface normal + bool reverse_normal = false; + // Scaffold mesh bool scaffold_flag = false; mshType scaffold_mesh; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 41790adb5..546f8cace 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2720,6 +2720,8 @@ URISMeshParameters::URISMeshParameters() set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); set_parameter("Scaffold_file_path", "", !required, scaffold_file_path); + set_parameter("Pressurization_time", 0.0, !required, pressurization_time); + set_parameter("Reverse_surface_normal", false, !required, reverse_surface_normal); } void URISMeshParameters::print_parameters() diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 556946a60..79bc5ae09 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1616,7 +1616,8 @@ class URISMeshParameters : public ParameterLists Parameter valve_starts_as_closed; // Whether the valve starts as closed Parameter positive_flow_normal_file_path; // File path for the positive flow normal Parameter scaffold_file_path; // File path for the scaffold - + Parameter pressurization_time; // Time for pressurization + Parameter reverse_surface_normal; // Whether to reverse the surface normal vector }; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index c7f0186e1..cd682cfbd 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1171,12 +1171,14 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].sdf_deps_scaffold); cm.bcast(cm_mod, &uris[iUris].resistance); cm.bcast(cm_mod, &uris[iUris].resistance_close); + cm.bcast(cm_mod, &uris[iUris].pressurization_time); cm.bcast(cm_mod, &uris[iUris].clsFlg); cm.bcast(cm_mod, &uris[iUris].cnt); cm.bcast(cm_mod, &uris[iUris].scF); cm.bcast(cm_mod, uris[iUris].nrm); cm.bcast(cm_mod, &uris[iUris].scaffold_flag); + cm.bcast(cm_mod, &uris[iUris].reverse_normal); cm.bcast(cm_mod, &uris[iUris].use_valve_velocity); if (uris[iUris].scaffold_flag) { cm.bcast(cm_mod, &uris[iUris].scaffold_mesh.lShl); diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index d55a99df1..7def80fcb 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -130,8 +130,10 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } // If the uris has passed the closing state - int close_cnt = 1; // The uris is considered to be closed if within this number of states - if (uris_obj.cnt > close_cnt*uris_obj.DxClose.nrows()) { + int stab_cycle = 1; // The uris is considered to be closed if within this number of stages + // std::cout << "rank: " << cm.idcm() << " pressurization time: " << uris_obj.pressurization_time << std::endl; + if (uris_obj.cnt > stab_cycle*uris_obj.DxClose.nrows() + && com_mod.cTS*com_mod.dt > uris_obj.pressurization_time) { if (uris_obj.meanPD > uris_obj.meanPU) { uris_obj.cnt = 1; uris_obj.clsFlg = false; @@ -227,8 +229,9 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } // If the uris has passed the open state - int open_cnt = 1; // The uris is considered to be open if within this number of states - if (uris_obj.cnt > open_cnt*uris_obj.DxOpen.nrows()) { + int stab_cycle = 1; // The uris is considered to be open if within this number of stages + if (uris_obj.cnt > stab_cycle*uris_obj.DxOpen.nrows() + && com_mod.cTS*com_mod.dt > uris_obj.pressurization_time) { if (meanV < 0.0) { uris_obj.cnt = 1; uris_obj.clsFlg = true; @@ -492,7 +495,6 @@ void uris_read_msh(Simulation* simulation) { uris_obj.nFa = param->URIS_face_parameters.size(); uris_obj.msh.resize(uris_obj.nFa); uris_obj.nrm.resize(nsd); - Array gX(0,0); std::string positive_flow_normal_file_path = param->positive_flow_normal_file_path(); @@ -548,6 +550,8 @@ void uris_read_msh(Simulation* simulation) { uris_obj.resistance_close = param->resistance_close(); uris_obj.clsFlg = param->valve_starts_as_closed(); uris_obj.use_valve_velocity = param->use_valve_velocity(); + uris_obj.reverse_normal = param->reverse_surface_normal(); + uris_obj.pressurization_time = param->pressurization_time(); for (int iM = 0; iM < uris_obj.nFa; iM++) { auto mesh_param = param->URIS_face_parameters[iM]; @@ -1107,7 +1111,11 @@ void uris_calc_sdf(ComMod& com_mod) { auto nV = utils::cross(xXi); auto Jac = sqrt(utils::norm(nV)); - nV = nV / Jac; + if (uris_obj.reverse_normal) { + nV = -nV / Jac; + } else { + nV = nV / Jac; + } auto dotP = utils::norm(xp-xb, nV); // if (dotP < 0.0) { From 614b9c60691e12a1fd595a734f2a6aaca1bd27cc Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Wed, 11 Feb 2026 02:28:28 -0800 Subject: [PATCH 17/19] Update URIS test case result --- .../pipe_RCR_weak_dir_3d/lumen_inlet.fcs | 21 ------------------- tests/cases/uris/pipe_uris_3d/result_005.vtu | 4 ++-- tests/cases/uris/pipe_uris_3d/solver.xml | 12 ++++++----- 3 files changed, 9 insertions(+), 28 deletions(-) delete mode 100755 tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs diff --git a/tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs b/tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs deleted file mode 100755 index 844483140..000000000 --- a/tests/cases/fluid/pipe_RCR_weak_dir_3d/lumen_inlet.fcs +++ /dev/null @@ -1,21 +0,0 @@ - 0.0000 - 1.0000 - 0.000000E+00 - 0.000000E+00 - 16 - -6.283185E+01 0.000000E+00 - 6.263025E+01 -1.868984E-07 - -1.851028E-07 1.979808E-15 - -3.293147E-07 1.821408E-07 - 2.967629E-07 4.499564E-17 - -7.134611E-07 -1.729158E-07 - -1.667927E-07 5.199497E-16 - -1.920733E-06 1.597801E-07 - 5.066059E-08 -2.587250E-16 - 9.364629E-07 -1.435117E-07 - -1.344913E-07 8.999129E-17 - -9.410540E-07 1.250477E-07 - 1.921845E-07 2.874722E-17 - 9.540060E-07 -1.054107E-07 - -9.547580E-08 -1.106526E-16 - 1.016783E-06 8.563008E-08 diff --git a/tests/cases/uris/pipe_uris_3d/result_005.vtu b/tests/cases/uris/pipe_uris_3d/result_005.vtu index 025a881eb..c32413555 100644 --- a/tests/cases/uris/pipe_uris_3d/result_005.vtu +++ b/tests/cases/uris/pipe_uris_3d/result_005.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:884dedfcfb99696dfac3f47521f52b3eea2663d9b22a7564eddce4a9035ab44b -size 734173 +oid sha256:b3f9903f54a7ce1ae0bc535af18448398a86fc068bb08c7d7dbff02ed1d2bd5c +size 734957 diff --git a/tests/cases/uris/pipe_uris_3d/solver.xml b/tests/cases/uris/pipe_uris_3d/solver.xml index ccd61e375..2719e23f9 100644 --- a/tests/cases/uris/pipe_uris_3d/solver.xml +++ b/tests/cases/uris/pipe_uris_3d/solver.xml @@ -82,6 +82,7 @@ 1.0 0.25 + 1.0e5 meshes/normal.dat @@ -91,7 +92,9 @@ true 1 10 - 1e-4 + 1e-10 + true + fluid @@ -106,8 +109,7 @@ fsils - 1e-4 - 1e-4 + 1e-10 100 50 @@ -147,14 +149,14 @@ true 1 10 - 1e-6 + 1e-10 0. fsils - 1e-4 + 1e-10 From 216eb0e09082f67243226656da188aa63ed76bac Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Wed, 11 Feb 2026 02:43:59 -0800 Subject: [PATCH 18/19] Fixed bugs when restarting RIS valve simulation --- Code/Source/solver/ComMod.h | 11 ++--- Code/Source/solver/distribute.cpp | 2 - Code/Source/solver/fluid.cpp | 28 +------------ Code/Source/solver/fsi.cpp | 37 +--------------- Code/Source/solver/initialize.cpp | 2 + Code/Source/solver/uris.cpp | 70 ++++++------------------------- 6 files changed, 21 insertions(+), 129 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 1353bfdd2..2cb925269 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1209,7 +1209,7 @@ class eqType /// @brief URIS: Outputs std::vector outURIS; - /// @brief Explicit geometry coupling + /// @brief Explicit geometry coupling for FSI problem bool expl_geom_cpl = false; /// @brief Body force associated with this equation @@ -1465,9 +1465,6 @@ class urisType // Position coordinates of last time step. Array x_prev; - // Position coordinates of next time step. - Array x_next; - // Time derivative of the position coordinates (2D array). Array v; @@ -1495,6 +1492,9 @@ class urisType // Whether to use the valve velocity. bool use_valve_velocity = false; + // Reverse the surface normal vector + bool reverse_normal = false; + // Displacements of the valve when it opens (3D array). Array3 DxOpen; @@ -1544,9 +1544,6 @@ class urisType // IB meshes std::vector msh; - // Reserved leaflet surface normal - bool reverse_normal = false; - // Scaffold mesh bool scaffold_flag = false; mshType scaffold_mesh; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index daae3a686..819cbd85e 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1360,7 +1360,6 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { uris[iUris].DxClose.resize(CloseN(iUris), com_mod.nsd, uris[iUris].tnNo); uris[iUris].x.resize(com_mod.nsd, uris[iUris].tnNo); uris[iUris].x_prev.resize(com_mod.nsd, uris[iUris].tnNo); - uris[iUris].x_next.resize(com_mod.nsd, uris[iUris].tnNo); uris[iUris].v.resize(com_mod.nsd, uris[iUris].tnNo); uris[iUris].Yd.resize(com_mod.nsd, uris[iUris].tnNo); DxOpenFlat[iUris].resize(DxOpenFlatSize(iUris)); @@ -1371,7 +1370,6 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { for (int iUris = 0; iUris < com_mod.nUris; iUris++) { cm.bcast(cm_mod, uris[iUris].x); cm.bcast(cm_mod, uris[iUris].x_prev); - cm.bcast(cm_mod, uris[iUris].x_next); cm.bcast(cm_mod, uris[iUris].v); cm.bcast(cm_mod, uris[iUris].Yd); cm.bcast(cm_mod, DxOpenFlat[iUris]); diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index a4f40127b..ac1d1e56e 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -654,7 +654,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag dmsg << "w: " << w; #endif - // Plot the coordinates of the quad point in the current configuration + // Compute the resistance factor for the URIS if (com_mod.urisFlag) { Vector distSrf(com_mod.nUris); Vector distSrf_scaffold(com_mod.nUris); @@ -687,45 +687,19 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag if (distSrf_scaffold(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ (2*sdf_deps_temp*sdf_deps_temp); - // if (DDirTmp > DDir) {DDir = DDirTmp;} } } if (com_mod.uris[iUris].clsFlg) { sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; ris_resistance = com_mod.uris[iUris].resistance_close; - // Gradual transition using quadratic interpolation - // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { - // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps - // double progress = static_cast(com_mod.uris[iUris].cnt) / - // static_cast(com_mod.uris[iUris].DxClose.nrows()); - // // Quadratic interpolation: smooth transition - // double quad_progress = progress * progress; - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + - // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); - // } else { - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; - // } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; ris_resistance = com_mod.uris[iUris].resistance; - // // Gradual transition using quadratic interpolation - // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxOpen.nrows()) { - // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxOpen.nrows() steps - // double progress = static_cast(com_mod.uris[iUris].DxOpen.nrows() - com_mod.uris[iUris].cnt) / - // static_cast(com_mod.uris[iUris].DxOpen.nrows()); - // // Quadratic interpolation: smooth transition - // double quad_progress = progress * progress; - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + - // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); - // } else { - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps; - // } } if (distSrf(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ (2*sdf_deps_temp*sdf_deps_temp); - // if (DDirTmp > DDir) {DDir = DDirTmp;} } DDir = DDirTmp; } diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 61cf3ea8b..46a7d1663 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -62,15 +62,6 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar dmsg << "pS0.size(): " << pS0.size(); #endif - // if (com_mod.cm.idcm() == 0) { - // std::cout << "tnNo: " << com_mod.tnNo << std::endl; - // std::cout << "tDof: " << com_mod.tDof << std::endl; - // std::cout << "dof: " << com_mod.dof << std::endl; - // std::cout << "mesh name: " << lM.name << std::endl; - // std::cout << "mesh gnno: " << lM.gnNo << std::endl; - // std::cout << "mesh nEl: " << lM.nEl << std::endl; - // } - Vector ptr(eNoN); Array3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN); Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), @@ -203,7 +194,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar double w = fs_1[0].w(g) * Jac; - // Plot the coordinates of the quad point in the current configuration + // Compute the resistance factor for the URIS if (com_mod.urisFlag) { Vector distSrf(com_mod.nUris); Vector distSrf_scaffold(com_mod.nUris); @@ -236,45 +227,19 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar if (distSrf_scaffold(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf_scaffold(iUris)/sdf_deps_temp))/ (2*sdf_deps_temp*sdf_deps_temp); - // if (DDirTmp > DDir) {DDir = DDirTmp;} } } if (com_mod.uris[iUris].clsFlg) { sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; ris_resistance = com_mod.uris[iUris].resistance_close; - // Gradual transition using quadratic interpolation - // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxClose.nrows()) { - // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxClose.nrows() steps - // double progress = static_cast(com_mod.uris[iUris].cnt) / - // static_cast(com_mod.uris[iUris].DxClose.nrows()); - // // Quadratic interpolation: smooth transition - // double quad_progress = progress * progress; - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + - // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); - // } else { - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; - // } } else { sdf_deps_temp = com_mod.uris[iUris].sdf_deps; ris_resistance = com_mod.uris[iUris].resistance; - // // Gradual transition using quadratic interpolation - // if (com_mod.uris[iUris].cnt < com_mod.uris[iUris].DxOpen.nrows()) { - // // Quadratic interpolation: sdf_deps -> sdf_deps_close over DxOpen.nrows() steps - // double progress = static_cast(com_mod.uris[iUris].DxOpen.nrows() - com_mod.uris[iUris].cnt) / - // static_cast(com_mod.uris[iUris].DxOpen.nrows()); - // // Quadratic interpolation: smooth transition - // double quad_progress = progress * progress; - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps + - // quad_progress * (com_mod.uris[iUris].sdf_deps_close - com_mod.uris[iUris].sdf_deps); - // } else { - // sdf_deps_temp = com_mod.uris[iUris].sdf_deps; - // } } if (distSrf(iUris) <= sdf_deps_temp) { DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ (2*sdf_deps_temp*sdf_deps_temp); - // if (DDirTmp > DDir) {DDir = DDirTmp;} } DDir = DDirTmp; } diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 29e4b4909..2491350bb 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -852,9 +852,11 @@ void initialize(Simulation* simulation, Vector& timeP) for (int iUris = 0; iUris < com_mod.nUris; iUris++) { auto& uris_obj = com_mod.uris[iUris]; uris_obj.sdf.resize(com_mod.tnNo); + uris_obj.sdf = -1.0; uris_obj.sdf_t.resize(nsd, com_mod.tnNo); if (uris_obj.scaffold_flag) { uris_obj.sdf_scaffold.resize(com_mod.tnNo); + uris_obj.sdf_scaffold = -1.0; } } } diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 7def80fcb..32228bbaa 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -37,11 +37,6 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { auto& uris_obj = uris[iUris]; const int nsd = com_mod.nsd; - // const int cEq = com_mod.cEq; - - // auto& An = com_mod.An; - // auto& Ad = com_mod.Ad; - // auto& Dn = com_mod.Dn; auto& Yn = com_mod.Yn; // Let's conpute the mean pressure in the two regions of the fluid mesh @@ -59,9 +54,6 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { sUPS = 0.0; for (size_t j = 0; j < sUPS.size(); j++) { if (uris_obj.sdf(j) >= 0.0 && uris_obj.sdf(j) <= Deps) { - // Reverse the sdf distance for aortic valve - // [HZ] Adjust this value to be more flexible about the box - // if (uris_obj.sdf(j) < 0.0 && uris_obj.sdf(j) >= -Deps) { sUPS(0,j) = 1.0; } } @@ -76,8 +68,6 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { sDST = 0.0; for (size_t j = 0; j < sDST.size(); j++) { if (uris_obj.sdf(j) < 0.0 && uris_obj.sdf(j) >= -Deps) { - // Reverse the sdf distance for aortic valve - // if (uris_obj.sdf(j) >= 0.0 && uris_obj.sdf(j) <= Deps) { sDST(0,j) = 1.0; } } @@ -98,7 +88,6 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { int iEq = 0; int m = 1; int s = eq[iEq].s + nsd; - // int e = s + m - 1; Array tmpV(maxNSD, com_mod.tnNo); for (int j = 0; j < Yn.ncols(); j++) { @@ -130,8 +119,7 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } // If the uris has passed the closing state - int stab_cycle = 1; // The uris is considered to be closed if within this number of stages - // std::cout << "rank: " << cm.idcm() << " pressurization time: " << uris_obj.pressurization_time << std::endl; + int stab_cycle = 1; // The uris is considered to be closed if within this number of states if (uris_obj.cnt > stab_cycle*uris_obj.DxClose.nrows() && com_mod.cTS*com_mod.dt > uris_obj.pressurization_time) { if (uris_obj.meanPD > uris_obj.meanPU) { @@ -170,11 +158,6 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { auto& uris_obj = uris[iUris]; const int nsd = com_mod.nsd; - // const int cEq = com_mod.cEq; - - // auto& An = com_mod.An; - // auto& Ad = com_mod.Ad; - // auto& Dn = com_mod.Dn; auto& Yn = com_mod.Yn; // Let's compute the neighboring region below the valve normal. When @@ -188,8 +171,6 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { for (int i = 0; i < com_mod.tnNo; i++) { if (uris_obj.sdf(i) <= -Deps) { - // Reverse the sdf distance for aortic valve - // if (uris_obj.sdf(i) >= Deps) { sImm(0,i) = 1.0; } } @@ -203,7 +184,6 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { int m = nsd; int s = eq[iEq].s; - // int e = s + m - 1; Array tmpV(maxNSD, com_mod.tnNo); for (int i = 0; i < nsd; i++) { @@ -229,7 +209,7 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } // If the uris has passed the open state - int stab_cycle = 1; // The uris is considered to be open if within this number of stages + int stab_cycle = 1; // The uris is considered to be open if within this number of states if (uris_obj.cnt > stab_cycle*uris_obj.DxOpen.nrows() && com_mod.cTS*com_mod.dt > uris_obj.pressurization_time) { if (meanV < 0.0) { @@ -330,8 +310,6 @@ void uris_update_disp(ComMod& com_mod, CmMod& cm_mod) { for (int nd = 0; nd < uris_obj.tnNo; nd++) { double divisor = static_cast(std::max(1, uris_obj.elemCounter(nd))); - // Vector Yd_vec = uris_obj.Yd.col(nd) / div; - // uris_obj.Yd.set_col(nd, Yd_vec); for (int i = 0; i < nsd; i++) { uris_obj.Yd(i,nd) /= divisor; } @@ -563,7 +541,6 @@ void uris_read_msh(Simulation* simulation) { // Read mesh nodal coordinates and element connectivity. uris_read_sv(simulation, mesh, mesh_param); - // // [HZ] eType is not NA, neglecting this for now // IF (uris(iUris)%msh(iM)%eType .EQ. eType_NA) THEN // CALL READCCNE(lPN, uris(iUris)%msh(iM)) // END IF @@ -691,20 +668,15 @@ void uris_read_msh(Simulation* simulation) { } } uris_obj.tnNo = a; - // mesh.x.clear(); - // dispOpen.clear(); - // dispClose.clear(); } uris_obj.x.resize(nsd, uris_obj.tnNo); uris_obj.x_prev.resize(nsd, uris_obj.tnNo); - uris_obj.x_next.resize(nsd, uris_obj.tnNo); if (!uris_obj.clsFlg) { int nrowsDxOpen = uris_obj.DxOpen.nrows(); for (int i = 0; i < uris_obj.x.nrows(); i++) { for (int j = 0; j < uris_obj.x.ncols(); j++) { uris_obj.x(i,j) = uris_obj.DxOpen(nrowsDxOpen-1,i,j); uris_obj.x_prev(i,j) = uris_obj.DxOpen(nrowsDxOpen-1,i,j); - uris_obj.x_next(i,j) = uris_obj.DxOpen(nrowsDxOpen-1,i,j); } } } else { @@ -713,7 +685,6 @@ void uris_read_msh(Simulation* simulation) { for (int j = 0; j < uris_obj.x.ncols(); j++) { uris_obj.x(i,j) = uris_obj.DxClose(nrowsDxClose-1,i,j); uris_obj.x_prev(i,j) = uris_obj.DxClose(nrowsDxClose-1,i,j); - uris_obj.x_next(i,j) = uris_obj.DxClose(nrowsDxClose-1,i,j); } } } @@ -721,7 +692,6 @@ void uris_read_msh(Simulation* simulation) { uris_obj.v = 0.0; uris_obj.Yd.resize(nsd, uris_obj.tnNo); uris_obj.Yd = 0.0; - // gX.clear(); // Setting mesh.gN, mesh.lN parameter int b = 0; @@ -857,7 +827,6 @@ void uris_write_vtus(ComMod& com_mod) { } } - // [HZ] Need to check this if it's 1 or 0 outNames[0] = "coordinates"; outNames[1] = "URIS_displacement"; outNames[2] = "URIS_velocity"; @@ -879,7 +848,6 @@ void uris_write_vtus(ComMod& com_mod) { // Writing the position data int iOut = 0; int s = outS(iOut); - // int e = outS(iOut+1) - 1; int nSh = 0; Array tmpV(maxNSD, nNo); tmpV = 0.0; @@ -972,11 +940,6 @@ void uris_calc_sdf(ComMod& com_mod) { for (int i = 0; i < nsd; i++) { for (int j = 0; j < uris_obj.tnNo; j++) { if (uris_obj.use_valve_velocity) { - // uris_obj.v(i,j) = (uris_obj.x_next(i,j) - uris_obj.x_prev(i,j)) / (2*com_mod.dt); - - // uris_obj.v(i,j) = (-eq.af * (1 - eq.af) * uris_obj.x_prev(i,j) - // + eq.af * (1 - eq.af) * uris_obj.x_next(i,j) - // + (1 - 2 * eq.af) * uris_obj.x(i,j)) / com_mod.dt; uris_obj.v(i,j) = (uris_obj.x(i,j) - uris_obj.x_prev(i,j)) / com_mod.dt; } else { uris_obj.v(i,j) = 0.0; @@ -984,11 +947,14 @@ void uris_calc_sdf(ComMod& com_mod) { uris_obj.x_prev(i,j) = uris_obj.x(i,j); } } - + if (com_mod.cTS > 1 && cnt < uris_obj.cnt) { - // std::cout << "cts:" << com_mod.cTS << std::endl; uris_obj.sdf_t = 0.0; - continue; + if (com_mod.stFileFlag && uris_obj.sdf[0] == -1.0) { + uris_obj.sdf = uris_obj.sdf_default; + } else { + continue; + } } else { uris_obj.sdf = uris_obj.sdf_default; uris_obj.sdf_t = 0.0; @@ -1118,13 +1084,6 @@ void uris_calc_sdf(ComMod& com_mod) { } auto dotP = utils::norm(xp-xb, nV); - // if (dotP < 0.0) { - // dotP = -1.0; - // } else { - // dotP = 1.0; - // } - // uris_obj.sdf[ca] = dotP * minS; - // [HZ] Improved implementation for SDF sign double sdf_sign = 1.0; if (uris_obj.clsFlg) { @@ -1142,7 +1101,6 @@ void uris_calc_sdf(ComMod& com_mod) { } } - // uris_obj.sdf[com_mod.msh[mesh_ind].gN(ca)] = sdf_sign * minS; uris_obj.sdf[ca] = sdf_sign * minS; if (uris_obj.use_valve_velocity) { @@ -1194,7 +1152,11 @@ void uris_calc_sdf(ComMod& com_mod) { auto& uris_obj = uris[iUris]; if (!uris_obj.scaffold_flag || com_mod.cTS > 1) { - continue; + if (com_mod.stFileFlag && uris_obj.sdf_scaffold[0] == -1.0) { + uris_obj.sdf_scaffold = uris_obj.sdf_default; + } else { + continue; + } } auto& scaffold_mesh = uris_obj.scaffold_mesh; @@ -1222,9 +1184,6 @@ void uris_calc_sdf(ComMod& com_mod) { // this is a simplifying assumption. Vector xp_scaf(nsd); for (int ca = 0; ca < com_mod.tnNo; ca++) { - // int mesh_ind = 0; - // for (int ca = 0; ca < com_mod.msh[mesh_ind].nNo; ca++) { - // std::cout << "URIS scaffold SDF index: " << ca << std::endl; double minS_scaf = std::numeric_limits::max(); for (int i = 0; i < nsd; i++) { xp_scaf(i) = com_mod.x(i,ca); @@ -1265,7 +1224,6 @@ void uris_calc_sdf(ComMod& com_mod) { } } uris_obj.sdf_scaffold[ca] = minS_scaf; - // uris_obj.sdf_scaffold[com_mod.msh[mesh_ind].gN(ca)] = minS_scaf; } } } @@ -1399,8 +1357,6 @@ int same_side(Vector& v1, Vector& v2, Vector& v3, double dotV4 = utils::norm(N, v41); double dotP = utils::norm(N, vp1); // check if P and P4 are from the same side - // int sn = utils::sign(dotP); - // if (sn == 1) { if ((dotP >= 0 && dotV4 >= 0) || (dotP < 0 && dotV4 < 0)) { sameside = 1; } From 0560c0fc9d9ebe49ef1de26ffc5309ca22af7888 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Wed, 11 Feb 2026 10:34:09 -0800 Subject: [PATCH 19/19] Minor code clean up --- Code/Source/solver/uris.cpp | 4 ---- Code/Source/solver/vtk_xml.cpp | 6 +----- 2 files changed, 1 insertion(+), 9 deletions(-) diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 32228bbaa..cf9b146f8 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -15,7 +15,6 @@ #include "read_msh.h" #include "VtkData.h" - namespace uris { /// @brief This subroutine computes the mean pressure and flux on the @@ -62,7 +61,6 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { volU += all_fun::integ(com_mod, cm_mod, iM, sUPS); } - // Let's compute right side Array sDST(1,com_mod.tnNo); sDST = 0.0; @@ -238,8 +236,6 @@ void uris_update_disp(ComMod& com_mod, CmMod& cm_mod) { dmsg.banner(); #endif - // using namespace consts; - auto& cm = com_mod.cm; // auto& eq = com_mod.eq; auto& uris = com_mod.uris; diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 7911d1204..29fd7289a 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1290,6 +1290,7 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, const Array& lA, const Arrayset_point_data(outNames[iOut], tmpV); }