From: Sascha Hofstetter Date: Sun, 29 Jun 2025 20:12:14 +0000 (+0200) Subject: dealii_tria_to_cgal_polygon implementation now takes mapping into account X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=603631cf732393322556154ec079a919743fa895;p=dealii.git dealii_tria_to_cgal_polygon implementation now takes mapping into account --- diff --git a/source/cgal/polygon.cc b/source/cgal/polygon.cc index c450f1fd68..6e5cfeabe6 100644 --- a/source/cgal/polygon.cc +++ b/source/cgal/polygon.cc @@ -84,35 +84,57 @@ namespace CGALWrappers // Counterclockwise first vertex index on first position, // counterclockwise second vertex index on second position. std::map face_vertex_indices; + std::map> vertex_to_point; // Iterate over all active cells at the boundary for (const auto &cell : tria.active_cell_iterators()) { - for (const unsigned int i : cell->face_indices()) + for (const unsigned int f : cell->face_indices()) { - const typename Triangulation<2, 2>::face_iterator &face = - cell->face(i); - // Ensure that vertex indices of the face are in - // counterclockwise order inserted in the map. - if (face->at_boundary()) + if (cell->face(f)->at_boundary()) { - if (cell->reference_cell() == ReferenceCells::Quadrilateral && - (i == 0 || i == 3)) + // get mapped vertices of the cell + const auto v_mapped = mapping.get_vertices(cell); + const unsigned int v0 = cell->face(f)->vertex_index(0); + const unsigned int v1 = cell->face(f)->vertex_index(1); + + if (cell->reference_cell() == ReferenceCells::Triangle) + { + // add indices and first mapped vertex of the face + vertex_to_point[v0] = v_mapped[f]; + face_vertex_indices[v0] = v1; + } + else if (cell->reference_cell() == + ReferenceCells::Quadrilateral) { - face_vertex_indices[face->vertex_index(1)] = - face->vertex_index(0); + // Ensure that vertex indices of the face are in + // counterclockwise order inserted in the map. + if (f == 0 || f == 3) + { + // add indices and first mapped vertex of the face + vertex_to_point[v1] = + v_mapped[GeometryInfo<2>::face_to_cell_vertices(f, + 1)]; + face_vertex_indices[v1] = v0; + } + else + { + // add indices and first mapped vertex of the face + vertex_to_point[v0] = + v_mapped[GeometryInfo<2>::face_to_cell_vertices(f, + 0)]; + face_vertex_indices[v0] = v1; + } } else { - face_vertex_indices[face->vertex_index(0)] = - face->vertex_index(1); + DEAL_II_ASSERT_UNREACHABLE(); } } } } CGAL::Polygon_2 polygon; - const auto &vertices = tria.get_vertices(); // Vertex to start counterclockwise insertion const unsigned int start_index = face_vertex_indices.begin()->first; @@ -122,12 +144,16 @@ namespace CGALWrappers // find next vertex index in counterclockwise order while (face_vertex_indices.size() > 0) { + const auto vertex_it = vertex_to_point.find(current_index); + Assert(vertex_it != vertex_to_point.end(), + ExcMessage("This should not occur, please report bug")); + polygon.push_back( dealii_point_to_cgal_point, 2>( - vertices[current_index])); + vertex_it->second)); + vertex_to_point.erase(vertex_it); const auto it = face_vertex_indices.find(current_index); - // If the boundary is one closed loop, the next vertex index // must exist as key until the map is empty. Assert(it != face_vertex_indices.end(), @@ -136,7 +162,9 @@ namespace CGALWrappers current_index = it->second; face_vertex_indices.erase(it); - // Ensure that last vertex index is the start index + // Ensure that last vertex index is the start index. + // This condition can be used to extend the code for + // trianglations with holes Assert( !(face_vertex_indices.size() == 0 && current_index != start_index), ExcMessage(