From f82c5750e55eca366542e9adb92303f15d4a175d Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Thu, 18 Jun 2026 18:11:55 -0700 Subject: [PATCH 1/5] Added support for scaffold surface in URIS valve --- Code/Source/solver/ComMod.h | 9 ++ Code/Source/solver/Parameters.cpp | 1 + Code/Source/solver/Parameters.h | 1 + Code/Source/solver/distribute.cpp | 27 ++++ Code/Source/solver/initialize.cpp | 5 + Code/Source/solver/uris.cpp | 232 ++++++++++++++++++++++-------- Code/Source/solver/vtk_xml.cpp | 23 ++- 7 files changed, 237 insertions(+), 61 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 86db290f4..c37c91a85 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -1553,6 +1553,15 @@ class urisType // IB meshes std::vector msh; + // Scaffold mesh flag + bool scaffold_flag = false; + // Scaffold mesh data + mshType scaffold_msh; + // Unsigned distance function (UDF) for scaffold mesh + Vector scaffold_udf; + // Flag to indicate if the UDF for scaffold mesh is computed + bool scaffold_udf_computed = false; + }; /// @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 e20dcf2f7..34644b2fe 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2990,6 +2990,7 @@ URISMeshParameters::URISMeshParameters() set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); set_parameter("Invert_normal", false, !required, invert_normal); set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); + set_parameter("Scaffold_file_path", "", !required, scaffold_file_path); set_parameter("Include_URIS_velocity", false, !required, include_uris_velocity); } diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index de383e80f..362c23643 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1786,6 +1786,7 @@ class URISMeshParameters : public ParameterLists Parameter valve_starts_as_closed; // Whether the valve starts as closed Parameter invert_normal; // Whether to invert the valve surface normal vector Parameter positive_flow_normal_file_path; // File path for the positive flow normal + Parameter scaffold_file_path; // File path for the valve scaffold mesh Parameter include_uris_velocity; // Whether to include the RIS velocity }; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 64acce552..6791e2d3b 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1180,10 +1180,37 @@ void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { cm.bcast(cm_mod, &uris[iUris].clsFlg); cm.bcast(cm_mod, &uris[iUris].invert_normal); cm.bcast(cm_mod, &uris[iUris].sdf_computed); + cm.bcast(cm_mod, &uris[iUris].scaffold_udf_computed); cm.bcast(cm_mod, &uris[iUris].include_uris_velocity); 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_msh.lShl); + cm.bcast(cm_mod, &uris[iUris].scaffold_msh.nEl); + cm.bcast(cm_mod, &uris[iUris].scaffold_msh.gnEl); + cm.bcast(cm_mod, &uris[iUris].scaffold_msh.eNoN); + cm.bcast(cm_mod, &uris[iUris].scaffold_msh.nNo); + cm.bcast(cm_mod, &uris[iUris].scaffold_msh.gnNo); + } + } + + if (cm.slv(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (uris[iUris].scaffold_flag) { + uris[iUris].scaffold_msh.x.resize(com_mod.nsd, uris[iUris].scaffold_msh.gnNo); + uris[iUris].scaffold_msh.IEN.resize(uris[iUris].scaffold_msh.eNoN, uris[iUris].scaffold_msh.gnEl); + } + } + } + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (uris[iUris].scaffold_flag) { + cm.bcast(cm_mod, uris[iUris].scaffold_msh.x); + cm.bcast(cm_mod, uris[iUris].scaffold_msh.IEN); + } } std::vector> lM_gN_flat(com_mod.nUris); diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 178b73fc1..07b924c66 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -888,6 +888,11 @@ void initialize(Simulation* simulation, Vector& timeP) uris_obj.sdf.resize(com_mod.tnNo); uris_obj.sdf = uris_obj.sdf_default; uris_obj.sdf_computed = false; + if (uris_obj.scaffold_flag && !uris_obj.scaffold_udf.allocated()) { + uris_obj.scaffold_udf.resize(com_mod.tnNo); + uris_obj.scaffold_udf = uris_obj.sdf_default; + uris_obj.scaffold_udf_computed = false; + } if (uris_obj.include_uris_velocity && !uris_obj.valve_velocity_fluid.allocated()) { uris_obj.valve_velocity_fluid.resize(nsd, com_mod.tnNo); uris_obj.valve_velocity_fluid = 0.0; diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index fe7c0fa90..1048f0c9a 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -18,6 +18,10 @@ namespace uris { +void uris_find_closest_mesh_centroid(const mshType& mesh, const Vector& xp, + const int nsd, double& minS, int& Ec, + Vector& xb); + /// @brief This subroutine computes the mean pressure and flux on the /// immersed surface void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionStates& solutions) { @@ -542,6 +546,38 @@ void uris_build_fluid_node_mask(ComMod& com_mod) { } } +/// @brief Build an expanded bounding box around coordinates. +void uris_compute_expanded_bbox(const Array& x, const int nsd, const double expansion, + Vector& minb, Vector& maxb) { + minb.resize(nsd); + maxb.resize(nsd); + + for (int i = 0; i < nsd; i++) { + double min_val = std::numeric_limits::max(); + double max_val = std::numeric_limits::lowest(); + for (int j = 0; j < x.ncols(); j++) { + const double val = x(i,j); + if (val < min_val) { min_val = val; } + if (val > max_val) { max_val = val; } + } + + const double extra = (max_val - min_val) * expansion; + minb(i) = min_val - extra; + maxb(i) = max_val + extra; + } +} + +/// @brief Check if a point lies inside a bounding box. +bool uris_point_in_bbox(const Vector& xp, const Vector& minb, + const Vector& maxb, const int nsd) { + for (int i = 0; i < nsd; ++i) { + if (xp(i) < minb(i) || xp(i) > maxb(i)) { + return false; + } + } + return true; +} + /// @brief Read the URIS mesh separately void uris_read_msh(Simulation* simulation) { #define n_debug_uris_read_msh @@ -607,6 +643,47 @@ void uris_read_msh(Simulation* simulation) { uris_obj.clsFlg = param->valve_starts_as_closed(); uris_obj.invert_normal = param->invert_normal(); uris_obj.include_uris_velocity = param->include_uris_velocity(); + std::string scaffold_file_path = param->scaffold_file_path(); + + if (scaffold_file_path != "") { + uris_obj.scaffold_flag = true; + uris_obj.scaffold_msh.lShl = true; + vtk_xml::read_vtu(scaffold_file_path, uris_obj.scaffold_msh); + // Scale the scaffold mesh coordinates by scF to match the URIS mesh scale + for (int a = 0; a < uris_obj.scaffold_msh.gnNo; a++) { + uris_obj.scaffold_msh.x.rcol(a) = uris_obj.scaffold_msh.x.rcol(a) * uris_obj.scF; + } + nn::select_ele(com_mod, uris_obj.scaffold_msh); + read_msh_ns::check_ien(simulation, uris_obj.scaffold_msh); + + int b = 0; + auto& scaffold_mesh = uris_obj.scaffold_msh; + 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_msh.gnNo << std::endl; + std::cout << "Scaffold mesh elements: " << uris_obj.scaffold_msh.gnEl << std::endl; + } // uris_obj.tnNo = 0; for (int iM = 0; iM < uris_obj.nFa; iM++) { @@ -1037,38 +1114,32 @@ void uris_calc_sdf(ComMod& com_mod) { uris_obj.valve_velocity_fluid = 0.0; } - if (cnt < uris_obj.cnt && uris_obj.sdf_computed) { + const bool compute_valve_sdf = !(cnt < uris_obj.cnt && uris_obj.sdf_computed); + const bool compute_scaffold_udf = uris_obj.scaffold_flag && !uris_obj.scaffold_udf_computed; + if (!compute_valve_sdf && !compute_scaffold_udf) { continue; } - uris_obj.sdf = uris_obj.sdf_default; - - if (cm.idcm() == 0) { - std::cout << "Recomputing SDF for " << uris_obj.name << std::endl; - } - // Each time when the URIS moves (open/close), we need to - // recompute the signed distance function. - // Find the bounding box of the valve, the BBox will be 10% larger - // than the actual valve. + const double bbox_expansion = 0.1; Vector minb(nsd); Vector maxb(nsd); - Vector extra(nsd); - for (int i = 0; i < nsd; i++) { - minb(i) = std::numeric_limits::max(); - maxb(i) = std::numeric_limits::lowest(); - } + if (compute_valve_sdf) { + uris_obj.sdf = uris_obj.sdf_default; - // 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; + if (cm.idcm() == 0) { + std::cout << "Recomputing SDF for " << uris_obj.name << std::endl; } - extra(i) = (maxb(i) - minb(i)) * extra_val; + + // The valve BBox is 10% larger than the current valve coordinates. + uris_compute_expanded_bbox(uris_obj.x, nsd, bbox_expansion, minb, maxb); + } + + auto& scaffold_mesh = uris_obj.scaffold_msh; + Vector minb_scaf(nsd); + Vector maxb_scaf(nsd); + if (compute_scaffold_udf) { + uris_obj.scaffold_udf = uris_obj.sdf_default; + uris_compute_expanded_bbox(scaffold_mesh.x, nsd, bbox_expansion, minb_scaf, maxb_scaf); } // The SDF is computed on the reference configuration, which @@ -1078,47 +1149,45 @@ void uris_calc_sdf(ComMod& com_mod) { // this is a simplifying assumption. Vector xp(nsd); for (int ca = 0; ca < com_mod.tnNo; ca++) { - double minS = std::numeric_limits::max(); xp = com_mod.x.rcol(ca); - // Check whether the node lies inside the expanded bounding box - bool inside_bbox = true; - for (int i = 0; i < nsd; ++i) { - const double lower = minb(i) - extra(i); - const double upper = maxb(i) + extra(i); - if (xp(i) < lower || xp(i) > upper) { - inside_bbox = false; - break; - } - } - if (!inside_bbox) { - continue; - } - if (!com_mod.urisFluidNodeMask[ca]) { continue; } - // This point is in the fluid domain and inside the BBox - // Find the closest URIS face centroid - int Ec = -1; - int jM = -1; - Vector xb(nsd); - Vector unitNormal(nsd); - uris_find_closest_face_centroid(uris_obj, xp, nsd, minS, Ec, jM, xb); - uris_face_unit_normal(uris_obj, nsd, jM, Ec, unitNormal); - const double dotp = (xp - xb) * unitNormal; - double sdf_sign = uris_compute_sdf_sign(uris_obj, xp, xb, dotp); - - uris_obj.sdf[ca] = sdf_sign * minS; - - if (uris_obj.include_uris_velocity) { - Vector interp_valve_vel(nsd); - uris_interp_valve_velocity(uris_obj, xp, nsd, jM, Ec, dotp, unitNormal, interp_valve_vel); - uris_obj.valve_velocity_fluid.rcol(ca) = interp_valve_vel; + if (compute_valve_sdf && uris_point_in_bbox(xp, minb, maxb, nsd)) { + double minS = std::numeric_limits::max(); + int Ec = -1; + int jM = -1; + Vector xb(nsd); + Vector unitNormal(nsd); + uris_find_closest_face_centroid(uris_obj, xp, nsd, minS, Ec, jM, xb); + uris_face_unit_normal(uris_obj, nsd, jM, Ec, unitNormal); + const double dotp = (xp - xb) * unitNormal; + const double sdf_sign = uris_compute_sdf_sign(uris_obj, xp, xb, dotp); + uris_obj.sdf[ca] = sdf_sign * minS; + + if (uris_obj.include_uris_velocity) { + Vector interp_valve_vel(nsd); + uris_interp_valve_velocity(uris_obj, xp, nsd, jM, Ec, dotp, unitNormal, interp_valve_vel); + uris_obj.valve_velocity_fluid.rcol(ca) = interp_valve_vel; + } } + if (compute_scaffold_udf && uris_point_in_bbox(xp, minb_scaf, maxb_scaf, nsd)) { + double minS_scaf = std::numeric_limits::max(); + int Ec = -1; + Vector xb(nsd); + uris_find_closest_mesh_centroid(scaffold_mesh, xp, nsd, minS_scaf, Ec, xb); + uris_obj.scaffold_udf[ca] = minS_scaf; + } } // ca: loop - uris_obj.sdf_computed = true; + + if (compute_valve_sdf) { + uris_obj.sdf_computed = true; + } + if (compute_scaffold_udf) { + uris_obj.scaffold_udf_computed = true; + } } // iUris: loop } @@ -1264,6 +1333,16 @@ void surface_element_barycenter(const urisType& uris_obj, int jM, int Ec, Vector xb = xb / mesh.eNoN; } +/// @brief Barycenter of a fixed shell element using mesh coordinates. +void mesh_element_barycenter(const mshType& mesh, const int Ec, Vector& xb) { + xb = 0.0; + for (int a = 0; a < mesh.eNoN; a++) { + const int Ac = mesh.IEN(a, Ec); + xb = xb + mesh.x.rcol(Ac); + } + xb = xb / mesh.eNoN; +} + /// @brief Find the closest URIS shell-element centroid to xp; writes that centroid to xb. void uris_find_closest_face_centroid(const urisType& uris_obj, const Vector& xp, const int nsd, double& minS, int& Ec, int& jM, @@ -1291,6 +1370,22 @@ void uris_find_closest_face_centroid(const urisType& uris_obj, const Vector& xp, + const int nsd, double& minS, int& Ec, + Vector& xb) { + Vector elem_centroid(nsd); + for (int e = 0; e < mesh.nEl; e++) { + mesh_element_barycenter(mesh, e, elem_centroid); + const double dS = std::sqrt((xp - elem_centroid) * (xp - elem_centroid)); + if (dS < minS) { + minS = dS; + Ec = e; + xb = elem_centroid; + } + } +} + /// @brief Unit normal of URIS face element (jM, Ec) from parametric tangents /// (optionally flipped by uris_obj.invert_normal). void uris_face_unit_normal(const urisType& uris_obj, const int nsd, const int jM, const int Ec, @@ -1421,15 +1516,20 @@ void eval_uris_ris_factors_quadrature(const ComMod& com_mod, const mshType& lM, } Vector dist_srf(nUris); + Vector dist_scaffold(nUris); Array valve_velocity(com_mod.nsd, nUris); for (int g = 0; g < fs.nG; g++) { dist_srf = 0.0; + dist_scaffold = 0.0; valve_velocity = 0.0; for (int a = 0; a < fs.eNoN; a++) { int Ac = lM.IEN(a,e); for (int iUris = 0; iUris < nUris; iUris++) { dist_srf(iUris) += fs.N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + if (com_mod.uris[iUris].scaffold_flag) { + dist_scaffold(iUris) += fs.N(a,g) * std::fabs(com_mod.uris[iUris].scaffold_udf(Ac)); + } if (com_mod.uris[iUris].include_uris_velocity) { valve_velocity.rcol(iUris) = valve_velocity.rcol(iUris) + fs.N(a,g) * com_mod.uris[iUris].valve_velocity_fluid.rcol(Ac); @@ -1439,9 +1539,11 @@ void eval_uris_ris_factors_quadrature(const ComMod& com_mod, const mshType& lM, double sdf_deps; double delta_eps; + double delta_eps_scaffold; for (int iUris = 0; iUris < nUris; iUris++) { sdf_deps = 0.0; delta_eps = 0.0; + delta_eps_scaffold = 0.0; double start_deps, end_deps; int n_steps; if (com_mod.uris[iUris].clsFlg) { @@ -1471,7 +1573,17 @@ void eval_uris_ris_factors_quadrature(const ComMod& com_mod, const mshType& lM, if (dist_srf(iUris) < sdf_deps && sdf_deps > 0.0) { delta_eps = (1 + cos(consts::pi*dist_srf(iUris)/sdf_deps))/(2*sdf_deps*sdf_deps); } - uris_factor_total_el(g) += com_mod.uris[iUris].resistance * delta_eps; + + if (com_mod.uris[iUris].scaffold_flag) { + // Compute the scaffold resistance factor based on the unsigned distance function (UDF) + // The scaffold surface uses the same thickness parameter as the closed valve surface + const double scaffold_deps = com_mod.uris[iUris].sdf_deps_close; + if (dist_scaffold(iUris) < scaffold_deps && scaffold_deps > 0.0) { + delta_eps_scaffold = (1 + cos(consts::pi*dist_scaffold(iUris)/scaffold_deps))/(2*scaffold_deps*scaffold_deps); + } + } + + uris_factor_total_el(g) += com_mod.uris[iUris].resistance * (delta_eps + delta_eps_scaffold); if (com_mod.uris[iUris].include_uris_velocity) { uris_valve_vel_term_total_el.rcol(g) = uris_valve_vel_term_total_el.rcol(g) @@ -1481,4 +1593,4 @@ void eval_uris_ris_factors_quadrature(const ComMod& com_mod, const mshType& lM, } // g: loop } -} \ No newline at end of file +} diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 0119a7bb2..4468cf95d 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1006,6 +1006,10 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b nOut = nOut + com_mod.nUris; outDof = outDof + com_mod.nUris; for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (com_mod.uris[iUris].scaffold_flag) { + nOut = nOut + 1; + outDof = outDof + 1; + } if (com_mod.uris[iUris].include_uris_velocity) { nOut = nOut + 1; outDof = outDof + nsd; @@ -1353,19 +1357,35 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b } if (com_mod.urisFlag) { + // SDF for each URIS for (int iUris = 0; iUris < com_mod.nUris; iUris++) { cOut = cOut + 1; int is = outS[cOut]; int ie = is; outS[cOut+1] = ie + 1; outNames[cOut] = "URIS_SDF_" + 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(Ac)); } } + // SDF for scaffold + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + 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_SCAF_UDF_" + 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].scaffold_udf(Ac)); + } + } + } + // Valve velocity for each URIS for (int iUris = 0; iUris < com_mod.nUris; iUris++) { if (com_mod.uris[iUris].include_uris_velocity) { cOut = cOut + 1; @@ -1373,6 +1393,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b int ie = is + nsd - 1; outS[cOut+1] = ie + 1; outNames[cOut] = "URIS_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++) { From a82d065b0ffee1fc1b6d75ceaea228d7321c57c4 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Fri, 19 Jun 2026 01:52:54 -0700 Subject: [PATCH 2/5] Updated URIS valve CFD test case to include valve scaffold surface mesh --- tests/cases/uris/pipe_uris_cfd/README.md | 1 + tests/cases/uris/pipe_uris_cfd/meshes/mv_scaffold.vtu | 3 +++ tests/cases/uris/pipe_uris_cfd/result_003.vtu | 4 ++-- tests/cases/uris/pipe_uris_cfd/solver.xml | 1 + 4 files changed, 7 insertions(+), 2 deletions(-) create mode 100644 tests/cases/uris/pipe_uris_cfd/meshes/mv_scaffold.vtu diff --git a/tests/cases/uris/pipe_uris_cfd/README.md b/tests/cases/uris/pipe_uris_cfd/README.md index d2a296f72..cdbf95f1d 100644 --- a/tests/cases/uris/pipe_uris_cfd/README.md +++ b/tests/cases/uris/pipe_uris_cfd/README.md @@ -18,6 +18,7 @@ The model parameters are specified in the `Add_URIS_mesh` sub-section false true meshes/normal.dat + meshes/mv_scaffold.vtu ``` diff --git a/tests/cases/uris/pipe_uris_cfd/meshes/mv_scaffold.vtu b/tests/cases/uris/pipe_uris_cfd/meshes/mv_scaffold.vtu new file mode 100644 index 000000000..42a8730bc --- /dev/null +++ b/tests/cases/uris/pipe_uris_cfd/meshes/mv_scaffold.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:763a2c73ebadfcb93586e4c4ea8d95edd9558506178cc37537fa5cf9347068e3 +size 14625 diff --git a/tests/cases/uris/pipe_uris_cfd/result_003.vtu b/tests/cases/uris/pipe_uris_cfd/result_003.vtu index 7a88a0dbc..49e6e6b21 100644 --- a/tests/cases/uris/pipe_uris_cfd/result_003.vtu +++ b/tests/cases/uris/pipe_uris_cfd/result_003.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:98385819e19b1e0244f83f4a8ccccba7dc7e1d94a74eb3e406ddedb70086b8c6 -size 1547683 +oid sha256:449a6ef3e203c65860eb50ed4d1b079b0d93641cbb25dac2b798408146974f16 +size 1556147 diff --git a/tests/cases/uris/pipe_uris_cfd/solver.xml b/tests/cases/uris/pipe_uris_cfd/solver.xml index bdf9dce88..bad576138 100644 --- a/tests/cases/uris/pipe_uris_cfd/solver.xml +++ b/tests/cases/uris/pipe_uris_cfd/solver.xml @@ -62,6 +62,7 @@ false true meshes/normal.dat + meshes/mv_scaffold.vtu From a6d5d9ab3ac818ed3c6a6ab7086060a8dc99c479 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Fri, 19 Jun 2026 01:56:47 -0700 Subject: [PATCH 3/5] Updated mesh-scale padding for degenerate URIS bounding box --- Code/Source/solver/uris.cpp | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 1048f0c9a..9d8d9e735 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -551,19 +551,33 @@ void uris_compute_expanded_bbox(const Array& x, const int nsd, const dou Vector& minb, Vector& maxb) { minb.resize(nsd); maxb.resize(nsd); + Vector min_val(nsd); + Vector max_val(nsd); for (int i = 0; i < nsd; i++) { - double min_val = std::numeric_limits::max(); - double max_val = std::numeric_limits::lowest(); + min_val(i) = std::numeric_limits::max(); + max_val(i) = std::numeric_limits::lowest(); for (int j = 0; j < x.ncols(); j++) { const double val = x(i,j); - if (val < min_val) { min_val = val; } - if (val > max_val) { max_val = val; } + if (val < min_val(i)) { min_val(i) = val; } + if (val > max_val(i)) { max_val(i) = val; } } + } + // Compute the diagonal length of the bounding box for use in degenerate cases where max_val == min_val + double diag_length = std::sqrt((max_val - min_val) * (max_val - min_val)); - const double extra = (max_val - min_val) * expansion; - minb(i) = min_val - extra; - maxb(i) = max_val + extra; + for (int i = 0; i < nsd; i++) { + double extra = 0.0; + if (max_val(i) > min_val(i)) { + extra = (max_val(i) - min_val(i)) * expansion; + } else if (max_val(i) == min_val(i)) { + // When points are all the same in this dimension, use the diagonal length of the bounding box + extra = expansion * (diag_length > 0.0 ? diag_length : 1.0); + } else { + throw std::runtime_error("Invalid bounding box: max_val < min_val for dimension " + std::to_string(i)); + } + minb(i) = min_val(i) - extra; + maxb(i) = max_val(i) + extra; } } @@ -1569,11 +1583,9 @@ void eval_uris_ris_factors_quadrature(const ComMod& com_mod, const mshType& lM, double progress = static_cast(com_mod.uris[iUris].cnt) / static_cast(n_steps); sdf_deps = start_deps + progress * (end_deps - start_deps); } - if (dist_srf(iUris) < sdf_deps && sdf_deps > 0.0) { delta_eps = (1 + cos(consts::pi*dist_srf(iUris)/sdf_deps))/(2*sdf_deps*sdf_deps); } - if (com_mod.uris[iUris].scaffold_flag) { // Compute the scaffold resistance factor based on the unsigned distance function (UDF) // The scaffold surface uses the same thickness parameter as the closed valve surface From 529c670d50c08b8cf856c996086f7d7fc79ec30a Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Fri, 19 Jun 2026 23:25:15 -0700 Subject: [PATCH 4/5] Updated exception check when reading scaffold meshes and debug info --- Code/Source/solver/uris.cpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 9d8d9e735..15fee28c7 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -662,7 +662,13 @@ void uris_read_msh(Simulation* simulation) { if (scaffold_file_path != "") { uris_obj.scaffold_flag = true; uris_obj.scaffold_msh.lShl = true; - vtk_xml::read_vtu(scaffold_file_path, uris_obj.scaffold_msh); + try { + vtk_xml::read_vtu(scaffold_file_path, uris_obj.scaffold_msh); + } catch (const std::exception& e) { + throw std::runtime_error("Failed to read URIS scaffold mesh for '" + uris_obj.name + "': " + e.what()); + } catch (...) { + throw std::runtime_error("Failed to read URIS scaffold mesh for '" + uris_obj.name + "'."); + } // Scale the scaffold mesh coordinates by scF to match the URIS mesh scale for (int a = 0; a < uris_obj.scaffold_msh.gnNo; a++) { uris_obj.scaffold_msh.x.rcol(a) = uris_obj.scaffold_msh.x.rcol(a) * uris_obj.scF; @@ -682,7 +688,7 @@ void uris_read_msh(Simulation* simulation) { scaffold_mesh.lN(b) = a; b++; } - // Remap msh%gIEN array + // Remap scaffold_mesh.gIEN to scaffold_mesh.IEN 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++) { @@ -907,7 +913,7 @@ void uris_read_msh(Simulation* simulation) { throw std::runtime_error("Mismatch in uris.tnNo. Correction needed."); } - // Remap msh%gIEN array + // Remap mesh.gIEN to mesh.IEN for (int iM = 0; iM < uris_obj.nFa; iM++) { auto& mesh = uris_obj.msh[iM]; mesh.nEl = mesh.gnEl; @@ -1138,12 +1144,10 @@ void uris_calc_sdf(ComMod& com_mod) { Vector minb(nsd); Vector maxb(nsd); if (compute_valve_sdf) { + #ifdef debug_uris_calc_sdf + dmsg << "Recomputing SDF for " << uris_obj.name; + #endif uris_obj.sdf = uris_obj.sdf_default; - - if (cm.idcm() == 0) { - std::cout << "Recomputing SDF for " << uris_obj.name << std::endl; - } - // The valve BBox is 10% larger than the current valve coordinates. uris_compute_expanded_bbox(uris_obj.x, nsd, bbox_expansion, minb, maxb); } From c1703ebb13b3a3453520a58c97ea325ce45ebed3 Mon Sep 17 00:00:00 2001 From: hanzhao2020 Date: Mon, 22 Jun 2026 17:08:54 -0700 Subject: [PATCH 5/5] Replaced URIS printed messages with DebugMsg --- Code/Source/solver/main.cpp | 5 +- Code/Source/solver/uris.cpp | 104 ++++++++++++++++++++++-------------- 2 files changed, 68 insertions(+), 41 deletions(-) diff --git a/Code/Source/solver/main.cpp b/Code/Source/solver/main.cpp index 4e778f03e..4cfb35fb4 100644 --- a/Code/Source/solver/main.cpp +++ b/Code/Source/solver/main.cpp @@ -523,9 +523,12 @@ void iterate_solution(Simulation* simulation) } else { uris::uris_meanv(com_mod, cm_mod, iUris, solutions); } + # ifdef debug_iterate_solution if (cm.mas(cm_mod)) { - std::cout << " URIS surface: " << com_mod.uris[iUris].name << ", count: " << com_mod.uris[iUris].cnt << std::endl; + dmsg << " URIS surface: " + com_mod.uris[iUris].name + ", count: " + + std::to_string(com_mod.uris[iUris].cnt) << std::endl; } + # endif } if (com_mod.mvMsh) { diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index 15fee28c7..e318574b6 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -58,12 +58,14 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionS double volU = 0.0; double volD = 0.0; - // if (cm.mas(cm_mod)) { - // std::cout << "Computing upstream region from SDF -" << meanp_sdf_outer_limit << " to -" - // << uris_obj.sdf_deps_close << " for: " << uris_obj.name << std::endl; - // std::cout << "Computing downstream region from SDF " << uris_obj.sdf_deps_close - // << " to " << meanp_sdf_outer_limit << " for: " << uris_obj.name << std::endl; - // } + # ifdef debug_uris_meanp + if (cm.mas(cm_mod)) { + dmsg << "Computing upstream region from SDF -" + std::to_string(meanp_sdf_outer_limit) + " to -" + + std::to_string(uris_obj.sdf_deps_close) + " for: " + uris_obj.name << std::endl; + dmsg << "Computing downstream region from SDF " + std::to_string(uris_obj.sdf_deps_close) + + " to " + std::to_string(meanp_sdf_outer_limit) + " for: " + uris_obj.name << std::endl; + } + # endif // Compute the upstream region: negative sdf side (opposite to valve normal), // outside resistance region @@ -96,10 +98,12 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionS } // Print volume messages. + # ifdef debug_uris_meanp if (cm.mas(cm_mod)) { - std::cout << "volume upstream " << volU << " for: " << uris_obj.name << std::endl; - std::cout << "volume downstream " << volD << " for: " << uris_obj.name << std::endl; + dmsg << "volume upstream " + std::to_string(volU) + " for: " + uris_obj.name << std::endl; + dmsg << "volume downstream " + std::to_string(volD) + " for: " + uris_obj.name << std::endl; } + # endif double meanPU = 0.0; double meanPD = 0.0; @@ -131,12 +135,14 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionS uris_obj.meanPD = uris_obj.relax_factor * meanPD + (1.0 - uris_obj.relax_factor) * uris_obj.meanPD; + # ifdef debug_uris_meanp if (cm.mas(cm_mod)) { - std::cout << "mean P upstream " << meanPU << " " << uris_obj.meanPU - << " for: " << uris_obj.name << std::endl; - std::cout << "mean P downstream " << meanPD << " " << uris_obj.meanPD - << " for: " << uris_obj.name << std::endl; + dmsg << "mean P upstream " + std::to_string(meanPU) + " " + + std::to_string(uris_obj.meanPU) + " for: " + uris_obj.name << std::endl; + dmsg << "mean P downstream " + std::to_string(meanPD) + " " + + std::to_string(uris_obj.meanPD) + " for: " + uris_obj.name << std::endl; } + #endif // If the uris has passed the closing state if (uris_obj.cnt > uris_obj.DxClose.nslices()) { @@ -144,17 +150,18 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionS uris_obj.cnt = 1; uris_obj.clsFlg = false; com_mod.urisActFlag = true; + # ifdef debug_uris_meanp if (cm.mas(cm_mod)) { - std::cout << "** Set urisCloseFlag to FALSE for: " - << uris_obj.name << std::endl; + dmsg << "** Set urisCloseFlag to FALSE for: " + uris_obj.name << std::endl; } + # endif } } + # ifdef debug_uris_meanp if (cm.mas(cm_mod)) { - std::cout << "urisCloseFlag is: " << uris_obj.clsFlg << " for: " - << uris_obj.name << std::endl; + dmsg << "urisCloseFlag is: " + std::to_string(uris_obj.clsFlg) + " for: " + uris_obj.name << std::endl; } - + # endif } /// @brief This subroutine computes the mean velocity in the fluid elements @@ -199,9 +206,11 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionS for (int iM = 0; iM < com_mod.nMsh; iM++) { volI += all_fun::integ(com_mod, cm_mod, iM, sImm, solutions); } + # ifdef debug_uris_meanv if (cm.mas(cm_mod)) { - std::cout << "volume inside " << volI << " for: " << uris_obj.name << std::endl; + dmsg << "volume inside " + std::to_string(volI) + " for: " + uris_obj.name << std::endl; } + # endif int m = nsd; int s = eq[iEq].s; @@ -225,9 +234,11 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionS meanV += all_fun::integ(com_mod, cm_mod, iM, tmpVNrm, solutions)/volI; } + # ifdef debug_uris_meanv if (cm.mas(cm_mod)) { - std::cout << "mean velocity: " << meanV << " for: " << uris_obj.name << std::endl; + dmsg << "mean velocity: " + std::to_string(meanV) + " for: " + uris_obj.name << std::endl; } + # endif // If the uris has passed the open state if (uris_obj.cnt > uris_obj.DxOpen.nslices()) { @@ -235,17 +246,18 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionS uris_obj.cnt = 1; uris_obj.clsFlg = true; com_mod.urisActFlag = true; + # ifdef debug_uris_meanv if (cm.mas(cm_mod)) { - std::cout << "** Set urisCloseFlag to TRUE for: " - << uris_obj.name << std::endl; + dmsg << "** Set urisCloseFlag to TRUE for: " + uris_obj.name << std::endl; } + # endif } } + # ifdef debug_uris_meanv if (cm.mas(cm_mod)) { - std::cout << "urisCloseFlag is: " << uris_obj.clsFlg << " for: " - << uris_obj.name << std::endl; + dmsg << "urisCloseFlag is: " + std::to_string(uris_obj.clsFlg) + " for: " + uris_obj.name << std::endl; } - + # endif } /// @brief This subroutine computes the displacement of the immersed @@ -609,15 +621,18 @@ void uris_read_msh(Simulation* simulation) { int nUris = simulation->parameters.URIS_mesh_parameters.size(); com_mod.nUris = nUris; - - std::cout << "Number of immersed surfaces for uris: " << nUris << std::endl; + # ifdef debug_uris_read_msh + dmsg << "Number of immersed surfaces for uris: " + std::to_string(nUris) << std::endl; + # endif uris.resize(nUris); for (int iUris = 0; iUris < nUris; iUris++) { auto param = simulation->parameters.URIS_mesh_parameters[iUris]; auto& uris_obj = uris[iUris]; uris_obj.name = param->name(); - std::cout << "** Reading URIS mesh: " << uris_obj.name << std::endl; + # ifdef debug_uris_read_msh + dmsg << "** Reading URIS mesh: " + uris_obj.name << std::endl; + # endif uris_obj.scF = param->mesh_scale_factor(); uris_obj.nFa = param->URIS_face_parameters.size(); @@ -699,10 +714,12 @@ 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_msh.gnNo << std::endl; - std::cout << "Scaffold mesh elements: " << uris_obj.scaffold_msh.gnEl << std::endl; + + # ifdef debug_uris_read_msh + dmsg << "Scaffold mesh is included for: " + uris_obj.name << std::endl; + dmsg << "Scaffold mesh nodes: " + std::to_string(uris_obj.scaffold_msh.gnNo) << std::endl; + dmsg << "Scaffold mesh elements: " + std::to_string(uris_obj.scaffold_msh.gnEl) << std::endl; + # endif } // uris_obj.tnNo = 0; @@ -712,7 +729,9 @@ void uris_read_msh(Simulation* simulation) { auto& mesh = uris_obj.msh[iM]; mesh.lShl = true; mesh.name = mesh_param->name(); - std::cout << "-- Reading URIS face: " << mesh.name << std::endl; + # ifdef debug_uris_read_msh + dmsg << "-- Reading URIS face: " + mesh.name << std::endl; + # endif // Read mesh nodal coordinates and element connectivity. uris_read_sv(simulation, mesh, mesh_param); @@ -731,8 +750,10 @@ void uris_read_msh(Simulation* simulation) { // err = " Failed to identify format of the uris mesh" // END IF - std::cout << "Number of uris nodes: " << mesh.gnNo << std::endl; - std::cout << "Number of uris elements: " << mesh.gnEl << std::endl; + # ifdef debug_uris_read_msh + dmsg << "Number of uris nodes: " + std::to_string(mesh.gnNo) << std::endl; + dmsg << "Number of uris elements: " + std::to_string(mesh.gnEl) << std::endl; + # endif // Read valve motion: note that this motion is defined on the // reference configuration @@ -929,18 +950,21 @@ void uris_read_msh(Simulation* simulation) { } if (uris_obj.nFa > 0) { - std::string msg = "Total number of uris nodes: " + std::to_string(uris_obj.tnNo); - std::cout << msg << std::endl; + # ifdef debug_uris_read_msh + dmsg << "Total number of uris nodes: " + std::to_string(uris_obj.tnNo) << std::endl; + # endif int total_nel = 0; for (int iM = 0; iM < uris_obj.nFa; iM++) { total_nel += uris_obj.msh[iM].nEl; } - msg = "Total number of uris elements: " + std::to_string(total_nel); - std::cout << msg << std::endl; + # ifdef debug_uris_read_msh + dmsg << "Total number of uris elements: " + std::to_string(total_nel) << std::endl; + # endif } } - std::cout << "URIS mesh data imported successfully." << std::endl; - + # ifdef debug_uris_read_msh + dmsg << "URIS mesh data imported successfully." << std::endl; + # endif } /// @brief Write URIS solution to a vtu file