From 32736c4bd973cac5ddaab498e3159772d2bf15e7 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Mon, 8 Dec 2025 13:00:37 +0100 Subject: [PATCH 01/22] add a convenience function for floating point equality --- src/t8.h | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/t8.h b/src/t8.h index 2b20d618c9..0e16b021c7 100644 --- a/src/t8.h +++ b/src/t8.h @@ -295,4 +295,24 @@ t8_sc_array_index_locidx (const sc_array_t *array, const t8_locidx_t index); /* call this at the end of a header file to match T8_EXTERN_C_BEGIN (). */ T8_EXTERN_C_END (); +#ifdef __cplusplus + +/** + * Check two floating point values for equality. + * \tparam TFloating1 Type of the first floating point value. + * \tparam TFloating2 Type of the second floating point value. + * \param [in] floating1 First floating point value to check for equality. + * \param [in] floating2 Second floating point value to check for equality. + * \param [in] tolerance The tolerance. + * \return True if floating values are equal, false otherwise. + */ +template +inline bool +t8_floating_is_equal (TFloating1 floating1, TFloating2 floating2, double tolerance = T8_PRECISION_EPS) noexcept +{ + return std::abs (floating1 - floating2) <= tolerance; +} + +#endif /* !__cplusplus */ + #endif /* !T8_H */ From d00d9b15526192dd65644a5bb38e189ac634b524 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Thu, 8 Jan 2026 15:53:33 +0100 Subject: [PATCH 02/22] major bugfix and restructuring needed to find the bug --- src/t8_cad/t8_cad.cxx | 304 +++++++++--- src/t8_cad/t8_cad.hxx | 171 +++++-- .../t8_cmesh_io/t8_cmesh_readmshfile.cxx | 452 +++++++++--------- .../t8_geometry_cad.cxx | 24 +- 4 files changed, 596 insertions(+), 355 deletions(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index f1a8f9d991..05e4a56ed3 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -32,6 +32,8 @@ #include #include #include +#include +#include t8_cad::t8_cad (std::string fileprefix) { @@ -53,6 +55,7 @@ t8_cad::t8_cad (std::string fileprefix) TopExp::MapShapes (cad_shape, TopAbs_FACE, cad_shape_face_map); TopExp::MapShapesAndUniqueAncestors (cad_shape, TopAbs_VERTEX, TopAbs_EDGE, cad_shape_vertex2edge_map); TopExp::MapShapesAndUniqueAncestors (cad_shape, TopAbs_EDGE, TopAbs_FACE, cad_shape_edge2face_map); + TopExp::MapShapesAndUniqueAncestors (cad_shape, TopAbs_VERTEX, TopAbs_FACE, cad_shape_vertex2face_map); } t8_cad::t8_cad (const TopoDS_Shape cad_shape) @@ -88,26 +91,44 @@ t8_cad::t8_geom_is_plane (const int surface_index) const return surface_adaptor.GetType () == GeomAbs_Plane; } +const TopoDS_Vertex +t8_cad::t8_geom_get_cad_vertex (const int index) const +{ + T8_ASSERT (index <= cad_shape_vertex_map.Size ()); + return TopoDS::Vertex (cad_shape_vertex_map.FindKey (index)); +} + +const TopoDS_Edge +t8_cad::t8_geom_get_cad_edge (const int index) const +{ + T8_ASSERT (index <= cad_shape_edge_map.Size ()); + return TopoDS::Edge (cad_shape_edge_map.FindKey (index)); +} + +const TopoDS_Face +t8_cad::t8_geom_get_cad_face (const int index) const +{ + T8_ASSERT (index <= cad_shape_face_map.Size ()); + return TopoDS::Face (cad_shape_face_map.FindKey (index)); +} + const gp_Pnt t8_cad::t8_geom_get_cad_point (const int index) const { - T8_ASSERT (index <= cad_shape_vertex_map.Size ()); - return BRep_Tool::Pnt (TopoDS::Vertex (cad_shape_vertex_map.FindKey (index))); + return BRep_Tool::Pnt (t8_cad::t8_geom_get_cad_vertex (index)); } const Handle_Geom_Curve t8_cad::t8_geom_get_cad_curve (const int index) const { - T8_ASSERT (index <= cad_shape_edge_map.Size ()); Standard_Real first, last; - return BRep_Tool::Curve (TopoDS::Edge (cad_shape_edge_map.FindKey (index)), first, last); + return BRep_Tool::Curve (t8_cad::t8_geom_get_cad_edge (index), first, last); } const Handle_Geom_Surface t8_cad::t8_geom_get_cad_surface (const int index) const { - T8_ASSERT (index <= cad_shape_face_map.Size ()); - return BRep_Tool::Surface (TopoDS::Face (cad_shape_face_map.FindKey (index))); + return BRep_Tool::Surface (t8_cad::t8_geom_get_cad_face (index)); } const TopTools_IndexedMapOfShape @@ -129,7 +150,7 @@ t8_cad::t8_geom_get_cad_shape_face_map () const } int -t8_cad::t8_geom_get_common_edge (const int vertex1_index, const int vertex2_index) const +t8_cad::t8_geom_get_common_edge_of_vertices (const int vertex1_index, const int vertex2_index) const { const TopTools_ListOfShape collection1 = cad_shape_vertex2edge_map.FindFromIndex (vertex1_index); const TopTools_ListOfShape collection2 = cad_shape_vertex2edge_map.FindFromIndex (vertex2_index); @@ -145,7 +166,7 @@ t8_cad::t8_geom_get_common_edge (const int vertex1_index, const int vertex2_inde } int -t8_cad::t8_geom_get_common_face (const int edge1_index, const int edge2_index) const +t8_cad::t8_geom_get_common_face_of_edges (const int edge1_index, const int edge2_index) const { const TopTools_ListOfShape collection1 = cad_shape_edge2face_map.FindFromIndex (edge1_index); const TopTools_ListOfShape collection2 = cad_shape_edge2face_map.FindFromIndex (edge2_index); @@ -160,6 +181,35 @@ t8_cad::t8_geom_get_common_face (const int edge1_index, const int edge2_index) c return 0; } +int +t8_cad::t8_geom_get_common_face_of_vertex_and_edge (const int vertex_index, const int edge_index) const +{ + const TopTools_ListOfShape edge_collection = cad_shape_edge2face_map.FindFromIndex (edge_index); + for (auto face = edge_collection.begin (); face != edge_collection.end (); ++face) { + const size_t face_index = cad_shape_face_map.FindIndex (*face); + if (t8_geom_is_vertex_on_face (vertex_index, face_index)) { + return face_index; + } + } + return 0; +} + +int +t8_cad::t8_geom_get_common_face_of_vertices (const int vertex1_index, const int vertex2_index) const +{ + const TopTools_ListOfShape collection1 = cad_shape_vertex2face_map.FindFromIndex (vertex1_index); + const TopTools_ListOfShape collection2 = cad_shape_vertex2face_map.FindFromIndex (vertex2_index); + + for (auto face1 = collection1.begin (); face1 != collection1.end (); ++face1) { + for (auto face2 = collection2.begin (); face2 != collection2.end (); ++face2) { + if (face1->IsEqual (*face2)) { + return cad_shape_face_map.FindIndex (*face1); + } + } + } + return 0; +} + int t8_cad::t8_geom_is_vertex_on_edge (const int vertex_index, const int edge_index) const { @@ -187,79 +237,204 @@ t8_cad::t8_geom_is_vertex_on_face (const int vertex_index, const int face_index) return 0; } +bool +t8_cad::t8_geom_vertex_is_on_seam (const int vertex_index, const int face_index) const +{ + const auto face = t8_cad::t8_geom_get_cad_face (face_index); + const TopTools_ListOfShape edge_collection = cad_shape_vertex2edge_map.FindFromIndex (vertex_index); + const ShapeAnalysis_Edge edge_analyzer; + for (auto edge = edge_collection.begin (); edge != edge_collection.end (); ++edge) { + if (edge_analyzer.IsSeam (TopoDS::Edge (*edge), face)) + return true; + } + return false; +} + +bool +t8_cad::t8_geom_edge_is_seam (const int edge_index, const int face_index) const +{ + const auto face = t8_cad::t8_geom_get_cad_face (face_index); + const auto edge = t8_cad::t8_geom_get_cad_edge (edge_index); + const ShapeAnalysis_Edge edge_analyzer; + return edge_analyzer.IsSeam (edge, face); +} + void -t8_cad::t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const int edge_index, double *edge_param) const +t8_cad::t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const int edge_index, double *edge_param, + std::optional reference_edge_param) const { T8_ASSERT (t8_cad::t8_geom_is_vertex_on_edge (vertex_index, edge_index)); TopoDS_Vertex vertex = TopoDS::Vertex (cad_shape_vertex_map.FindKey (vertex_index)); TopoDS_Edge edge = TopoDS::Edge (cad_shape_edge_map.FindKey (edge_index)); - *edge_param = BRep_Tool::Parameter (vertex, edge); + + const bool edge_is_closed = t8_cad::t8_geom_edge_is_closed (edge_index); + if (reference_edge_param && edge_is_closed) { + bool first_point = true; + for (TopExp_Explorer dora (edge, TopAbs_VERTEX); dora.More (); dora.Next ()) { + const TopoDS_Vertex current_vertex = TopoDS::Vertex (dora.Current ()); + if (first_point) { + *edge_param = BRep_Tool::Parameter (current_vertex, edge); + first_point = false; + } + else { + double other_param = BRep_Tool::Parameter (current_vertex, edge); + if (std::abs (other_param - reference_edge_param.value ()) + < std::abs (*edge_param - reference_edge_param.value ())) + *edge_param = other_param; + } + } + } + else { + *edge_param = BRep_Tool::Parameter (vertex, edge); + } } void -t8_cad::t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const int face_index, - double *face_params) const +t8_cad::t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const int face_index, double face_params[2], + std::optional> reference_face_params) const { + /* DISCLAIMER: This function is overly complicated and I do not understand why the simpler versions to not work. + The overly complicated part is only the edge case where the vertex lies on a seam and reference parameters are provided. + The (simpler) plan was as follows: We have the vertex and face. But in the topology of the face, the vertex exists two times, + on the one side of the face and on the other. Bot sides are on the same physical location, because it is a seam. Thats + why both vertices have the same coordinates. But they do not have the same parameters on the surface. The simple plan + is to use a TopExp_Explorer to iterate over all vertices of the face and to look at all vertices with the right coords. + From all those vertices we choose the one, which has the closest parameters to our reference parameters. + But for whatever reason the surface in my test case only has vertices with one parameterset. Not with the other set. + And thats why this algorithm does not work. What you see here is the more complicated workaround: + */ + T8_ASSERT (t8_cad::t8_geom_is_vertex_on_face (vertex_index, face_index)); - gp_Pnt2d uv; + std::optional uv = std::nullopt; /** Final parameters on the surface */ TopoDS_Vertex vertex = TopoDS::Vertex (cad_shape_vertex_map.FindKey (vertex_index)); TopoDS_Face face = TopoDS::Face (cad_shape_face_map.FindKey (face_index)); - uv = BRep_Tool::Parameters (vertex, face); - face_params[0] = uv.X (); - face_params[1] = uv.Y (); + + /* If the surface is not closed or the user did not provide any reference params we just query the parameters. */ + if (!t8_cad::t8_geom_surface_is_closed (face_index) || !reference_face_params) { + uv.emplace (BRep_Tool::Parameters (vertex, face)); + } + + /* If the vertex is not on the seam we can also just query the parameters. */ + else { + const bool is_on_seam = t8_cad::t8_geom_vertex_is_on_seam (vertex_index, face_index); + if (!is_on_seam) { + uv.emplace (BRep_Tool::Parameters (vertex, face)); + } + /* Now the overcomplicated edge case: The user provided reference params and the vertex is on the seam of the surface. */ + else { + /* Convert reference parameters to a OCCT point */ + gp_Pnt2d reference_point (reference_face_params.value ()[0], reference_face_params.value ()[1]); + + /* The vertex is on a seam and so we have to get all vertices and check which one is closer. + But since the vertices we are searching for are not in the topology of the face we have to take a detour. + We iterate over all edges of the surface to find the seam edge. Yes, the seam EDGE. Not the seam EDGES. + I would think that there are two seam edges. On both sides of the face there is one. And for a cylindrical + surface the TopExp_Explorer gives us four edges. Two circles and two lines. I would say, that both lines are + seams, because they are the same line. But OCCT only marks one of the two as seam. We now try to find this seam. + */ + + /** List of all edges connected to the vertex */ + const TopTools_ListOfShape vertex_edges = cad_shape_vertex2edge_map.FindFromIndex (vertex_index); + const ShapeAnalysis_Edge edge_analyzer; + std::optional seam = std::nullopt; + + for (TopExp_Explorer dora (face, TopAbs_EDGE); dora.More (); dora.Next ()) { + const TopoDS_Edge current_edge = TopoDS::Edge (dora.Current ()); + /* The edge has to be a seam and has to be connected to our vertex. */ + const bool current_edge_connected_to_vertex = vertex_edges.Contains (current_edge); + const bool current_edge_is_seam = edge_analyzer.IsSeam (current_edge, face); + if (current_edge_connected_to_vertex && current_edge_is_seam) { + seam.emplace (current_edge); + break; + } + } + /* Hopefully we found a seam. If not something has gone catastrophically wrong or I overlooked something. + Probably the latter. */ + SC_CHECK_ABORTF (seam.has_value (), "Error: Could not find seam of periodic surface."); + + /* We now have the seam and can iterate over all edges of the face again. Every edge matching our seam is also a seam (in my eyes). + In the case of our cylinder this would be the two lines. */ + for (TopExp_Explorer dora (face, TopAbs_EDGE); dora.More (); dora.Next ()) { + const TopoDS_Edge current_edge = TopoDS::Edge (dora.Current ()); + if (seam->IsSame (current_edge)) { + /* If the edge matches our seam we can query the parameter of the vertex on the edge on the face. + Why also on the face and not only on the edge? I really don't know but I had all information to call this function and it seems to work. */ + Standard_Real first, last; + const double curve_param = BRep_Tool::Parameter (vertex, current_edge, face); + const Handle_Geom2d_Curve pcurve_on_surface = BRep_Tool::CurveOnSurface (current_edge, face, first, last); + /* We can now insert the queried parameter on the curve into the pcurve (the curve on the surface). + If this is the first point we found uv is empty and we fill it. Otherwise, we define another point to + stor our parameters and we save the parameters in iv which are closer to the reference parameters. */ + if (!uv) { + uv.emplace (); + pcurve_on_surface->D0 (curve_param, *uv); + } + /* If it is the second one we check which point is close. */ + else { + gp_Pnt2d other_point; + pcurve_on_surface->D0 (curve_param, other_point); + if (other_point.Distance (reference_point) < uv->Distance (reference_point)) + uv = other_point; + } + } + } + } + } + + /* Lastly, we put the parameters into our output array. */ + face_params[0] = uv->X (); + face_params[1] = uv->Y (); } void -t8_cad::t8_geom_edge_parameter_to_face_parameters (const int edge_index, const int face_index, const int num_face_nodes, - const double edge_param, const double *surface_params, - double *face_params) const +t8_cad::t8_geom_edge_parameter_to_face_parameters ( + const int edge_index, const int face_index, const double edge_param, double face_params_out[2], + std::optional> reference_face_params) const { T8_ASSERT (t8_cad::t8_geom_is_edge_on_face (edge_index, face_index)); Standard_Real first, last; gp_Pnt2d uv; - TopoDS_Edge edge = TopoDS::Edge (cad_shape_edge_map.FindKey (edge_index)); TopoDS_Face face = TopoDS::Face (cad_shape_face_map.FindKey (face_index)); - Handle_Geom2d_Curve curve_on_surface = BRep_Tool::CurveOnSurface (edge, face, first, last); - Handle_Geom_Surface surface = BRep_Tool::Surface (face); - curve_on_surface->D0 (edge_param, uv); - face_params[0] = uv.X (); - face_params[1] = uv.Y (); - - /* Check for right conversion of edge to surface parameter and correct if needed */ - /* Checking u parameter */ - if (surface_params != NULL) { - double parametric_bounds[4]; - surface->Bounds (parametric_bounds[0], parametric_bounds[1], parametric_bounds[2], parametric_bounds[3]); - if (surface->IsUClosed ()) { - for (int i_face_node = 0; i_face_node < num_face_nodes; ++i_face_node) { - if (surface_params[i_face_node * 2] == parametric_bounds[0]) { - if (face_params[0] == parametric_bounds[1]) { - face_params[0] = parametric_bounds[0]; - } - } - else if (surface_params[i_face_node * 2] == parametric_bounds[1]) { - if (face_params[0] == parametric_bounds[0]) { - face_params[0] = parametric_bounds[1]; - } - } - } - } - /* Checking v parameter */ - if (surface->IsVClosed ()) { - for (int i_face_node = 0; i_face_node < num_face_nodes; ++i_face_node) { - if (surface_params[i_face_node * 2 + 1] == parametric_bounds[0]) { - if (face_params[1] == parametric_bounds[1]) { - face_params[1] = parametric_bounds[0]; - } + TopoDS_Edge edge = TopoDS::Edge (cad_shape_edge_map.FindKey (edge_index)); + + const bool is_seam = t8_cad::t8_geom_edge_is_seam (edge_index, face_index); + if (is_seam && reference_face_params) { + /* Convert reference parameters to OCCT point */ + gp_Pnt2d reference_point (reference_face_params.value ()[0], reference_face_params.value ()[1]); + + /* If the edge is a seam we have to get both points and check which one is closer. + We iterate over all edges of the face and check if the edges are the same as the + input edge. Then we check if the converted parameters are closer to the + already computed parameters. */ + bool first_point = true; + for (TopExp_Explorer dora (face, TopAbs_EDGE); dora.More (); dora.Next ()) { + const TopoDS_Edge current_edge = TopoDS::Edge (dora.Current ()); + /* Check if edge is one of the seams. */ + if (edge.IsSame (current_edge)) { + Handle_Geom2d_Curve curve_on_surface = BRep_Tool::CurveOnSurface (edge, face, first, last); + /* If it is the first seam we compute the parameters. */ + if (first_point) { + curve_on_surface->D0 (edge_param, uv); + first_point = false; } - else if (surface_params[i_face_node * 2 + 1] == parametric_bounds[1]) { - if (face_params[1] == parametric_bounds[0]) { - face_params[1] = parametric_bounds[1]; - } + /* If it is the second one we check which point is close. */ + else { + gp_Pnt2d other_point; + curve_on_surface->D0 (edge_param, other_point); + if (other_point.Distance (reference_point) < uv.Distance (reference_point)) + uv = other_point; } } } } + else { + /* If the edge is not a seam or if no reference parameters were provided we just + use the normal edge. */ + Handle_Geom2d_Curve curve_on_surface = BRep_Tool::CurveOnSurface (edge, face, first, last); + curve_on_surface->D0 (edge_param, uv); + } + face_params_out[0] = uv.X (); + face_params_out[1] = uv.Y (); } void @@ -277,18 +452,25 @@ t8_cad::t8_geom_get_edge_parametric_bounds (const int edge_index, double *bounds bounds[1] = cad_edge->LastParameter (); } -int -t8_cad::t8_geom_is_edge_closed (int edge_index) const +bool +t8_cad::t8_geom_edge_is_closed (int edge_index) const { const Handle_Geom_Curve cad_edge = t8_geom_get_cad_curve (edge_index); return cad_edge->IsClosed (); } -int -t8_cad::t8_geom_is_surface_closed (int geometry_index, int parameter) const +bool +t8_cad::t8_geom_surface_is_closed (int geometry_index) const +{ + const Handle_Geom_Surface cad_surface = t8_geom_get_cad_surface (geometry_index); + return cad_surface->IsUClosed () || cad_surface->IsVClosed (); +} + +bool +t8_cad::t8_geom_surface_is_closed (int geometry_index, int direction) const { const Handle_Geom_Surface cad_surface = t8_geom_get_cad_surface (geometry_index); - switch (parameter) { + switch (direction) { case 0: return cad_surface->IsUClosed (); break; diff --git a/src/t8_cad/t8_cad.hxx b/src/t8_cad/t8_cad.hxx index 25ebd82f89..e17fa18fa9 100644 --- a/src/t8_cad/t8_cad.hxx +++ b/src/t8_cad/t8_cad.hxx @@ -20,7 +20,7 @@ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -/** \file t8_geometry_cad.hxx +/** \file t8_cad.hxx * This file implements the t8_cad class. It manages OpenCASCADE shapes and implements * helper functions for working with the shapes. */ @@ -32,6 +32,11 @@ #include #include #include +#include +#include +#include +#include +#include /** * This class manages OpenCASCADE shapes and implements helper functions for working with the shapes. @@ -80,21 +85,42 @@ class t8_cad { int t8_geom_is_plane (const int surface_index) const; - /** Get an cad point from the cad_shape. + /** Get a cad vertex from the cad_shape. + * \param [in] index The index of the vertex in the cad_shape. + * \return The cad vertex. + */ + const TopoDS_Vertex + t8_geom_get_cad_vertex (const int index) const; + + /** Get a cad edge from the cad_shape. + * \param [in] index The index of the edge in the cad_shape. + * \return The cad edge. + */ + const TopoDS_Edge + t8_geom_get_cad_edge (const int index) const; + + /** Get a cad face from the cad_shape. + * \param [in] index The index of the face in the cad_shape. + * \return The cad face. + */ + const TopoDS_Face + t8_geom_get_cad_face (const int index) const; + + /** Get a cad point from the cad_shape. * \param [in] index The index of the point in the cad_shape. * \return The cad point. */ const gp_Pnt t8_geom_get_cad_point (const int index) const; - /** Get an cad curve from the cad_shape. + /** Get a cad curve from the cad_shape. * \param [in] index The index of the curve in the cad_shape. * \return The cad curve. */ const Handle_Geom_Curve t8_geom_get_cad_curve (const int index) const; - /** Get an cad surface from the cad_shape. + /** Get a cad surface from the cad_shape. * \param [in] index The index of the surface in the cad_shape. * \return The cad surface. */ @@ -125,7 +151,7 @@ class t8_cad { * \return Index of the shared edge. 0 if there is no shared edge. */ int - t8_geom_get_common_edge (const int vertex1_index, const int vertex2_index) const; + t8_geom_get_common_edge_of_vertices (const int vertex1_index, const int vertex2_index) const; /** Check if two cad edges share a common cad face. * \param [in] edge1_index The index of the first cad edge. @@ -133,7 +159,23 @@ class t8_cad { * \return Index of the shared face. 0 if there is no shared face. */ int - t8_geom_get_common_face (const int edge1_index, const int edge2_index) const; + t8_geom_get_common_face_of_edges (const int edge1_index, const int edge2_index) const; + + /** Check if a cad vertex and cad edge share a common face. + * \param [in] vertex_index The index of the cad vertex. + * \param [in] edge_index The index of the cad edge. + * \return Index of the shared face. 0 if there is no shared face. + */ + int + t8_geom_get_common_face_of_vertex_and_edge (const int vertex_index, const int edge_index) const; + + /** Check if two cad vertices share a common cad face. + * \param [in] vertex1_index The index of the first cad edge. + * \param [in] vertex2_index The index of the second cad edge. + * \return Index of the shared face. 0 if there is no shared face. + */ + int + t8_geom_get_common_face_of_vertices (const int vertex1_index, const int vertex2_index) const; /** Check if a cad vertex lies on an cad edge. * \param [in] vertex_index The index of the cad vertex. @@ -159,42 +201,76 @@ class t8_cad { int t8_geom_is_vertex_on_face (const int vertex_index, const int face_index) const; + /** Returns true if \a vertex_index is on a seam of \a face_index. + * A seam is an edge which connects a surface to itself. + * + * \param [in] vertex_index + * \param [in] face_index + * \return true if the edge is a seam. false otherwise. + */ + bool + t8_geom_vertex_is_on_seam (const int vertex_index, const int face_index) const; + + /** Returns true if \a edge_index is a seam of \a face_index. + * A seam is an edge which connects a surface to itself. + * + * \param [in] edge_index + * \param [in] face_index + * \return true if the edge is a seam. false otherwise. + */ + bool + t8_geom_edge_is_seam (const int edge_index, const int face_index) const; + /** Retrieves the parameter of an cad vertex on an cad edge. - * The vertex has to lie on the edge. + * The vertex has to lie on the edge. + * \warning If the edge is closed in any direction and the vertex is the closing bound, + * it is random which side of the closed edge the parameter is from. The parameter + * is correct, but for curved elements it has to be checked if the parameter has to be + * converted onto the other bound. * \param [in] vertex_index The index of the cad vertex. * \param [in] edge_index The index of the cad edge. * \param [out] edge_param The parameter of the vertex on the edge. + * \param [in] reference_face_params Reference parameters on the surface. */ void - t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const int edge_index, double *edge_param) const; + t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const int edge_index, double *edge_param, + std::optional reference_edge_param = std::nullopt) const; /** Retrieves the parameters of an cad vertex on a cad face. - * The vertex has to lie on the face. - * \param [in] vertex_index The index of the cad vertex. - * \param [in] face_index The index of the cad face. - * \param [out] face_params The parameters of the vertex on the face. + * The vertex has to lie on the face. + * If the vertex is in a seam, the vertex closer to the \a reference_face_params is chosen. + * \warning If the face is closed in any direction and the vertex is on the seam and no reference + * parameters are provided it is random which side of the closed face the parameters are from. + * The parameters are correct, but for curved elements it has to be checked if the parameters have to be + * converted onto the other bound. + * \param [in] vertex_index The index of the cad vertex. + * \param [in] face_index The index of the cad face. + * \param [out] face_params The parameters of the vertex on the face. + * \param [in] reference_face_params Reference parameters on the surface. */ void - t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const int face_index, double *face_params) const; + t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const int face_index, double face_params[2], + std::optional> reference_face_params + = std::nullopt) const; - /** Converts the parameters of an cad edge to the corresponding parameters on an cad face. + /** Converts the parameters of a cad edge to the corresponding parameters on a cad face. * The edge has to lie on the face. - * For the conversion of edge parameters of mesh elements to topological face parameters of a closed surface, it is additionally - * checked, whether the conversion was correct, to prevent disorted elements. - * \param [in] edge_index The index of the cad edge, which parameters should be converted to face parameters. - * \param [in] face_index The index of the cad face, on to which the edge parameters should be converted. - * \param [in] num_face_nodes The number of the face nodes of the evaluated element. Only needed for closed surface check, otherwise NULL. - * \param [in] edge_param The parameter on the edge. - * \param [in] surface_params The parameters of the surface nodes. - * When provided, there are additional checks for closed geometries. - * If there are no surface parameter - * to pass in to the function, you can pass NULL. - * \param [out] face_params The corresponding parameters on the face. + * If the edge is a seam, the edge closer to the \a reference_face_params is chosen. + * \warning If the face is closed in any direction and the edge is the seam and no reference + * parameters are provided it is random which side of the closed face the parameters are from. + * The parameters are correct, but for curved elements it has to be checked if the parameters have to be + * converted onto the other side of the seam. + * \param [in] edge_index The index of the cad edge, which parameters should be converted to face parameters. + * \param [in] face_index The index of the cad face, on to which the edge parameters should be converted. + * \param [in] edge_param The parameter on the edge. + * \param [out] face_params_out The corresponding parameters on the face. + * \param [in] reference_face_params Reference parameters on the surface. */ void - t8_geom_edge_parameter_to_face_parameters (const int edge_index, const int face_index, const int num_face_nodes, - const double edge_param, const double *surface_params, - double *face_params) const; + t8_geom_edge_parameter_to_face_parameters (const int edge_index, const int face_index, const double edge_param, + double face_params_out[2], + std::optional> reference_face_params + = std::nullopt) const; /** Finds the parametric bounds of an cad face. * \param [in] surface_index The index of the cad face. @@ -210,31 +286,40 @@ class t8_cad { void t8_geom_get_edge_parametric_bounds (const int edge_index, double *bounds) const; - /** Checks if an edge is closed in its U parameter. + /** Checks if an edge is closed. * \param [in] edge_index The index of the closed edge. - * \return 1 if edge is closed in U. 0 if edge is not closed in U. + * \return true if edge is closed */ - int - t8_geom_is_edge_closed (int edge_index) const; + bool + t8_geom_edge_is_closed (int edge_index) const; - /** Checks if a surface is closed in its U parameter or V parameter. + /** Checks if a surface is closed in any direction. * \param [in] geometry_index The index of the closed geometry. - * \param [in] parameter The parameter, which should be check for closeness. - * 0 stands for the U parameter and 1 for the V parameter. - * \return 1 if geometry is closed in U. 0 if geometry is not closed in U. + * \return true if geometry is closed in any direction. */ - int - t8_geom_is_surface_closed (int geometry_index, int parameter) const; + bool + t8_geom_surface_is_closed (int geometry_index) const; + + /** Checks if a surface is closed in its U direction or V direction. + * \param [in] geometry_index The index of the closed geometry. + * \param [in] direction The direction, which should be check for closeness. + * 0 stands for the U direction and 1 for the V direction. + * \return true if geometry is closed in the given direction. + */ + bool + t8_geom_surface_is_closed (int geometry_index, int direction) const; private: TopoDS_Shape cad_shape; /**< cad geometry */ - TopTools_IndexedMapOfShape cad_shape_vertex_map; /**< Map of all TopoDS_Vertex in shape. */ - TopTools_IndexedMapOfShape cad_shape_edge_map; /**< Map of all TopoDS_Edge in shape. */ - TopTools_IndexedMapOfShape cad_shape_face_map; /**< Map of all TopoDS_Face in shape. */ + TopTools_IndexedMapOfShape cad_shape_vertex_map; /**< Map of all TopoDS_Vertex. */ + TopTools_IndexedMapOfShape cad_shape_edge_map; /**< Map of all TopoDS_Edge. */ + TopTools_IndexedMapOfShape cad_shape_face_map; /**< Map of all TopoDS_Face. */ + TopTools_IndexedDataMapOfShapeListOfShape + cad_shape_vertex2edge_map; /**< Maps all TopoDS_Vertex to all its connected TopoDS_Edge */ TopTools_IndexedDataMapOfShapeListOfShape - cad_shape_vertex2edge_map; /**< Maps all TopoDS_Vertex of shape to all its connected TopoDS_Edge */ + cad_shape_edge2face_map; /**< Maps all TopoDS_Edge to all its connected TopoDS_Face */ TopTools_IndexedDataMapOfShapeListOfShape - cad_shape_edge2face_map; /**< Maps all TopoDS_Edge of shape to all its connected TopoDS_Face */ + cad_shape_vertex2face_map; /**< Maps all TopoDS_Vertex to all its connected TopoDS_Face */ }; #endif /* !T8_CAD_HXX */ diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 690bb85a11..a4b4f088b9 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -33,6 +33,8 @@ #include #include #include +#include +#include #ifdef _WIN32 #include "t8_windows.h" @@ -641,7 +643,7 @@ t8_cmesh_msh_file_2_read_eles (t8_cmesh_t cmesh, FILE *fp, const t8_msh_node_tab t8_locidx_t num_trees, tree_loop; t8_gloidx_t tree_count; t8_eclass_t eclass; - t8_msh_file_node Node; + t8_msh_file_node node; long lnum_trees; int retval; int ele_type, num_tags; @@ -751,8 +753,8 @@ t8_cmesh_msh_file_2_read_eles (t8_cmesh_t cmesh, FILE *fp, const t8_msh_node_tab } /* Get node from the hashtable */ - Node.index = node_indices[t8_vertex_num]; - const auto found_node = vertices.find (Node); + node.index = node_indices[t8_vertex_num]; + const auto found_node = vertices.find (node); if (found_node == vertices.end ()) { t8_global_errorf ("Could not find Node %li.\n", node_indices[t8_vertex_num]); free (line); @@ -786,110 +788,6 @@ t8_cmesh_msh_file_2_read_eles (t8_cmesh_t cmesh, FILE *fp, const t8_msh_node_tab return std::make_optional (vertex_indices); } -#if T8_ENABLE_OCC -/** Corrects the parameters on closed geometries to prevent disorted elements. - * \param [in] geometry_dim The dimension of the geometry. - * 1 for edges, 2 for surfaces. - * \param [in] geometry_index The index of the geometry. - * \param [in] num_face_nodes The number of the nodes of the surface. - * NULL if the geometry is an edge. - * \param [in] geometry_cad The cad_geometry. - * \param [in,out] parameters The parameters to be corrected. - */ -static void -t8_cmesh_correct_parameters_on_closed_geometry (const int geometry_dim, const int geometry_index, - const int num_face_nodes, const t8_geometry_cad *geometry_cad, - double *parameters) -{ - switch (geometry_dim) { - /* Check for closed U parameter in case of an edge. */ - case 1: - /* Only correct the U parameter if the edge is closed. */ - if (geometry_cad->get_cad_manager ()->t8_geom_is_edge_closed (geometry_index)) { - /* Get the parametric bounds of the closed geometry - * edge -> [Umin, Umax] - */ - double parametric_bounds[2]; - /* Get the parametric edge bounds. */ - geometry_cad->get_cad_manager ()->t8_geom_get_edge_parametric_bounds (geometry_index, parametric_bounds); - /* Check the upper an the lower parametric bound. */ - for (int bound = 0; bound < 2; ++bound) { - /* Iterate over both nodes of the edge. */ - for (int i_nodes = 0; i_nodes < 2; ++i_nodes) { - /* Check if one of the U parameters lies on one of the parametric bounds. */ - if (std::abs (parameters[i_nodes] - parametric_bounds[bound]) <= T8_PRECISION_EPS) { - /* Check the U parameter of the other node ((i_node + 1) % 2) to find out - * to which parametric bound the tree is closer. - */ - if (std::abs (parameters[(i_nodes + 1) % 2] - parametric_bounds[bound]) > T8_PRECISION_EPS) { - /* Now check if the difference of the parameters of both nodes are bigger than the half parametric range. - * In this case, the parameter at i_nodes has to be changed to the other parametric bound ((bound + 1) % 2). - */ - if (std::abs (parameters[(i_nodes + 1) % 2] - parameters[i_nodes]) - > ((parametric_bounds[1] - parametric_bounds[0]) / 2)) { - /* Switch to the other parametric bound. */ - parameters[i_nodes] = parametric_bounds[(bound + 1) % 2]; - break; - } - } - } - } - } - } - break; - - /* Check for closed U parameter and closed V parameter in case of a surface. */ - case 2: - /* Iterate over both parameters. 0 stands for the U parameter an 1 for the V parameter. */ - for (int param_dim = 0; param_dim < 2; ++param_dim) { - /* Only correct the surface parameters if they are closed */ - if (geometry_cad->get_cad_manager ()->t8_geom_is_surface_closed (geometry_index, param_dim)) { - /* Get the parametric bounds of the closed geometry - * surface -> [Umin, Umax, Vmin, Vmax] - */ - double parametric_bounds[4]; - geometry_cad->get_cad_manager ()->t8_geom_get_face_parametric_bounds (geometry_index, parametric_bounds); - /* Check the upper an the lower parametric bound. */ - for (int bound = 0; bound < 2; ++bound) { - /* Iterate over every corner node of the tree. */ - for (int i_nodes = 0; i_nodes < num_face_nodes; ++i_nodes) { - /* Check if one of the U parameters lies on one of the parametric bounds. */ - if (std::abs (parameters[2 * i_nodes + param_dim] - parametric_bounds[bound + 2 * param_dim]) - <= T8_PRECISION_EPS) { - /* Iterate over every corner node of the tree again. */ - for (int j_nodes = 0; j_nodes < num_face_nodes; ++j_nodes) { - /* Search for a U parameter that is non of the parametric bounds. To check - * whether the tree is closer to the lower or the upper parametric bound. - */ - if (std::abs (parameters[2 * j_nodes + param_dim] - parametric_bounds[bound + 2 * param_dim]) - > T8_PRECISION_EPS - && std::abs (parameters[2 * j_nodes + param_dim] - - parametric_bounds[((bound + 1) % 2) + 2 * param_dim]) - > T8_PRECISION_EPS) { - /* Now check if the difference of the parameters of both nodes are bigger than the half parametric range. - * In this case, the parameter at i_nodes has to be changed to the other parametric bound ((bound + 1) % 2). - */ - if (std::abs (parameters[2 * j_nodes + param_dim] - parameters[2 * i_nodes + param_dim]) - > ((parametric_bounds[1 + 2 * param_dim] - parametric_bounds[0 + 2 * param_dim]) / 2)) { - /* Switch to the other parametric bound. */ - parameters[2 * i_nodes + param_dim] = parametric_bounds[((bound + 1) % 2) + 2 * param_dim]; - break; - } - } - } - } - } - } - } - } - break; - default: - SC_ABORT_NOT_REACHED (); - break; - } -} -#endif /* T8_ENABLE_OCC */ - /** * This function stores the entity dimensions, tags, and parameters of each node in a given * tree in arrays. These arrays are then passed to the tree as attributes in cmesh. @@ -948,11 +846,10 @@ t8_store_element_node_data (t8_cmesh_t cmesh, t8_gloidx_t tree_count, * */ static bool -t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t8_gloidx_t tree_count, - const t8_geometry_c *cad_geometry_base, const t8_geometry_c *linear_geometry_base, - std::array tree_nodes, - std::array face_nodes, - std::array edge_nodes) +t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass, const int dim, + const t8_gloidx_t tree_count, const t8_geometry_c *cad_geometry_base, + const t8_geometry_c *linear_geometry_base, + const std::array &tree_nodes) { /* Calculate the parametric geometries of the tree */ T8_ASSERT (cad_geometry_base->t8_geom_get_type () == T8_GEOMETRY_TYPE_CAD); @@ -987,27 +884,34 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t default: SC_ABORTF ("Invalid dimension of tree. Dimension: %i\n", dim); } + for (int i_tree_faces = 0; i_tree_faces < num_faces; ++i_tree_faces) { const int face_eclass = dim == 2 ? eclass : t8_eclass_face_types[eclass][i_tree_faces]; const int num_face_nodes = t8_eclass_num_vertices[face_eclass]; const int num_face_edges = t8_eclass_num_faces[face_eclass]; + std::vector face_nodes; + face_nodes.reserve (T8_ECLASS_MAX_CORNERS_2D); + + /*----------------------------------------- Start of face-surface linkage -----------------------------------------*/ /* Save each node of face separately. Face nodes of 2D elements are also tree nodes. - * Face nodes of 3D elements need to be translated to tree nodes. */ + * Face nodes of 3D elements need to be translated to tree nodes. */ for (int i_face_node = 0; i_face_node < num_face_nodes; ++i_face_node) { if (dim == 2) { - face_nodes[i_face_node] = tree_nodes[i_face_node]; + face_nodes.push_back (tree_nodes[i_face_node]); + T8_ASSERT (num_faces == 1); } else { - face_nodes[i_face_node] = tree_nodes[t8_face_vertex_to_tree_vertex[eclass][i_tree_faces][i_face_node]]; + face_nodes.push_back (tree_nodes[t8_face_vertex_to_tree_vertex[eclass][i_tree_faces][i_face_node]]); } } + T8_ASSERT (face_nodes.size () == (size_t) num_face_nodes); - /* A face can only be linked to an cad surface if all nodes of the face are parametric or on a vertex - * (gmsh labels nodes on vertices as not parametric) */ + /* A face can only be linked to a cad surface if all nodes of the face are parametric or on a vertex + * (gmsh labels nodes on vertices as not parametric) */ int all_parametric = 1; - for (int i_face_nodes = 0; i_face_nodes < num_face_nodes; ++i_face_nodes) { - if (!face_nodes[i_face_nodes].parametric && face_nodes[i_face_nodes].entity_dim != 0) { + for (auto &face_node : face_nodes) { + if (!face_node.parametric && face_node.entity_dim != 0) { all_parametric = 0; break; } @@ -1017,47 +921,50 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t continue; } /* Now we can check if the face is connected to a surface */ - int surface_index = 0; + int surface_index = 0; /** Index of the cad surface we link the element face to. */ /* If one node is already on a surface we can check if the rest lies also on the surface. */ - for (int i_face_nodes = 0; i_face_nodes < num_face_nodes; ++i_face_nodes) { - if (face_nodes[i_face_nodes].entity_dim == 2) { - surface_index = face_nodes[i_face_nodes].entity_tag; + for (auto &face_node : face_nodes) { + if (face_node.entity_dim == 2) { + surface_index = face_node.entity_tag; break; } } /* If not we can take two curves and look if they share a surface and then use this surface */ if (!surface_index) { /* To do this we can look if there are two curves, otherwise we have to check which vertices - * share the same curve. */ + * share the same curve. */ int edge1_index = 0; int edge2_index = 0; /* We search for 2 different curves */ - for (int i_face_nodes = 0; i_face_nodes < num_face_nodes; ++i_face_nodes) { - if (face_nodes[i_face_nodes].entity_dim == 1) { + for (auto &face_node : face_nodes) { + if (face_node.entity_dim == 1) { if (edge1_index == 0) { - edge1_index = face_nodes[i_face_nodes].entity_tag; + /* Fill the first edge index. */ + edge1_index = face_node.entity_tag; } - else if (face_nodes[i_face_nodes].entity_tag != edge1_index) { - edge2_index = face_nodes[i_face_nodes].entity_tag; + else if (face_node.entity_tag != edge1_index) { + /* Fill the second index if the first one is already filled. */ + edge2_index = face_node.entity_tag; break; } } } - /* If there are less than 2 curves we can look at the vertices and check, - * if two of them are on the same curve */ + /* We now have either 0 curves, 1 curve or 2 curves. But we need at least 2 curves to continue. + * So if we have less than 2 curves we can check if two vertices share the same curve + * and use these to fill our two curves. */ if (edge2_index == 0) { /* For each edge of face */ for (int i_face_edges = 0; i_face_edges < num_face_edges; ++i_face_edges) { /* Save nodes separately */ - const int node1_number = t8_face_vertex_to_tree_vertex[face_eclass][i_tree_faces][0]; + const int node1_number = t8_face_vertex_to_tree_vertex[face_eclass][i_face_edges][0]; const t8_msh_file_node node1 = face_nodes[node1_number]; - const int node2_number = t8_face_vertex_to_tree_vertex[face_eclass][i_tree_faces][1]; + const int node2_number = t8_face_vertex_to_tree_vertex[face_eclass][i_face_edges][1]; const t8_msh_file_node node2 = face_nodes[node2_number]; /* If both nodes are on a vertex we look if both vertices share an edge */ if (node1.entity_dim == 0 && node2.entity_dim == 0) { - int common_edge - = cad_geometry->get_cad_manager ()->t8_geom_get_common_edge (node1.entity_tag, node2.entity_tag); + int common_edge = cad_geometry->get_cad_manager ()->t8_geom_get_common_edge_of_vertices (node1.entity_tag, + node2.entity_tag); if (common_edge > 0) { if (edge1_index == 0) { edge1_index = common_edge; @@ -1070,8 +977,9 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t } } } + /* If the face is linked to a surface we should now have two curves and we can check if both share a surface. */ if (edge2_index > 0) { - surface_index = cad_geometry->get_cad_manager ()->t8_geom_get_common_face (edge1_index, edge2_index); + surface_index = cad_geometry->get_cad_manager ()->t8_geom_get_common_face_of_edges (edge1_index, edge2_index); } else { continue; @@ -1080,54 +988,91 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t /* Now we can check if every node lies on the surface and retrieve its parameters */ if (surface_index) { int all_nodes_on_surface = 1; - for (int i_face_nodes = 0; i_face_nodes < num_face_nodes; ++i_face_nodes) { - /* We check if the node is on the right surface */ - if (face_nodes[i_face_nodes].entity_dim == 2) { - /* Check if node is on the right surface */ - if (face_nodes[i_face_nodes].entity_tag != surface_index) { - all_nodes_on_surface = 0; - break; + /* If the surface is closed in any direction we may need some reference parameters on the surface. + These reference parameters can be used to find the right parameters if we need to convert + seam edge parameters to the face parameters. */ + std::optional> reference_parameters_on_surface (std::nullopt); + /* We iterate twice over all face nodes. The first time we look, if all nodes are on the right surface and + search for the reference parameters. In the second loop we convert all coordinates which are not directly on that surface. */ + for (auto face_node = face_nodes.begin (); face_node != face_nodes.end () && all_nodes_on_surface; ++face_node) { + switch (face_node->entity_dim) { + case 2: + /* If the node is on the surface, we can use its parameters as reference parameters. */ + if (face_node->entity_tag == surface_index) { + if (!reference_parameters_on_surface) + reference_parameters_on_surface = face_node->parameters; } - } - else { - /* If it is on another geometry we retrieve its parameters */ - if (face_nodes[i_face_nodes].entity_dim == 0) { - if (cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_face (face_nodes[i_face_nodes].entity_tag, - surface_index)) { - cad_geometry->get_cad_manager ()->t8_geom_get_parameters_of_vertex_on_face ( - face_nodes[i_face_nodes].entity_tag, surface_index, face_nodes[i_face_nodes].parameters.data ()); - face_nodes[i_face_nodes].entity_dim = 2; - } - else { - all_nodes_on_surface = 0; - break; - } + /* If not we can stop the face linkage. */ + else { + all_nodes_on_surface = 0; } - if (face_nodes[i_face_nodes].entity_dim == 1) { - if (cad_geometry->get_cad_manager ()->t8_geom_is_edge_on_face (face_nodes[i_face_nodes].entity_tag, - surface_index)) { + break; + case 1: + /* If the edge is on the face we can check if it is a seam and if not we can use its parameters as reference. */ + if (cad_geometry->get_cad_manager ()->t8_geom_is_edge_on_face (face_node->entity_tag, surface_index)) { + if (!reference_parameters_on_surface + && !cad_geometry->get_cad_manager ()->t8_geom_edge_is_seam (face_node->entity_tag, surface_index)) { + reference_parameters_on_surface.emplace (); cad_geometry->get_cad_manager ()->t8_geom_edge_parameter_to_face_parameters ( - face_nodes[i_face_nodes].entity_tag, surface_index, num_face_nodes, - face_nodes[i_face_nodes].parameters[0], NULL, face_nodes[i_face_nodes].parameters.data ()); - face_nodes[i_face_nodes].entity_dim = 2; - } - else { - all_nodes_on_surface = 0; - break; + face_node->entity_tag, surface_index, face_node->parameters[0], + reference_parameters_on_surface.value ().data ()); } } + /* If the edge is not on the face, we can stop the face linkage. */ + else { + all_nodes_on_surface = 0; + } + break; + case 0: + /* If the vertex is not on the face, we can stop the face linkage. + We do not have to check for reference parameters since a closed face linked to four vertices makes no sense. */ + if (!cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_face (face_node->entity_tag, surface_index)) + all_nodes_on_surface = 0; + break; + default: + SC_ABORT_NOT_REACHED (); } } /* Abort if not all nodes are on the surface or if the surface is a plane */ if (!all_nodes_on_surface || cad_geometry->get_cad_manager ()->t8_geom_is_plane (surface_index)) { continue; } + + if (!reference_parameters_on_surface.has_value ()) { + t8_global_errorf ("Error during mesh-cad recombination: Reference parameters on surface not found.\n"); + return 0; + } + + /* We iterate again to convert all parameters. */ + for (auto &face_node : face_nodes) { + switch (face_node.entity_dim) { + case 2: + /* Nothing to do. */ + break; + case 1: + cad_geometry->get_cad_manager ()->t8_geom_edge_parameter_to_face_parameters ( + face_node.entity_tag, surface_index, face_node.parameters[0], face_node.parameters.data (), + reference_parameters_on_surface); + face_node.entity_dim = 2; + face_node.entity_tag = surface_index; + break; + case 0: + cad_geometry->get_cad_manager ()->t8_geom_get_parameters_of_vertex_on_face ( + face_node.entity_tag, surface_index, face_node.parameters.data (), reference_parameters_on_surface); + face_node.entity_dim = 2; + face_node.entity_tag = surface_index; + break; + default: + SC_ABORT_NOT_REACHED (); + } + } + /* If we have found a surface we link it to the face */ face_geometries[i_tree_faces] = surface_index; tree_is_linked = 1; for (int i_face_edges = 0; i_face_edges < num_face_edges; ++i_face_edges) { /* We lock the edges of the face for surfaces, so that we do not link the same surface again - * to the edges of the face */ + * to the edges of the face */ if (dim == 2) /* 2D */ { edge_geometries[i_face_edges + t8_eclass_num_edges[eclass]] = -1; @@ -1143,21 +1088,21 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t parameters[i_face_nodes * 2] = face_nodes[i_face_nodes].parameters[0]; parameters[i_face_nodes * 2 + 1] = face_nodes[i_face_nodes].parameters[1]; } - /* Corrects the parameters on the surface if it is closed to prevent disorted elements. */ - for (int param_dim = 0; param_dim < 2; ++param_dim) { - if (cad_geometry->get_cad_manager ()->t8_geom_is_surface_closed (surface_index, param_dim)) { - t8_cmesh_correct_parameters_on_closed_geometry (2, surface_index, num_face_nodes, cad_geometry, parameters); - } - } t8_cmesh_set_attribute (cmesh, tree_count, t8_get_package_id (), T8_CMESH_CAD_FACE_PARAMETERS_ATTRIBUTE_KEY + i_tree_faces, parameters, num_face_nodes * 2 * sizeof (double), 0); } } + /*----------------------------------------- End of face-surface linkage -----------------------------------------*/ + + /*----------------------------------------- Start of edge linkage -----------------------------------------*/ const int num_edges = t8_eclass_num_edges[eclass]; - /* Then we look for geometries linked to the edges */ + for (int i_tree_edges = 0; i_tree_edges < num_edges; ++i_tree_edges) { + std::array edge_nodes; + /* Save edge nodes separately. We have to distinct between eclass dimension because + * the geometrical edge of a 2D element is counted as face for flux computation. */ if (t8_eclass_to_dimension[eclass] == 3) { edge_nodes[0] = tree_nodes[t8_edge_vertex_to_tree_vertex[eclass][i_tree_edges][0]]; edge_nodes[1] = tree_nodes[t8_edge_vertex_to_tree_vertex[eclass][i_tree_edges][1]]; @@ -1172,7 +1117,7 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t continue; } /* An edge can be linked to a curve as well as a surface. - * Therefore, we have to save the geometry dim and tag */ + * Therefore, we have to save the geometry dim and tag */ int edge_geometry_dim = 0; int edge_geometry_tag = 0; /* We check which is the highest dim a node geometry has and what is its tag */ @@ -1194,8 +1139,8 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t continue; } - /* If one vertex lies on a geometry of a higher dim as the other, we have to check, - * if the geometry of lower dimension is on that geometry. */ + /* If one node lies on a geometry of a higher dim as the other, we have to check, + * if the geometry of lower dimension is also on that same geometry. */ { int is_on_geom = 1; for (int i_edge = 0; i_edge < 2; ++i_edge) { @@ -1221,47 +1166,91 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t } } } + /* It is still possible that one node is on an edge and the other on a vertex and if so, + * they could share a surface. + common + surface + | | + | x node on edge + | /| + | edge / | + | / | + | / | + | / | + node on vertex x----------x vertex + */ if (!is_on_geom) { - continue; + int common_face = 0; + if (edge_nodes[0].entity_dim == 1 && edge_nodes[1].entity_dim == 0) + common_face = cad_geometry->get_cad_manager ()->t8_geom_get_common_face_of_vertex_and_edge ( + edge_nodes[1].entity_tag, edge_nodes[0].entity_tag); + if (edge_nodes[1].entity_dim == 1 && edge_nodes[0].entity_dim == 0) + common_face = cad_geometry->get_cad_manager ()->t8_geom_get_common_face_of_vertex_and_edge ( + edge_nodes[0].entity_tag, edge_nodes[1].entity_tag); + /* If a common face exists we can save it. */ + if (common_face > 0) { + edge_geometry_dim = 2; + edge_geometry_tag = common_face; + } + else + /* Otherwise we discard the whole edge */ + continue; } } - /* If both nodes are on a vertex we still got no edge. - * But we can look if both vertices share an edge and use this edge. - * If not we can skip this edge. */ + /* If both nodes are on a vertex we still got no edge or face. + * But we can look if both vertices share an edge/face and use this edge/face. + * If not we can skip this edge. */ if (edge_geometry_dim == 0 && edge_geometry_tag == 0) { - int common_curve = cad_geometry->get_cad_manager ()->t8_geom_get_common_edge (edge_nodes[0].entity_tag, - edge_nodes[1].entity_tag); - if (common_curve > 0) { - edge_geometry_tag = common_curve; + int common_edge = cad_geometry->get_cad_manager ()->t8_geom_get_common_edge_of_vertices ( + edge_nodes[0].entity_tag, edge_nodes[1].entity_tag); + if (common_edge > 0) { + edge_geometry_tag = common_edge; edge_geometry_dim = 1; } else { - continue; + int common_face = cad_geometry->get_cad_manager ()->t8_geom_get_common_face_of_vertices ( + edge_nodes[0].entity_tag, edge_nodes[1].entity_tag); + if (common_face > 0) { + edge_geometry_tag = common_face; + edge_geometry_dim = 2; + } + else { + continue; + } } } /* If both nodes are on different edges we have to look if both edges share a surface. - * If not we can skip this edge */ + * If not we can skip this edge */ if (edge_nodes[0].entity_dim == 1 && edge_nodes[1].entity_dim == 1 && edge_nodes[0].entity_tag != edge_nodes[1].entity_tag) { - int common_surface = cad_geometry->get_cad_manager ()->t8_geom_get_common_face (edge_nodes[0].entity_tag, - edge_nodes[1].entity_tag); - if (common_surface > 0) { - edge_geometry_tag = common_surface; + int common_face = cad_geometry->get_cad_manager ()->t8_geom_get_common_face_of_edges (edge_nodes[0].entity_tag, + edge_nodes[1].entity_tag); + if (common_face > 0) { + edge_geometry_tag = common_face; edge_geometry_dim = 2; } else { continue; } } + + /*----------------------------------------- Start of edge-curve linkage -----------------------------------------*/ /* If we have found a curve we can look for the parameters */ if (edge_geometry_dim == 1) { + + /* Abort if the edge is a line. If the edge is a line, the curved computation would + * only result in more computational overhead. */ + if (cad_geometry->get_cad_manager ()->t8_geom_is_line (edge_geometry_tag)) { + continue; + } + /* Check if adjacent faces carry a surface and if this edge lies on the surface */ for (int i_adjacent_face = 0; i_adjacent_face < 2; ++i_adjacent_face) { if (face_geometries[t8_edge_to_face[eclass][i_tree_edges][i_adjacent_face]] > 0) { if (!cad_geometry->get_cad_manager ()->t8_geom_is_edge_on_face ( edge_geometry_tag, face_geometries[t8_edge_to_face[eclass][i_tree_edges][i_adjacent_face]])) { - t8_global_errorf ("Error: Adjacent edge and face of a tree carry " + t8_global_errorf ("Error during mesh-cad recombination: Adjacent edge and face of a tree carry " "incompatible geometries.\n"); return 0; } @@ -1270,13 +1259,13 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t for (int i_edge_node = 0; i_edge_node < 2; ++i_edge_node) { /* Some error checking */ if (edge_nodes[i_edge_node].entity_dim == 2) { - t8_global_errorf ("Error: Node %li should lie on a vertex or an edge, " + t8_global_errorf ("Error during mesh-cad recombination: Node %li should lie on a vertex or an edge, " "but it lies on a surface.\n", edge_nodes[i_edge_node].index); return 0; } if (edge_nodes[i_edge_node].entity_dim == 1 && edge_nodes[i_edge_node].entity_tag != edge_geometry_tag) { - t8_global_errorf ("Error: Node %li should lie on a specific edge, " + t8_global_errorf ("Error during mesh-cad recombination: Node %li should lie on a specific edge, " "but it lies on another edge.\n", edge_nodes[i_edge_node].index); return 0; @@ -1284,9 +1273,10 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t if (edge_nodes[i_edge_node].entity_dim == 0) { if (!cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_edge (edge_nodes[i_edge_node].entity_tag, edge_geometry_tag)) { - t8_global_errorf ("Error: Node %li should lie on a vertex which lies on an edge, " - "but the vertex does not lie on that edge.\n", - edge_nodes[i_edge_node].index); + t8_global_errorf ( + "Error during mesh-cad recombination: Node %li should lie on a vertex which lies on an edge, " + "but the vertex does not lie on that edge.\n", + edge_nodes[i_edge_node].index); return 0; } } @@ -1294,56 +1284,60 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t /* If the node lies on a vertex we retrieve its parameter on the curve */ if (edge_nodes[i_edge_node].entity_dim == 0) { cad_geometry->get_cad_manager ()->t8_geom_get_parameter_of_vertex_on_edge ( - edge_nodes[i_edge_node].entity_tag, edge_geometry_tag, edge_nodes[i_edge_node].parameters.data ()); + edge_nodes[i_edge_node].entity_tag, edge_geometry_tag, edge_nodes[i_edge_node].parameters.data (), + edge_nodes[(i_edge_node + 1) % 2].parameters[0]); edge_nodes[i_edge_node].entity_dim = 1; } } - /* Abort if the edge is a line */ - if (cad_geometry->get_cad_manager ()->t8_geom_is_line (edge_geometry_tag)) { - continue; - } - edge_geometries[i_tree_edges] = edge_geometry_tag; tree_is_linked = 1; parameters[0] = edge_nodes[0].parameters[0]; parameters[1] = edge_nodes[1].parameters[0]; - /* Corrects the parameters on the edge if it is closed to prevent disorted elements. */ - if (cad_geometry->get_cad_manager ()->t8_geom_is_edge_closed (edge_geometry_tag)) { - t8_cmesh_correct_parameters_on_closed_geometry (1, edge_geometry_tag, 2, cad_geometry, parameters); - } - t8_cmesh_set_attribute (cmesh, tree_count, t8_get_package_id (), T8_CMESH_CAD_EDGE_PARAMETERS_ATTRIBUTE_KEY + i_tree_edges, parameters, 2 * sizeof (double), 0); } + /*----------------------------------------- End of edge-curve linkage -----------------------------------------*/ + + /*----------------------------------------- Start of edge-surface linkage -----------------------------------------*/ /* If we have found a surface we can look for the parameters. - * If the edge is locked for edges on surfaces we have to skip this edge */ + * If the edge is locked for edges on surfaces we have to skip this edge */ else if (edge_geometry_dim == 2 && edge_geometries[i_tree_edges + num_edges] >= 0) { + + /* We also skip this exge if the edge is on a plane. Planes have no curvature and + * therefore only result in computational overhead. */ + if (cad_geometry->get_cad_manager ()->t8_geom_is_plane (edge_geometry_tag)) { + continue; + } + /* If the node lies on a geometry with a different dimension we try to retrieve the parameters */ for (int i_edge_node = 0; i_edge_node < 2; ++i_edge_node) { /* Some error checking */ if (edge_nodes[i_edge_node].entity_dim == 2 && edge_nodes[i_edge_node].entity_tag != edge_geometry_tag) { - t8_global_errorf ("Error: Node %li should lie on a specific face, but it lies on another face.\n", + t8_global_errorf ("Error during mesh-cad recombination: Node %li should lie on a specific face, but it lies " + "on another face.\n", edge_nodes[i_edge_node].index); return 0; } if (edge_nodes[i_edge_node].entity_dim == 0) { if (!cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_face (edge_nodes[i_edge_node].entity_tag, edge_geometry_tag)) { - t8_global_errorf ("Error: Node %li should lie on a vertex which lies on a face, " - "but the vertex does not lie on that face.\n", - edge_nodes[i_edge_node].index); + t8_global_errorf ( + "Error during mesh-cad recombination: Node %li should lie on a vertex which lies on a face, " + "but the vertex does not lie on that face.\n", + edge_nodes[i_edge_node].index); return 0; } } if (edge_nodes[i_edge_node].entity_dim == 1) { if (!cad_geometry->get_cad_manager ()->t8_geom_is_edge_on_face (edge_nodes[i_edge_node].entity_tag, edge_geometry_tag)) { - t8_global_errorf ("Error: Node %li should lie on an edge which lies on a face, " - "but the edge does not lie on that face.\n", - edge_nodes[i_edge_node].index); + t8_global_errorf ( + "Error during mesh-cad recombination: Node %li should lie on an edge which lies on a face, " + "but the edge does not lie on that face.\n", + edge_nodes[i_edge_node].index); return 0; } } @@ -1356,19 +1350,13 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t } /* If the node lies on an edge we have to do the same */ if (edge_nodes[i_edge_node].entity_dim == 1) { - const int num_face_nodes = t8_eclass_num_vertices[eclass]; cad_geometry->get_cad_manager ()->t8_geom_edge_parameter_to_face_parameters ( - edge_nodes[i_edge_node].entity_tag, edge_geometry_tag, num_face_nodes, - edge_nodes[i_edge_node].parameters[0], parameters, edge_nodes[i_edge_node].parameters.data ()); + edge_nodes[i_edge_node].entity_tag, edge_geometry_tag, edge_nodes[i_edge_node].parameters[0], + edge_nodes[i_edge_node].parameters.data ()); edge_nodes[i_edge_node].entity_dim = 2; } } - /* Abort if the edge is a line */ - if (cad_geometry->get_cad_manager ()->t8_geom_is_line (edge_geometry_tag)) { - continue; - } - edge_geometries[i_tree_edges + t8_eclass_num_edges[eclass]] = edge_geometry_tag; tree_is_linked = 1; parameters[0] = edge_nodes[0].parameters[0]; @@ -1376,13 +1364,6 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t parameters[2] = edge_nodes[1].parameters[0]; parameters[3] = edge_nodes[1].parameters[1]; - /* Corrects the parameters on the surface if it is closed to prevent disorted elements. */ - for (int param_dim = 0; param_dim < 2; ++param_dim) { - if (cad_geometry->get_cad_manager ()->t8_geom_is_surface_closed (edge_geometry_tag, param_dim)) { - t8_cmesh_correct_parameters_on_closed_geometry (2, edge_geometry_tag, 2, cad_geometry, parameters); - } - } - t8_cmesh_set_attribute (cmesh, tree_count, t8_get_package_id (), T8_CMESH_CAD_EDGE_PARAMETERS_ATTRIBUTE_KEY + i_tree_edges, parameters, 4 * sizeof (double), 0); @@ -1409,6 +1390,7 @@ t8_cmesh_process_tree_geometry (t8_cmesh_t cmesh, t8_eclass_t eclass, int dim, t t8_debugf ("Registering tree %li with geometry %s \n", tree_count, linear_geometry_base->t8_geom_get_name ().c_str ()); } + return 1; } #endif /* T8_ENABLE_OCC */ @@ -1432,10 +1414,6 @@ t8_cmesh_msh_file_4_read_eles (t8_cmesh_t cmesh, FILE *fp, const t8_msh_node_tab t8_gloidx_t tree_count; t8_eclass_t eclass; t8_msh_file_node Node; -#if T8_ENABLE_OCC - std::array face_nodes; - std::array edge_nodes; -#endif /* T8_ENABLE_OCC */ long lnum_trees, lnum_blocks, entity_tag; int retval; int ele_type; @@ -1603,7 +1581,7 @@ t8_cmesh_msh_file_4_read_eles (t8_cmesh_t cmesh, FILE *fp, const t8_msh_node_tab else { #if T8_ENABLE_OCC if (!t8_cmesh_process_tree_geometry (cmesh, eclass, dim, tree_count, cad_geometry_base, linear_geometry_base, - tree_nodes, face_nodes, edge_nodes)) { + tree_nodes)) { free (line); t8_cmesh_destroy (&cmesh); return std::nullopt; diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.cxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.cxx index 0dabe52934..3866b2720d 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.cxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.cxx @@ -213,7 +213,6 @@ t8_geometry_cad::t8_geom_evaluate_cad_tri (t8_cmesh_t cmesh, t8_gloidx_t gtreeid cmesh, t8_get_package_id (), T8_CMESH_CAD_EDGE_PARAMETERS_ATTRIBUTE_KEY + i_edge, ltreeid); T8_ASSERT (edge_parameters != NULL); - const int num_face_nodes = t8_eclass_num_vertices[active_tree_class]; /* Calculate the curve parameter at the intersection point for each reference coordinate */ for (size_t i_coord = 0; i_coord < num_coords; ++i_coord) { const int offset_2d = i_coord * 2; @@ -222,9 +221,9 @@ t8_geometry_cad::t8_geom_evaluate_cad_tri (t8_cmesh_t cmesh, t8_gloidx_t gtreeid t8_geom_linear_interpolation (&ref_intersection[(i_edge == 0) + offset_2d], edge_parameters, 1, 1, &interpolated_curve_parameter); /* Convert the interpolated edge parameter of each reference point to surface parameters */ - cad_manager->t8_geom_edge_parameter_to_face_parameters (edges[i_edge], *faces, num_face_nodes, - interpolated_curve_parameter, face_parameters, - converted_edge_surface_parameters + offset_2d); + cad_manager->t8_geom_edge_parameter_to_face_parameters (edges[i_edge], *faces, interpolated_curve_parameter, + converted_edge_surface_parameters + offset_2d, + std::span (face_parameters, 2)); } double edge_surface_parameters[4]; @@ -377,7 +376,6 @@ t8_geometry_cad::t8_geom_evaluate_cad_quad (t8_cmesh_t cmesh, t8_gloidx_t gtreei */ const int edge_orthogonal_direction = (i_edge >> 1); const int edge_direction = 1 - edge_orthogonal_direction; - const int num_face_nodes = t8_eclass_num_vertices[active_tree_class]; double temp_edge_parameters[2]; double temp_face_parameters[2]; /* Retrieve edge parameters and interpolate */ @@ -393,11 +391,10 @@ t8_geometry_cad::t8_geom_evaluate_cad_quad (t8_cmesh_t cmesh, t8_gloidx_t gtreei t8_geom_linear_interpolation (&ref_coords[edge_direction + offset_2d], edge_parameters, 1, 1, temp_edge_parameters); - pnt = process_curve (edges[i_edge], temp_edge_parameters[0]); - /* Convert curve parameter to surface parameters */ - cad_manager->t8_geom_edge_parameter_to_face_parameters ( - edges[i_edge], *faces, num_face_nodes, temp_edge_parameters[0], face_parameters, temp_edge_parameters); + cad_manager->t8_geom_edge_parameter_to_face_parameters (edges[i_edge], *faces, temp_edge_parameters[0], + temp_edge_parameters, + std::span (face_parameters, 2)); /* Interpolate between the surface parameters of the current edge */ double edge_surface_parameters[4]; @@ -618,7 +615,6 @@ t8_geometry_cad::t8_geom_evaluate_cad_tet (t8_cmesh_t cmesh, t8_gloidx_t gtreeid /* Linear interpolation between parameters */ t8_geom_linear_interpolation (&interpolation_coeff, parameters, 2, 1, interpolated_surface_params); - T8_ASSERT (edges[i_edge + num_edges] <= cad_manager->t8_geom_get_cad_shape_face_map ().Size ()); pnt = process_surface (edges[i_edge + num_edges], interpolated_surface_params); } @@ -735,7 +731,7 @@ t8_geometry_cad::t8_geom_evaluate_cad_tet (t8_cmesh_t cmesh, t8_gloidx_t gtreeid interpolated_coords[dim] += face_displacement_from_edges[dim]; } - pnt = process_surface (faces[i_faces], interpolated_surface_params); + pnt = process_surface (faces[i_faces], interpolated_surface_parameters); /* Compute the scaling factor. The scaling happens along the straight from * the opposite vertex of the face to the face_intersection. */ @@ -965,10 +961,10 @@ t8_geometry_cad::t8_geom_evaluate_cad_hex (t8_cmesh_t cmesh, t8_gloidx_t gtreeid } } /* Convert the interpolated parameter of the curve into the corresponding parameters on the surface */ - const int num_face_nodes = t8_eclass_num_vertices[T8_ECLASS_QUAD]; cad_manager->t8_geom_edge_parameter_to_face_parameters ( - edges[t8_face_edge_to_tree_edge[T8_ECLASS_HEX][i_faces][i_face_edge]], faces[i_faces], num_face_nodes, - interpolated_curve_param, surface_parameters, surface_parameters_from_curve); + edges[t8_face_edge_to_tree_edge[T8_ECLASS_HEX][i_faces][i_face_edge]], faces[i_faces], + interpolated_curve_param, surface_parameters_from_curve, + std::span (surface_parameters, 2)); /* Calculate the displacement between the interpolated parameters on the surface * and the parameters on the surface converted from the parameter of the curve From 83e95dc4fc4c6d3cfadd341d12084f3a0f0eb3bc Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Fri, 9 Jan 2026 10:11:14 +0100 Subject: [PATCH 03/22] add debug check --- src/t8_cad/t8_cad.cxx | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index 05e4a56ed3..81cf1deef3 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -35,6 +35,10 @@ #include #include +#if T8_ENABLE_DEBUG +#include +#endif + t8_cad::t8_cad (std::string fileprefix) { BRep_Builder builder; @@ -272,6 +276,14 @@ t8_cad::t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const i bool first_point = true; for (TopExp_Explorer dora (edge, TopAbs_VERTEX); dora.More (); dora.Next ()) { const TopoDS_Vertex current_vertex = TopoDS::Vertex (dora.Current ()); + +#if T8_ENABLE_DEBUG + /* Check of point is really the same. */ + const gp_Pnt debug_reference_point = BRep_Tool::Pnt (vertex); + const gp_Pnt debug_current_point = BRep_Tool::Pnt (current_vertex); + T8_ASSERT (debug_reference_point.Distance (debug_current_point) <= Precision::Confusion ()); +#endif + if (first_point) { *edge_param = BRep_Tool::Parameter (current_vertex, edge); first_point = false; From 94456c4ace8763278cef545858b0d962e2ba518c Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Fri, 9 Jan 2026 10:11:39 +0100 Subject: [PATCH 04/22] documentation --- src/t8_cad/t8_cad.cxx | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index 81cf1deef3..ece90871d3 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -271,8 +271,17 @@ t8_cad::t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const i TopoDS_Vertex vertex = TopoDS::Vertex (cad_shape_vertex_map.FindKey (vertex_index)); TopoDS_Edge edge = TopoDS::Edge (cad_shape_edge_map.FindKey (edge_index)); + /* If the edge is not closed or the user did not provide any reference coords we just query + the parameter of the vertex on the edge. But if the edge is closed and the user provided + said reference parameters it gets more sophisticated:*/ + const bool edge_is_closed = t8_cad::t8_geom_edge_is_closed (edge_index); if (reference_edge_param && edge_is_closed) { + /* Edge is closed and the user provided reference parameters. We iterate over all vertices of the edge. + Since the edge is closed, the start and end vertex should be in the same physical location. + We choose the point where the parameters are closer to the reference parameters. For debugging reasons + we also check, if the points are really in the same physical location. */ + bool first_point = true; for (TopExp_Explorer dora (edge, TopAbs_VERTEX); dora.More (); dora.Next ()) { const TopoDS_Vertex current_vertex = TopoDS::Vertex (dora.Current ()); From 634ef3778b506deb6d1a5cca2e141513aebd425b Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Tue, 13 Jan 2026 16:17:47 +0100 Subject: [PATCH 05/22] fix documentation --- src/t8_cad/t8_cad.hxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/t8_cad/t8_cad.hxx b/src/t8_cad/t8_cad.hxx index e17fa18fa9..ab2f8d24e4 100644 --- a/src/t8_cad/t8_cad.hxx +++ b/src/t8_cad/t8_cad.hxx @@ -227,10 +227,10 @@ class t8_cad { * it is random which side of the closed edge the parameter is from. The parameter * is correct, but for curved elements it has to be checked if the parameter has to be * converted onto the other bound. - * \param [in] vertex_index The index of the cad vertex. - * \param [in] edge_index The index of the cad edge. - * \param [out] edge_param The parameter of the vertex on the edge. - * \param [in] reference_face_params Reference parameters on the surface. + * \param [in] vertex_index The index of the cad vertex. + * \param [in] edge_index The index of the cad edge. + * \param [out] edge_param The parameter of the vertex on the edge. + * \param [in] reference_edge_param Reference parameters on the edge. */ void t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const int edge_index, double *edge_param, From 1bac59c68c0d78195efee11cd72e7fef753613a6 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Tue, 13 Jan 2026 16:18:15 +0100 Subject: [PATCH 06/22] fix documentation --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 8e0a0b0c32..600b291b4c 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -110,7 +110,7 @@ const int t8_msh_tree_vertex_to_t8_vertex_num[T8_ECLASS_COUNT][8] = { * the new number of bytes is stored in n. * \param [in] fp The file stream to read from. * \return The number of read arguments of the last line read. - * negative on failure + * negative on failure */ static int t8_cmesh_msh_read_next_line (char **line, size_t *n, FILE *fp) @@ -844,11 +844,7 @@ t8_store_element_node_data (t8_cmesh_t cmesh, t8_gloidx_t tree_count, * \param [in] cad_geometry_base A pointer to the CAD-based geometry object. * \param [in] linear_geometry_base A pointer to the linear geometry object. * \param [in] tree_nodes An array of nodes representing the vertices of the tree. - * \param [in] face_nodes An array of nodes representing the faces of the tree. - * \param [in] edge_nodes An array of nodes representing the edges of the tree. - * * \return True if the tree geometry was successfully processed; false otherwise. - * */ static bool t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass, const int dim, From cd00821d4e3d34b4735f55329bd4e2789dccbcb2 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Thu, 15 Jan 2026 14:43:32 +0100 Subject: [PATCH 07/22] remove unused function --- src/t8.h | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/src/t8.h b/src/t8.h index 5477970880..4ac4b947af 100644 --- a/src/t8.h +++ b/src/t8.h @@ -305,24 +305,4 @@ t8_sc_array_index_locidx (const sc_array_t *array, const t8_locidx_t index); /** Call this at the end of a header file to match T8_EXTERN_C_BEGIN (). */ T8_EXTERN_C_END (); -#ifdef __cplusplus - -/** - * Check two floating point values for equality. - * \tparam TFloating1 Type of the first floating point value. - * \tparam TFloating2 Type of the second floating point value. - * \param [in] floating1 First floating point value to check for equality. - * \param [in] floating2 Second floating point value to check for equality. - * \param [in] tolerance The tolerance. - * \return True if floating values are equal, false otherwise. - */ -template -inline bool -t8_floating_is_equal (TFloating1 floating1, TFloating2 floating2, double tolerance = T8_PRECISION_EPS) noexcept -{ - return std::abs (floating1 - floating2) <= tolerance; -} - -#endif /* !__cplusplus */ - #endif /* !T8_H */ From 02d743c06bf030a54f5d9da92d518b136ce04606 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Thu, 15 Jan 2026 16:18:00 +0100 Subject: [PATCH 08/22] add assertions --- src/t8_cad/t8_cad.cxx | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index ece90871d3..e71489a86c 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -156,6 +156,8 @@ t8_cad::t8_geom_get_cad_shape_face_map () const int t8_cad::t8_geom_get_common_edge_of_vertices (const int vertex1_index, const int vertex2_index) const { + T8_ASSERT (vertex1_index <= cad_shape_vertex2edge_map.Size ()); + T8_ASSERT (vertex2_index <= cad_shape_vertex2edge_map.Size ()); const TopTools_ListOfShape collection1 = cad_shape_vertex2edge_map.FindFromIndex (vertex1_index); const TopTools_ListOfShape collection2 = cad_shape_vertex2edge_map.FindFromIndex (vertex2_index); @@ -172,6 +174,8 @@ t8_cad::t8_geom_get_common_edge_of_vertices (const int vertex1_index, const int int t8_cad::t8_geom_get_common_face_of_edges (const int edge1_index, const int edge2_index) const { + T8_ASSERT (edge1_index <= cad_shape_edge2face_map.Size ()); + T8_ASSERT (edge2_index <= cad_shape_edge2face_map.Size ()); const TopTools_ListOfShape collection1 = cad_shape_edge2face_map.FindFromIndex (edge1_index); const TopTools_ListOfShape collection2 = cad_shape_edge2face_map.FindFromIndex (edge2_index); @@ -188,6 +192,8 @@ t8_cad::t8_geom_get_common_face_of_edges (const int edge1_index, const int edge2 int t8_cad::t8_geom_get_common_face_of_vertex_and_edge (const int vertex_index, const int edge_index) const { + T8_ASSERT (vertex_index <= cad_shape_vertex_map.Size ()); + T8_ASSERT (edge_index <= cad_shape_edge2face_map.Size ()); const TopTools_ListOfShape edge_collection = cad_shape_edge2face_map.FindFromIndex (edge_index); for (auto face = edge_collection.begin (); face != edge_collection.end (); ++face) { const size_t face_index = cad_shape_face_map.FindIndex (*face); @@ -201,6 +207,8 @@ t8_cad::t8_geom_get_common_face_of_vertex_and_edge (const int vertex_index, cons int t8_cad::t8_geom_get_common_face_of_vertices (const int vertex1_index, const int vertex2_index) const { + T8_ASSERT (vertex1_index <= cad_shape_vertex2face_map.Size ()); + T8_ASSERT (vertex2_index <= cad_shape_vertex2face_map.Size ()); const TopTools_ListOfShape collection1 = cad_shape_vertex2face_map.FindFromIndex (vertex1_index); const TopTools_ListOfShape collection2 = cad_shape_vertex2face_map.FindFromIndex (vertex2_index); @@ -217,6 +225,8 @@ t8_cad::t8_geom_get_common_face_of_vertices (const int vertex1_index, const int int t8_cad::t8_geom_is_vertex_on_edge (const int vertex_index, const int edge_index) const { + T8_ASSERT (vertex_index <= cad_shape_vertex2edge_map.Size ()); + T8_ASSERT (edge_index <= cad_shape_edge_map.Size ()); const TopTools_ListOfShape collection = cad_shape_vertex2edge_map.FindFromIndex (vertex_index); return collection.Contains (cad_shape_edge_map.FindKey (edge_index)); } @@ -224,6 +234,8 @@ t8_cad::t8_geom_is_vertex_on_edge (const int vertex_index, const int edge_index) int t8_cad::t8_geom_is_edge_on_face (const int edge_index, const int face_index) const { + T8_ASSERT (edge_index <= cad_shape_edge2face_map.Size ()); + T8_ASSERT (face_index <= cad_shape_face_map.Size ()); const TopTools_ListOfShape collection = cad_shape_edge2face_map.FindFromIndex (edge_index); return collection.Contains (cad_shape_face_map.FindKey (face_index)); } @@ -231,6 +243,8 @@ t8_cad::t8_geom_is_edge_on_face (const int edge_index, const int face_index) con int t8_cad::t8_geom_is_vertex_on_face (const int vertex_index, const int face_index) const { + T8_ASSERT (vertex_index <= cad_shape_vertex2edge_map.Size ()); + T8_ASSERT (face_index <= cad_shape_face_map.Size ()); const TopTools_ListOfShape edge_collection = cad_shape_vertex2edge_map.FindFromIndex (vertex_index); for (auto edge = edge_collection.begin (); edge != edge_collection.end (); ++edge) { const TopTools_ListOfShape face_collection = cad_shape_edge2face_map.FindFromKey (*edge); @@ -244,6 +258,8 @@ t8_cad::t8_geom_is_vertex_on_face (const int vertex_index, const int face_index) bool t8_cad::t8_geom_vertex_is_on_seam (const int vertex_index, const int face_index) const { + T8_ASSERT (vertex_index <= cad_shape_vertex2edge_map.Size ()); + T8_ASSERT (face_index <= cad_shape_face_map.Size ()); const auto face = t8_cad::t8_geom_get_cad_face (face_index); const TopTools_ListOfShape edge_collection = cad_shape_vertex2edge_map.FindFromIndex (vertex_index); const ShapeAnalysis_Edge edge_analyzer; @@ -257,6 +273,8 @@ t8_cad::t8_geom_vertex_is_on_seam (const int vertex_index, const int face_index) bool t8_cad::t8_geom_edge_is_seam (const int edge_index, const int face_index) const { + T8_ASSERT (edge_index <= cad_shape_edge_map.Size ()); + T8_ASSERT (face_index <= cad_shape_face_map.Size ()); const auto face = t8_cad::t8_geom_get_cad_face (face_index); const auto edge = t8_cad::t8_geom_get_cad_edge (edge_index); const ShapeAnalysis_Edge edge_analyzer; From acd1961dcf1b0f5e3a7458086c49e45c4164579a Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Thu, 15 Jan 2026 16:21:28 +0100 Subject: [PATCH 09/22] Apply suggestions from code review Typos Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> --- src/t8_cad/t8_cad.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index e71489a86c..f89ff5d423 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -332,10 +332,10 @@ void t8_cad::t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const int face_index, double face_params[2], std::optional> reference_face_params) const { - /* DISCLAIMER: This function is overly complicated and I do not understand why the simpler versions to not work. + /* DISCLAIMER: This function is overly complicated and I do not understand why the simpler versions do not work. The overly complicated part is only the edge case where the vertex lies on a seam and reference parameters are provided. The (simpler) plan was as follows: We have the vertex and face. But in the topology of the face, the vertex exists two times, - on the one side of the face and on the other. Bot sides are on the same physical location, because it is a seam. Thats + on the one side of the face and on the other. Both sides are on the same physical location, because it is a seam. Thats why both vertices have the same coordinates. But they do not have the same parameters on the surface. The simple plan is to use a TopExp_Explorer to iterate over all vertices of the face and to look at all vertices with the right coords. From all those vertices we choose the one, which has the closest parameters to our reference parameters. From 39d087ae3610807c40f2bc25102482c46f8658a4 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Thu, 15 Jan 2026 16:20:28 +0100 Subject: [PATCH 10/22] use optional.has_value() for better readability --- src/t8_cad/t8_cad.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index f89ff5d423..e0d291f269 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -294,7 +294,7 @@ t8_cad::t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const i said reference parameters it gets more sophisticated:*/ const bool edge_is_closed = t8_cad::t8_geom_edge_is_closed (edge_index); - if (reference_edge_param && edge_is_closed) { + if (reference_edge_param.has_value () && edge_is_closed) { /* Edge is closed and the user provided reference parameters. We iterate over all vertices of the edge. Since the edge is closed, the start and end vertex should be in the same physical location. We choose the point where the parameters are closer to the reference parameters. For debugging reasons @@ -349,7 +349,7 @@ t8_cad::t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const TopoDS_Face face = TopoDS::Face (cad_shape_face_map.FindKey (face_index)); /* If the surface is not closed or the user did not provide any reference params we just query the parameters. */ - if (!t8_cad::t8_geom_surface_is_closed (face_index) || !reference_face_params) { + if (!t8_cad::t8_geom_surface_is_closed (face_index) || !reference_face_params.has_value ()) { uv.emplace (BRep_Tool::Parameters (vertex, face)); } @@ -437,7 +437,7 @@ t8_cad::t8_geom_edge_parameter_to_face_parameters ( TopoDS_Edge edge = TopoDS::Edge (cad_shape_edge_map.FindKey (edge_index)); const bool is_seam = t8_cad::t8_geom_edge_is_seam (edge_index, face_index); - if (is_seam && reference_face_params) { + if (is_seam && reference_face_params.has_value ()) { /* Convert reference parameters to OCCT point */ gp_Pnt2d reference_point (reference_face_params.value ()[0], reference_face_params.value ()[1]); From c5401c3ad38f891032e7490ef2e5b29ef89b8df9 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Tue, 3 Feb 2026 13:19:42 +0100 Subject: [PATCH 11/22] fix another cad seam bug --- src/t8_cad/t8_cad.cxx | 39 +++-- src/t8_cad/t8_cad.hxx | 20 ++- .../t8_cmesh_io/t8_cmesh_readmshfile.cxx | 135 +++++++++++++----- 3 files changed, 148 insertions(+), 46 deletions(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index e0d291f269..06fc0ab20d 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -256,7 +256,30 @@ t8_cad::t8_geom_is_vertex_on_face (const int vertex_index, const int face_index) } bool -t8_cad::t8_geom_vertex_is_on_seam (const int vertex_index, const int face_index) const +t8_cad::t8_geom_vertex_is_seam (const int vertex_index, const int edge_index) const +{ + T8_ASSERT (vertex_index <= cad_shape_vertex_map.Size ()); + T8_ASSERT (edge_index <= cad_shape_edge_map.Size ()); + const auto vertex = t8_cad::t8_geom_get_cad_vertex (vertex_index); + const auto edge = t8_cad::t8_geom_get_cad_edge (edge_index); + double u = BRep_Tool::Parameter (vertex, edge); + double first = 0, last = 0; + const auto curve = BRep_Tool::Curve (edge, first, last); + + /* If curve is not periodic, there is no seam. */ + if (!curve->IsPeriodic ()) + return false; + + // Seam is at the start/end of the trimmed range + double tol = Precision::PConfusion (first > last ? first : last); + const bool on_first = (std::abs (u - first) <= tol); + const bool on_last = (std::abs (u - last) <= tol); + + return on_first || on_last; +} + +bool +t8_cad::t8_geom_vertex_is_on_seam_edge (const int vertex_index, const int face_index) const { T8_ASSERT (vertex_index <= cad_shape_vertex2edge_map.Size ()); T8_ASSERT (face_index <= cad_shape_face_map.Size ()); @@ -301,15 +324,14 @@ t8_cad::t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const i we also check, if the points are really in the same physical location. */ bool first_point = true; + const gp_Pnt vertex_pnt = BRep_Tool::Pnt (vertex); for (TopExp_Explorer dora (edge, TopAbs_VERTEX); dora.More (); dora.Next ()) { const TopoDS_Vertex current_vertex = TopoDS::Vertex (dora.Current ()); -#if T8_ENABLE_DEBUG - /* Check of point is really the same. */ - const gp_Pnt debug_reference_point = BRep_Tool::Pnt (vertex); - const gp_Pnt debug_current_point = BRep_Tool::Pnt (current_vertex); - T8_ASSERT (debug_reference_point.Distance (debug_current_point) <= Precision::Confusion ()); -#endif + /* The underlying curve can be closed, but since the active part of the curve (the edge) is not checked + there can be points which do not match our provided point. */ + if (!vertex_pnt.IsEqual (BRep_Tool::Pnt (current_vertex), Precision::Confusion ())) + continue; if (first_point) { *edge_param = BRep_Tool::Parameter (current_vertex, edge); @@ -322,6 +344,7 @@ t8_cad::t8_geom_get_parameter_of_vertex_on_edge (const int vertex_index, const i *edge_param = other_param; } } + T8_ASSERT (first_point == false); } else { *edge_param = BRep_Tool::Parameter (vertex, edge); @@ -355,7 +378,7 @@ t8_cad::t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const /* If the vertex is not on the seam we can also just query the parameters. */ else { - const bool is_on_seam = t8_cad::t8_geom_vertex_is_on_seam (vertex_index, face_index); + const bool is_on_seam = t8_cad::t8_geom_vertex_is_on_seam_edge (vertex_index, face_index); if (!is_on_seam) { uv.emplace (BRep_Tool::Parameters (vertex, face)); } diff --git a/src/t8_cad/t8_cad.hxx b/src/t8_cad/t8_cad.hxx index ab2f8d24e4..92e33b4ccb 100644 --- a/src/t8_cad/t8_cad.hxx +++ b/src/t8_cad/t8_cad.hxx @@ -201,21 +201,31 @@ class t8_cad { int t8_geom_is_vertex_on_face (const int vertex_index, const int face_index) const; + /** Returns true if \a vertex_index is on a seam of \a edge_index. + * A seam is a vertex which connects a curve to itself. + * + * \param [in] vertex_index The index of the cad vertex. + * \param [in] edge_index The index of the cad edge. + * \return true if the vertex is a seam. false otherwise. + */ + bool + t8_geom_vertex_is_seam (const int vertex_index, const int edge_index) const; + /** Returns true if \a vertex_index is on a seam of \a face_index. * A seam is an edge which connects a surface to itself. * - * \param [in] vertex_index - * \param [in] face_index + * \param [in] vertex_index The index of the cad vertex. + * \param [in] face_index The index of the cad face. * \return true if the edge is a seam. false otherwise. */ bool - t8_geom_vertex_is_on_seam (const int vertex_index, const int face_index) const; + t8_geom_vertex_is_on_seam_edge (const int vertex_index, const int face_index) const; /** Returns true if \a edge_index is a seam of \a face_index. * A seam is an edge which connects a surface to itself. * - * \param [in] edge_index - * \param [in] face_index + * \param [in] edge_index The index of the cad edge. + * \param [in] face_index The index of the cad face. * \return true if the edge is a seam. false otherwise. */ bool diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 600b291b4c..081a36ac07 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -1099,7 +1099,6 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass /*----------------------------------------- Start of edge linkage -----------------------------------------*/ const int num_edges = t8_eclass_num_edges[eclass]; - for (int i_tree_edges = 0; i_tree_edges < num_edges; ++i_tree_edges) { std::array edge_nodes; /* Save edge nodes separately. We have to distinct between eclass dimension because @@ -1257,37 +1256,69 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass } } } - for (int i_edge_node = 0; i_edge_node < 2; ++i_edge_node) { + + /* We now know that bot vertices are on the same edge. + Next, we iterate twice over both vertices. In the first iteration we do + some sanity checks and try to find a parameter which is not on a seam. + Because when one of the nodes is on a vertex which is on the seam of the edge, + there are two possible solutions for conversion. + In the second iteration, we convert the parameters using the reference param. */ + std::optional reference_param = std::nullopt; + for (const auto &edge_node : edge_nodes) { /* Some error checking */ - if (edge_nodes[i_edge_node].entity_dim == 2) { + if (edge_node.entity_dim == 2) { t8_global_errorf ("Error during mesh-cad recombination: Node %li should lie on a vertex or an edge, " "but it lies on a surface.\n", - edge_nodes[i_edge_node].index); + edge_node.index); return 0; } - if (edge_nodes[i_edge_node].entity_dim == 1 && edge_nodes[i_edge_node].entity_tag != edge_geometry_tag) { + if (edge_node.entity_dim == 1 && edge_node.entity_tag != edge_geometry_tag) { t8_global_errorf ("Error during mesh-cad recombination: Node %li should lie on a specific edge, " "but it lies on another edge.\n", - edge_nodes[i_edge_node].index); + edge_node.index); return 0; } - if (edge_nodes[i_edge_node].entity_dim == 0) { - if (!cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_edge (edge_nodes[i_edge_node].entity_tag, - edge_geometry_tag)) { + if (edge_node.entity_dim == 0) { + if (!cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_edge (edge_node.entity_tag, edge_geometry_tag)) { t8_global_errorf ( "Error during mesh-cad recombination: Node %li should lie on a vertex which lies on an edge, " "but the vertex does not lie on that edge.\n", - edge_nodes[i_edge_node].index); + edge_node.index); return 0; } } + if (!reference_param.has_value ()) { + switch (edge_node.entity_dim) { + case 1: + /* If the node is on the curve, we can directly use its parameter as reference. */ + reference_param = edge_node.parameters[0]; + break; + case 0: + /* If the node is on a vertex we can convert the parameters safely if the vertex os not the seam of the edge. */ + if (!cad_geometry->get_cad_manager ()->t8_geom_vertex_is_seam (edge_node.entity_tag, edge_geometry_tag)) { + reference_param.emplace (); + cad_geometry->get_cad_manager ()->t8_geom_get_parameter_of_vertex_on_edge ( + edge_node.entity_tag, edge_geometry_tag, &reference_param.value ()); + } + break; + default: + SC_ABORT_NOT_REACHED (); + } + } + } + + if (!reference_param.has_value ()) { + t8_global_errorf ("Error during mesh-cad recombination: Reference parameter on curve not found.\n"); + return 0; + } + + for (auto &edge_node : edge_nodes) { /* If the node lies on a vertex we retrieve its parameter on the curve */ - if (edge_nodes[i_edge_node].entity_dim == 0) { + if (edge_node.entity_dim == 0) { cad_geometry->get_cad_manager ()->t8_geom_get_parameter_of_vertex_on_edge ( - edge_nodes[i_edge_node].entity_tag, edge_geometry_tag, edge_nodes[i_edge_node].parameters.data (), - edge_nodes[(i_edge_node + 1) % 2].parameters[0]); - edge_nodes[i_edge_node].entity_dim = 1; + edge_node.entity_tag, edge_geometry_tag, edge_node.parameters.data (), reference_param); + edge_node.entity_dim = 1; } } @@ -1313,48 +1344,86 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass continue; } - /* If the node lies on a geometry with a different dimension we try to retrieve the parameters */ - for (int i_edge_node = 0; i_edge_node < 2; ++i_edge_node) { + /* We now know that bot vertices are on the same edge. + Next, we iterate twice over both vertices. In the first iteration we do + some sanity checks and try to find a parameter which is not on a seam. + Because when one of the nodes is on a vertex which is on the seam of the edge, + there are two possible solutions for conversion. + In the second iteration, we convert the parameters using the reference param. */ + std::optional> reference_params = std::nullopt; + for (const auto &edge_node : edge_nodes) { /* Some error checking */ - if (edge_nodes[i_edge_node].entity_dim == 2 && edge_nodes[i_edge_node].entity_tag != edge_geometry_tag) { + if (edge_node.entity_dim == 2 && edge_node.entity_tag != edge_geometry_tag) { t8_global_errorf ("Error during mesh-cad recombination: Node %li should lie on a specific face, but it lies " "on another face.\n", - edge_nodes[i_edge_node].index); + edge_node.index); return 0; } - if (edge_nodes[i_edge_node].entity_dim == 0) { - if (!cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_face (edge_nodes[i_edge_node].entity_tag, - edge_geometry_tag)) { + if (edge_node.entity_dim == 0) { + if (!cad_geometry->get_cad_manager ()->t8_geom_is_vertex_on_face (edge_node.entity_tag, edge_geometry_tag)) { t8_global_errorf ( "Error during mesh-cad recombination: Node %li should lie on a vertex which lies on a face, " "but the vertex does not lie on that face.\n", - edge_nodes[i_edge_node].index); + edge_node.index); return 0; } } - if (edge_nodes[i_edge_node].entity_dim == 1) { - if (!cad_geometry->get_cad_manager ()->t8_geom_is_edge_on_face (edge_nodes[i_edge_node].entity_tag, - edge_geometry_tag)) { + if (edge_node.entity_dim == 1) { + if (!cad_geometry->get_cad_manager ()->t8_geom_is_edge_on_face (edge_node.entity_tag, edge_geometry_tag)) { t8_global_errorf ( "Error during mesh-cad recombination: Node %li should lie on an edge which lies on a face, " "but the edge does not lie on that face.\n", - edge_nodes[i_edge_node].index); + edge_node.index); return 0; } } + if (!reference_params.has_value ()) { + switch (edge_node.entity_dim) { + case 2: + /* If the node is on the surface, we can directly use its parameter as reference. */ + reference_params = edge_node.parameters; + break; + case 1: + /* If the node is on a curve, we can use the parameters if the curve is not theseam of the surface. */ + if (!cad_geometry->get_cad_manager ()->t8_geom_edge_is_seam (edge_node.entity_tag, edge_geometry_tag)) { + reference_params.emplace (); + cad_geometry->get_cad_manager ()->t8_geom_edge_parameter_to_face_parameters ( + edge_node.entity_tag, edge_geometry_tag, edge_node.parameters[0], reference_params.value ().data ()); + } + break; + case 0: + /* If the node is on a vertex we can convert the parameters safely if the vertex os not the seam of the edge. */ + if (!cad_geometry->get_cad_manager ()->t8_geom_vertex_is_on_seam_edge (edge_node.entity_tag, + edge_geometry_tag)) { + reference_params.emplace (); + cad_geometry->get_cad_manager ()->t8_geom_get_parameters_of_vertex_on_face ( + edge_node.entity_tag, edge_geometry_tag, reference_params.value ().data ()); + } + break; + default: + SC_ABORT_NOT_REACHED (); + } + } + } + + if (!reference_params.has_value ()) { + t8_global_errorf ("Error during mesh-cad recombination: Reference parameter on curve not found.\n"); + return 0; + } + for (auto &edge_node : edge_nodes) { /* If the node lies on a vertex we retrieve its parameters on the surface */ - if (edge_nodes[i_edge_node].entity_dim == 0) { + if (edge_node.entity_dim == 0) { cad_geometry->get_cad_manager ()->t8_geom_get_parameters_of_vertex_on_face ( - edge_nodes[i_edge_node].entity_tag, edge_geometry_tag, edge_nodes[i_edge_node].parameters.data ()); - edge_nodes[i_edge_node].entity_dim = 2; + edge_node.entity_tag, edge_geometry_tag, edge_node.parameters.data (), reference_params); + edge_node.entity_dim = 2; } /* If the node lies on an edge we have to do the same */ - if (edge_nodes[i_edge_node].entity_dim == 1) { + if (edge_node.entity_dim == 1) { cad_geometry->get_cad_manager ()->t8_geom_edge_parameter_to_face_parameters ( - edge_nodes[i_edge_node].entity_tag, edge_geometry_tag, edge_nodes[i_edge_node].parameters[0], - edge_nodes[i_edge_node].parameters.data ()); - edge_nodes[i_edge_node].entity_dim = 2; + edge_node.entity_tag, edge_geometry_tag, edge_node.parameters[0], edge_node.parameters.data (), + reference_params); + edge_node.entity_dim = 2; } } From 7e2d2e9bdeece6f48843e0aeb9ea2880a6891958 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Tue, 3 Feb 2026 13:40:37 +0100 Subject: [PATCH 12/22] Apply suggestions from code review Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> --- src/t8_cad/t8_cad.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index 06fc0ab20d..7da873f820 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -426,7 +426,7 @@ t8_cad::t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const const Handle_Geom2d_Curve pcurve_on_surface = BRep_Tool::CurveOnSurface (current_edge, face, first, last); /* We can now insert the queried parameter on the curve into the pcurve (the curve on the surface). If this is the first point we found uv is empty and we fill it. Otherwise, we define another point to - stor our parameters and we save the parameters in iv which are closer to the reference parameters. */ + store our parameters and we save the parameters in iv which are closer to the reference parameters. */ if (!uv) { uv.emplace (); pcurve_on_surface->D0 (curve_param, *uv); From fff389d38d355c37da2322e1ecc01a2aa553ff87 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Tue, 3 Feb 2026 13:56:04 +0100 Subject: [PATCH 13/22] typo --- src/t8_cad/t8_cad.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index 7da873f820..9f99fde26f 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -426,7 +426,7 @@ t8_cad::t8_geom_get_parameters_of_vertex_on_face (const int vertex_index, const const Handle_Geom2d_Curve pcurve_on_surface = BRep_Tool::CurveOnSurface (current_edge, face, first, last); /* We can now insert the queried parameter on the curve into the pcurve (the curve on the surface). If this is the first point we found uv is empty and we fill it. Otherwise, we define another point to - store our parameters and we save the parameters in iv which are closer to the reference parameters. */ + store our parameters and we save the parameters in uv which are closer to the reference parameters. */ if (!uv) { uv.emplace (); pcurve_on_surface->D0 (curve_param, *uv); From 193bc85cce636de7afe0ee40861526056524902a Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Tue, 3 Feb 2026 13:56:30 +0100 Subject: [PATCH 14/22] use an optional instead of a bool --- src/t8_cad/t8_cad.cxx | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index 9f99fde26f..bd28999b54 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -455,7 +455,7 @@ t8_cad::t8_geom_edge_parameter_to_face_parameters ( { T8_ASSERT (t8_cad::t8_geom_is_edge_on_face (edge_index, face_index)); Standard_Real first, last; - gp_Pnt2d uv; + std::optional uv = std::nullopt; TopoDS_Face face = TopoDS::Face (cad_shape_face_map.FindKey (face_index)); TopoDS_Edge edge = TopoDS::Edge (cad_shape_edge_map.FindKey (edge_index)); @@ -464,26 +464,25 @@ t8_cad::t8_geom_edge_parameter_to_face_parameters ( /* Convert reference parameters to OCCT point */ gp_Pnt2d reference_point (reference_face_params.value ()[0], reference_face_params.value ()[1]); - /* If the edge is a seam we have to get both points and check which one is closer. + /* If the edge is a seam we have to get both points and check which one has closer parameters. We iterate over all edges of the face and check if the edges are the same as the input edge. Then we check if the converted parameters are closer to the already computed parameters. */ - bool first_point = true; for (TopExp_Explorer dora (face, TopAbs_EDGE); dora.More (); dora.Next ()) { const TopoDS_Edge current_edge = TopoDS::Edge (dora.Current ()); /* Check if edge is one of the seams. */ if (edge.IsSame (current_edge)) { Handle_Geom2d_Curve curve_on_surface = BRep_Tool::CurveOnSurface (edge, face, first, last); /* If it is the first seam we compute the parameters. */ - if (first_point) { - curve_on_surface->D0 (edge_param, uv); - first_point = false; + if (!uv.has_value ()) { + uv.emplace (); + curve_on_surface->D0 (edge_param, uv.value ()); } /* If it is the second one we check which point is close. */ else { gp_Pnt2d other_point; curve_on_surface->D0 (edge_param, other_point); - if (other_point.Distance (reference_point) < uv.Distance (reference_point)) + if (other_point.Distance (reference_point) < uv.value ().Distance (reference_point)) uv = other_point; } } @@ -492,11 +491,15 @@ t8_cad::t8_geom_edge_parameter_to_face_parameters ( else { /* If the edge is not a seam or if no reference parameters were provided we just use the normal edge. */ + uv.emplace (); Handle_Geom2d_Curve curve_on_surface = BRep_Tool::CurveOnSurface (edge, face, first, last); - curve_on_surface->D0 (edge_param, uv); + curve_on_surface->D0 (edge_param, uv.value ()); + } + if (!uv.has_value ()) { + SC_ABORTF ("Failed to convert edge parameters to Face parameters."); } - face_params_out[0] = uv.X (); - face_params_out[1] = uv.Y (); + face_params_out[0] = uv.value ().X (); + face_params_out[1] = uv.value ().Y (); } void From de6d303fb78a47de77f1bfbe8e4a50764dae28ca Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Tue, 3 Feb 2026 15:23:37 +0100 Subject: [PATCH 15/22] remove unnecessary includes --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 3d5898456e..a843b87a12 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -37,8 +37,6 @@ #include #include #include -#include -#include #include #ifdef _WIN32 From f3c885c8cf5eecb3450269a8f678a5b6ab7cca30 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Tue, 3 Feb 2026 15:26:21 +0100 Subject: [PATCH 16/22] typo --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index a843b87a12..bdb78678f3 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -1091,7 +1091,7 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass * If the edge is locked for edges on surfaces we have to skip this edge */ else if (edge_geometry_dim == 2 && edge_geometries[i_tree_edges + num_edges] >= 0) { - /* We also skip this exge if the edge is on a plane. Planes have no curvature and + /* We also skip this edge if the edge is on a plane. Planes have no curvature and * therefore only result in computational overhead. */ if (cad_geometry->get_cad_manager ()->t8_geom_is_plane (edge_geometry_tag)) { continue; From 3462d36f16264a6ff2d9ee4b0caa8aa1ee965710 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Wed, 4 Feb 2026 14:18:30 +0100 Subject: [PATCH 17/22] add braces --- src/t8_cad/t8_cad.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/t8_cad/t8_cad.cxx b/src/t8_cad/t8_cad.cxx index bd28999b54..dc319f08c9 100644 --- a/src/t8_cad/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -267,8 +267,9 @@ t8_cad::t8_geom_vertex_is_seam (const int vertex_index, const int edge_index) co const auto curve = BRep_Tool::Curve (edge, first, last); /* If curve is not periodic, there is no seam. */ - if (!curve->IsPeriodic ()) + if (!curve->IsPeriodic ()) { return false; + } // Seam is at the start/end of the trimmed range double tol = Precision::PConfusion (first > last ? first : last); From 17c44c2daa8a18269c524289731636e395e4383a Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer Date: Wed, 4 Feb 2026 14:19:51 +0100 Subject: [PATCH 18/22] use bool instead of int --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index bdb78678f3..1afc638aa5 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -741,7 +741,7 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass } /* Now we can check if every node lies on the surface and retrieve its parameters */ if (surface_index) { - int all_nodes_on_surface = 1; + bool all_nodes_on_surface = 1; /* If the surface is closed in any direction we may need some reference parameters on the surface. These reference parameters can be used to find the right parameters if we need to convert seam edge parameters to the face parameters. */ From f29e1af51c31e3cc391f37af2d60890f979fc93c Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Wed, 4 Feb 2026 14:21:14 +0100 Subject: [PATCH 19/22] Update src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 1afc638aa5..65d1ac8966 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -1097,7 +1097,7 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass continue; } - /* We now know that bot vertices are on the same edge. + /* We now know that both vertices are on the same edge. Next, we iterate twice over both vertices. In the first iteration we do some sanity checks and try to find a parameter which is not on a seam. Because when one of the nodes is on a vertex which is on the seam of the edge, From e309e862ffae044002c0d161e06e7e7a07fbb8d1 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Wed, 4 Feb 2026 14:21:25 +0100 Subject: [PATCH 20/22] Update src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 65d1ac8966..0b4cc831e6 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -1010,7 +1010,7 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass } } - /* We now know that bot vertices are on the same edge. + /* We now know that both vertices are on the same edge. Next, we iterate twice over both vertices. In the first iteration we do some sanity checks and try to find a parameter which is not on a seam. Because when one of the nodes is on a vertex which is on the seam of the edge, From 61baaf990a7637ed0f8556acac332c30145d48fa Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Wed, 4 Feb 2026 14:21:35 +0100 Subject: [PATCH 21/22] Update src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 0b4cc831e6..6ce0a320ee 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -1145,7 +1145,7 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass } break; case 0: - /* If the node is on a vertex we can convert the parameters safely if the vertex os not the seam of the edge. */ + /* If the node is on a vertex we can convert the parameters safely if the vertex is not the seam of the edge. */ if (!cad_geometry->get_cad_manager ()->t8_geom_vertex_is_on_seam_edge (edge_node.entity_tag, edge_geometry_tag)) { reference_params.emplace (); From 0b7c5113ab934521c0fb46ec4bea65b6814eb5ca Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Wed, 4 Feb 2026 14:21:43 +0100 Subject: [PATCH 22/22] Update src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx Co-authored-by: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> --- src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx index 6ce0a320ee..aaa39c5e92 100644 --- a/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx +++ b/src/t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx @@ -1137,7 +1137,7 @@ t8_cmesh_process_tree_geometry (const t8_cmesh_t cmesh, const t8_eclass_t eclass reference_params = edge_node.parameters; break; case 1: - /* If the node is on a curve, we can use the parameters if the curve is not theseam of the surface. */ + /* If the node is on a curve, we can use the parameters if the curve is not the seam of the surface. */ if (!cad_geometry->get_cad_manager ()->t8_geom_edge_is_seam (edge_node.entity_tag, edge_geometry_tag)) { reference_params.emplace (); cad_geometry->get_cad_manager ()->t8_geom_edge_parameter_to_face_parameters (