Add scaffold surface support for URIS valves#569
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #569 +/- ##
==========================================
+ Coverage 69.44% 69.53% +0.09%
==========================================
Files 181 181
Lines 34072 34162 +90
Branches 5930 5948 +18
==========================================
+ Hits 23662 23756 +94
+ Misses 10273 10268 -5
- Partials 137 138 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
ktbolt
left a comment
There was a problem hiding this comment.
@hanzhao2020 These changes look good !
However, there are still lots of messages printed to std::cout in uris.cpp; please replace those with DebugMsg calls.
@ktbolt Sounds good! Replaced all the URIS printed messages with |
ktbolt
left a comment
There was a problem hiding this comment.
@hanzhao2020 Good !
Note that if you would like to output useful run-time statistics for uris then you could create a log file for it, something like SimulationLogger.h, that could be activated from a parameter in the solver XML file.
|
@ktbolt Thanks for the suggestion, I'll keep that in mind! |
michelebucelli
left a comment
There was a problem hiding this comment.
Thanks @hanzhao2020! I left a few comments.
My one general concern is related to the comment I made to #566. What I meant by that comment is that, as far as I can see, in principle one could already have scaffolding with the current implementation (i.e. before this PR), by processing the input surface as follows:
- define a displacement field equal to zero on the scaffold mesh;
- append the scaffold mesh and the valve mesh (just append, no need to fix the connectivity in principle);
- give the resulting surface in input to URIS.
Is there some part of this that doesn't actually work as I would think?
| // Scaffold mesh flag | ||
| bool scaffold_flag = false; | ||
| // Scaffold mesh data | ||
| mshType scaffold_msh; | ||
| // Unsigned distance function (UDF) for scaffold mesh | ||
| Vector<double> scaffold_udf; | ||
| // Flag to indicate if the UDF for scaffold mesh is computed | ||
| bool scaffold_udf_computed = false; |
There was a problem hiding this comment.
I suggest to turn these comments into Doxygen-style documentation comments. I think it might be a good idea to do the same for other members of this class, even if they're otherwise untouched by this PR.
Additionally, I think it might be useful to document what the scaffold is and what its purpose is, for future reference and for other developers. I think a good way of doing that is to still use Doxygen comments, in one of the following ways:
- in the general documentation of the
urisTypeclass, e.g./// @brief Unfitted Resistive Immersed surface data type /// /// ...some generalities about the class... /// /// ### Scaffolding /// /// ...what a scaffold is, what is its purpose, what options are /// supported, and how to enable it...
- in a Doxygen section grouping these members, e.g.
/// @name Scaffold mesh /// ...what a scaffold is, what is its purpose, what options are /// supported, and how to enable it... /// @{ /// Scaffold mesh flag bool scaffold_flag = false; // ...other members... /// @}
I personally slightly prefer option 1, as it gives the chance to expand a bit on the general class documentation as well.
There was a problem hiding this comment.
@michelebucelli Thanks for catching that, I updated the comments to Doxygen-style comments.
| /// @brief Barycenter of a fixed shell element using mesh coordinates. | ||
| void mesh_element_barycenter(const mshType& mesh, const int Ec, Vector<double>& 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; | ||
| } |
There was a problem hiding this comment.
Since this function is only supposed to return a single value, I suggest to use that as the return value instead of an output argument, by changing the signature to
Vector<double>
mesh_element_barycenter(const mshType& msh, const int Ec)Additionally, I suggest renaming Ec to something more explicit like element_index, to make its meaning more obvious.
There was a problem hiding this comment.
@michelebucelli Updated. Thank you for the suggestions!
| if (scaffold_file_path != "") { | ||
| uris_obj.scaffold_flag = true; | ||
| uris_obj.scaffold_msh.lShl = true; | ||
| 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; | ||
| } | ||
| 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 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++) { | ||
| 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(); | ||
|
|
||
| # 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 | ||
| } |
There was a problem hiding this comment.
I feel like the readability of this piece of code would benefit from being encapsulated into a function with a meaningful name, such as load_scaffold_from_file. I imagine similar, if not identical, instructions are used to load the "normal" surfaces, in which case the function should be designed to be used in both cases.
In a dream world, this function would be a method of a class representing the surface itself (e.g. scaffold.load_from_file(...)).
There was a problem hiding this comment.
Yes, that's a good idea! Updated.
| /// @brief Find the closest fixed mesh element centroid to xp. | ||
| void uris_find_closest_mesh_centroid(const mshType& mesh, const Vector<double>& xp, | ||
| const int nsd, double& minS, int& Ec, | ||
| Vector<double>& xb) { |
There was a problem hiding this comment.
Here I suggest expanding the documentation to explain the meaning of each parameter using the Doxygen syntax @param[in] and @param[out].
I also suggest renaming the function to find_closest_element_centroid, since it seems more clear to me.
There was a problem hiding this comment.
Added comments to the input and output parameters with the Doxygen style comments.
| Vector<double> 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)); |
There was a problem hiding this comment.
This should be equivalent to
| const double dS = std::sqrt((xp - elem_centroid) * (xp - elem_centroid)); | |
| const double dS = utils::norm(xp - elem_centroid); |
@michelebucelli Thank you for the review and comments! Yes, with the current implementation on the main branch, it is already possible to include a scaffold by treating it as an additional valve leaflet and providing motion data that keeps it fixed throughout the simulation. However, this approach requires the user to provide motion data for both the opening and closing phases (even though they are identical for a stationary scaffold), and the motion files must have the same number of time steps as those of the valve leaflets. This adds unnecessary effort and could be confusing for users. Also, when the scaffold is treated as another valve leaflet, the SDF is recomputed during valve opening and closing, even though the scaffold itself does not move. These computations for scaffold SDF are unnecessary and introduce a bit additional computational cost. I think supporting the scaffold as an optional input, without requiring any motion data, provides a cleaner and more user-friendly interface while also avoiding these unnecessary SDF updates. |
Resolves #566
Release Notes
Scaffold_file_pathTesting
Code of Conduct & Contributing Guidelines