From: Daniel Arndt Date: Sat, 7 Jul 2018 23:30:54 +0000 (+0200) Subject: Allow copying manifold ids in GridGenerator::merge_triangulations X-Git-Tag: v9.1.0-rc1~916^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a94e811b1be4cbaa7e28665a74ac2643bfae663f;p=dealii.git Allow copying manifold ids in GridGenerator::merge_triangulations --- diff --git a/doc/news/changes/minor/20180708DanielArndt b/doc/news/changes/minor/20180708DanielArndt new file mode 100644 index 0000000000..af38198ba8 --- /dev/null +++ b/doc/news/changes/minor/20180708DanielArndt @@ -0,0 +1,3 @@ +New: GridGenerator_merge_triangulations gained an option to copy manifold ids. +
+(Daniel Arndt, 2018/07/08) diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index e9f7e9f214..b801236da0 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -984,11 +984,11 @@ namespace GridGenerator * refined cells. * * @note The function copies the material ids of the cells of the two input - * triangulations into the output triangulation but it currently makes no - * attempt to do the same for boundary ids. In other words, if the two - * coarse meshes have anything but the default boundary indicators, then you - * will currently have to set boundary indicators again by hand in the - * output triangulation. + * triangulations into the output triangulation. If @p copy_manifold_ids is + * activated, manifold ids will be copied. Currently, boundary_ids are never + * copied. In other words, if the two coarse meshes have anything but the + * default boundary indicators, then you will currently have to set boundary + * indicators again by hand in the output triangulation. * * @note Unlike most GridGenerator functions, this function does not attach * any manifolds to @p result, nor does it set any manifold ids. @@ -1002,7 +1002,8 @@ namespace GridGenerator merge_triangulations(const Triangulation &triangulation_1, const Triangulation &triangulation_2, Triangulation & result, - const double duplicated_vertex_tolerance = 1.0e-12); + const double duplicated_vertex_tolerance = 1.0e-12, + const bool copy_manifold_ids = false); /** * Given the two triangulations specified as the first two arguments, create diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 94f28d20fa..3a55b2cf04 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -3793,7 +3793,8 @@ namespace GridGenerator merge_triangulations(const Triangulation &triangulation_1, const Triangulation &triangulation_2, Triangulation & result, - const double duplicated_vertex_tolerance) + const double duplicated_vertex_tolerance, + const bool copy_manifold_ids) { Assert(triangulation_1.n_levels() == 1, ExcMessage("The input triangulations must be coarse meshes.")); @@ -3809,36 +3810,121 @@ namespace GridGenerator // now form the union of the set of cells std::vector> cells; cells.reserve(triangulation_1.n_cells() + triangulation_2.n_cells()); - for (typename Triangulation::cell_iterator cell = - triangulation_1.begin(); - cell != triangulation_1.end(); - ++cell) + for (const auto &cell : triangulation_1.cell_iterators()) { CellData this_cell; for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) this_cell.vertices[v] = cell->vertex_index(v); this_cell.material_id = cell->material_id(); - cells.push_back(this_cell); + this_cell.manifold_id = cell->manifold_id(); + cells.push_back(std::move(this_cell)); } // now do the same for the other other mesh. note that we have to // translate the vertex indices - for (typename Triangulation::cell_iterator cell = - triangulation_2.begin(); - cell != triangulation_2.end(); - ++cell) + for (const auto &cell : triangulation_2.cell_iterators()) { CellData this_cell; for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) this_cell.vertices[v] = cell->vertex_index(v) + triangulation_1.n_vertices(); this_cell.material_id = cell->material_id(); - cells.push_back(this_cell); + this_cell.manifold_id = cell->manifold_id(); + cells.push_back(std::move(this_cell)); + } + + // Now for the SubCellData + SubCellData subcell_data; + if (copy_manifold_ids) + { + if (dim == 2) + { + subcell_data.boundary_lines.reserve(triangulation_1.n_lines() + + triangulation_2.n_lines()); + + auto copy_line_manifold_ids = + [&](const Triangulation &tria, + const unsigned int offset) { + for (typename Triangulation::line_iterator line = + tria.begin_face(); + line != tria.end_face(); + ++line) + if (line->manifold_id() != numbers::flat_manifold_id) + { + CellData<1> boundary_line; + boundary_line.vertices[0] = + line->vertex_index(0) + offset; + boundary_line.vertices[1] = + line->vertex_index(1) + offset; + boundary_line.manifold_id = line->manifold_id(); + subcell_data.boundary_lines.push_back( + std::move(boundary_line)); + } + }; + + copy_line_manifold_ids(triangulation_1, 0); + copy_line_manifold_ids(triangulation_2, + triangulation_1.n_vertices()); + } + + if (dim == 3) + { + subcell_data.boundary_quads.reserve(triangulation_1.n_quads() + + triangulation_2.n_quads()); + // we can't do better here than to loop over all the lines bounding + // a face. For regular meshes an (interior) line in 3D is part of + // four cells. So this should be an appropriate guess. + subcell_data.boundary_lines.reserve(triangulation_1.n_cells() * 4 + + triangulation_2.n_cells() * 4); + + auto copy_face_and_line_manifold_ids = + [&](const Triangulation &tria, + const unsigned int offset) { + for (typename Triangulation::face_iterator face = + tria.begin_face(); + face != tria.end_face(); + ++face) + if (face->manifold_id() != numbers::flat_manifold_id) + { + CellData<2> boundary_quad; + boundary_quad.vertices[0] = + face->vertex_index(0) + offset; + boundary_quad.vertices[1] = + face->vertex_index(1) + offset; + boundary_quad.vertices[2] = + face->vertex_index(2) + offset; + boundary_quad.vertices[3] = + face->vertex_index(3) + offset; + + boundary_quad.manifold_id = face->manifold_id(); + subcell_data.boundary_quads.push_back( + std::move(boundary_quad)); + } + for (const auto &cell : tria.cell_iterators()) + for (unsigned int l = 0; l < 12; ++l) + if (cell->line(l)->manifold_id() != + numbers::flat_manifold_id) + { + CellData<1> boundary_line; + boundary_line.vertices[0] = + cell->line(l)->vertex_index(0) + offset; + boundary_line.vertices[1] = + cell->line(l)->vertex_index(1) + offset; + boundary_line.manifold_id = + cell->line(l)->manifold_id(); + subcell_data.boundary_lines.push_back( + std::move(boundary_line)); + } + }; + + copy_face_and_line_manifold_ids(triangulation_1, 0); + copy_face_and_line_manifold_ids(triangulation_2, + triangulation_1.n_vertices()); + } } // throw out duplicated vertices from the two meshes, reorder vertices as // necessary and create the triangulation - SubCellData subcell_data; std::vector considered_vertices; GridTools::delete_duplicated_vertices(vertices, cells, diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in index f581285c88..b64c8f60a4 100644 --- a/source/grid/grid_generator.inst.in +++ b/source/grid/grid_generator.inst.in @@ -76,7 +76,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const Triangulation &triangulation_2, Triangulation &result, - const double duplicated_vertex_tolerance); + const double duplicated_vertex_tolerance, + const bool copy_manifold_ids); template void create_union_triangulation( diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index a62d82612c..188aafd91e 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -586,16 +586,19 @@ namespace GridTools } } - // now we got a renumbering list. simply - // renumber all vertices (non-duplicate - // vertices get renumbered to themselves, so - // nothing bad happens). after that, the - // duplicate vertices will be unused, so call - // delete_unused_vertices() to do that part - // of the job. - for (unsigned int c = 0; c < cells.size(); ++c) - for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v) - cells[c].vertices[v] = new_vertex_numbers[cells[c].vertices[v]]; + // now we got a renumbering list. simply renumber all vertices + // (non-duplicate vertices get renumbered to themselves, so nothing bad + // happens). after that, the duplicate vertices will be unused, so call + // delete_unused_vertices() to do that part of the job. + for (auto &cell : cells) + for (auto &vertex_index : cell.vertices) + vertex_index = new_vertex_numbers[vertex_index]; + for (auto &quad : subcelldata.boundary_quads) + for (auto &vertex_index : quad.vertices) + vertex_index = new_vertex_numbers[vertex_index]; + for (auto &line : subcelldata.boundary_lines) + for (auto &vertex_index : line.vertices) + vertex_index = new_vertex_numbers[vertex_index]; delete_unused_vertices(vertices, cells, subcelldata); } diff --git a/source/grid/tria.cc b/source/grid/tria.cc index c9b2e2bfd0..90a1bf268f 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -3106,28 +3106,25 @@ namespace internal ExcLineInexistant(line_vertices.first, line_vertices.second)); } - // Assert that only exterior - // lines are given a boundary - // indicator - AssertThrow( - line->at_boundary(), - ExcInteriorLineCantBeBoundary(line->vertex_index(0), - line->vertex_index(1), - boundary_line->boundary_id)); - - // and make sure that we don't - // attempt to reset the - // boundary indicator to a - // different than the - // previously set value - if (line->boundary_id() != 0) - AssertThrow(line->boundary_id() == boundary_line->boundary_id, - ExcMessage( - "Duplicate boundary lines are only allowed " - "if they carry the same boundary indicator.")); + // Only exterior lines can be given a boundary indicator + if (line->at_boundary()) + { + // make sure that we don't attempt to reset the boundary + // indicator to a different than the previously set value + if (line->boundary_id() != 0) + AssertThrow(line->boundary_id() == boundary_line->boundary_id, + ExcMessage( + "Duplicate boundary lines are only allowed " + "if they carry the same boundary indicator.")); - line->set_boundary_id_internal(boundary_line->boundary_id); + line->set_boundary_id_internal(boundary_line->boundary_id); + } // Set manifold id if given + if (line->manifold_id() != numbers::flat_manifold_id) + AssertThrow(line->manifold_id() == boundary_line->manifold_id, + ExcMessage( + "Duplicate boundary lines are only allowed " + "if they carry the same manifold indicator.")); line->set_manifold_id(boundary_line->manifold_id); } @@ -3262,26 +3259,25 @@ namespace internal // check whether this face is // really an exterior one - AssertThrow( - quad->at_boundary(), - ExcInteriorQuadCantBeBoundary(quad->vertex_index(0), - quad->vertex_index(1), - quad->vertex_index(2), - quad->vertex_index(3), - boundary_quad->boundary_id)); - - // and make sure that we don't - // attempt to reset the - // boundary indicator to a - // different than the - // previously set value - if (quad->boundary_id() != 0) - AssertThrow(quad->boundary_id() == boundary_quad->boundary_id, + if (quad->at_boundary()) + { + // and make sure that we don't attempt to reset the boundary + // indicator to a different than the previously set value + if (quad->boundary_id() != 0) + AssertThrow(quad->boundary_id() == boundary_quad->boundary_id, + ExcMessage( + "Duplicate boundary quads are only allowed " + "if they carry the same boundary indicator.")); + + quad->set_boundary_id_internal(boundary_quad->boundary_id); + } + // Set manifold id if given + if (quad->manifold_id() != numbers::flat_manifold_id) + AssertThrow(quad->manifold_id() == boundary_quad->manifold_id, ExcMessage( "Duplicate boundary quads are only allowed " - "if they carry the same boundary indicator.")); + "if they carry the same manifold indicator.")); - quad->set_boundary_id_internal(boundary_quad->boundary_id); quad->set_manifold_id(boundary_quad->manifold_id); }