From e315e3b7972d4226457025e698bb7f509ea6defc Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 24 May 2018 14:36:03 -0400 Subject: [PATCH] Add a tolerance for merge_triangulations. The default setting is not robust with respect to roundoff errors (or other minor numerical errors). --- include/deal.II/grid/grid_generator.h | 27 +++- source/grid/grid_generator.cc | 10 +- source/grid/grid_generator.inst.in | 3 +- tests/grid/merge_triangulations_04.cc | 168 ++++++++++++++++++++++ tests/grid/merge_triangulations_04.output | 122 ++++++++++++++++ 5 files changed, 325 insertions(+), 5 deletions(-) create mode 100644 tests/grid/merge_triangulations_04.cc create mode 100644 tests/grid/merge_triangulations_04.output diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 284d63e49b..02e1139e47 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -975,6 +975,30 @@ namespace GridGenerator * mesh building blocks are GridTools::transform, GridTools::rotate, and * GridTools::scale). * + * Vertices that are less than @p duplicated_vertex_tolerance apart will be merged + * together. It is usually necessary to set this value to something that + * depends on the input triangulations in some way. One reasonable choice is + * to use the minimum distance between all adjacent vertices of the input + * mesh divided by some constant: + * + * @code + * auto min_line_length = [](const Triangulation &tria) -> double + * { + * double length = std::numeric_limits::max(); + * for (const auto cell : tria.active_cell_iterators()) + * for (unsigned int n = 0; n < GeometryInfo::lines_per_cell; ++n) + * length = std::min(length, (cell->line(n)->vertex(0) - + * cell->line(n)->vertex(1)).norm()); + * return length; + * }; + * + * const double tolerance = std::min(min_line_length(triangulation_1), + * min_line_length(triangulation_2)) / 2.0; + * @endcode + * + * This will merge any vertices that are closer than any pair of vertices on + * the input meshes. + * * @note The two input triangulations must be coarse meshes that have no * refined cells. * @@ -996,7 +1020,8 @@ namespace GridGenerator void merge_triangulations(const Triangulation &triangulation_1, const Triangulation &triangulation_2, - Triangulation & result); + Triangulation & result, + const double duplicated_vertex_tolerance = 1.0e-12); /** * 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 f7fa27710c..29df390576 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -3784,7 +3784,8 @@ namespace GridGenerator void merge_triangulations(const Triangulation &triangulation_1, const Triangulation &triangulation_2, - Triangulation & result) + Triangulation & result, + const double duplicated_vertex_tolerance) { Assert(triangulation_1.n_levels() == 1, ExcMessage("The input triangulations must be coarse meshes.")); @@ -3831,8 +3832,11 @@ namespace GridGenerator // necessary and create the triangulation SubCellData subcell_data; std::vector considered_vertices; - GridTools::delete_duplicated_vertices( - vertices, cells, subcell_data, considered_vertices); + GridTools::delete_duplicated_vertices(vertices, + cells, + subcell_data, + considered_vertices, + duplicated_vertex_tolerance); // reorder the cells to ensure that they satisfy the convention for // edge and face directions diff --git a/source/grid/grid_generator.inst.in b/source/grid/grid_generator.inst.in index 4fd6e2958f..136480afb9 100644 --- a/source/grid/grid_generator.inst.in +++ b/source/grid/grid_generator.inst.in @@ -75,7 +75,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) &triangulation_1, const Triangulation &triangulation_2, - Triangulation &result); + Triangulation &result, + const double duplicated_vertex_tolerance); template void create_union_triangulation( diff --git a/tests/grid/merge_triangulations_04.cc b/tests/grid/merge_triangulations_04.cc new file mode 100644 index 0000000000..04cc6ef1c6 --- /dev/null +++ b/tests/grid/merge_triangulations_04.cc @@ -0,0 +1,168 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +int +main() +{ + initlog(); + + // set up the bulk triangulation + Triangulation<2> bulk_triangulation; + GridGenerator::subdivided_hyper_rectangle( + bulk_triangulation, {22u, 4u}, Point<2>(0.0, 0.0), Point<2>(2.2, 0.41)); + std::set::active_cell_iterator> cells_to_remove; + Tensor<1, 2> cylinder_triangulation_offset; + for (const auto cell : bulk_triangulation.active_cell_iterators()) + { + if ((cell->center() - Point<2>(0.2, 0.2)).norm() < 0.15) + cells_to_remove.insert(cell); + + if (cylinder_triangulation_offset == Point<2>()) + { + for (unsigned int vertex_n = 0; + vertex_n < GeometryInfo<2>::vertices_per_cell; + ++vertex_n) + if (cell->vertex(vertex_n) == Point<2>()) + { + // skip two cells in the bottom left corner + cylinder_triangulation_offset = + 2.0 * (cell->vertex(3) - Point<2>()); + break; + } + } + } + Triangulation<2> result_1; + GridGenerator::create_triangulation_with_removed_cells( + bulk_triangulation, cells_to_remove, result_1); + + // set up the cylinder triangulation + Triangulation<2> cylinder_triangulation; + GridGenerator::hyper_cube_with_cylindrical_hole( + cylinder_triangulation, 0.05, 0.41 / 4.0); + GridTools::shift(cylinder_triangulation_offset, cylinder_triangulation); + // dumb hack + for (const auto cell : cylinder_triangulation.active_cell_iterators()) + cell->set_material_id(2); + + // merge them together + auto minimal_line_length = [](const Triangulation<2> &tria) -> double { + double min_line_length = std::numeric_limits::max(); + for (const auto cell : tria.active_cell_iterators()) + for (unsigned int line_n = 0; line_n < GeometryInfo<2>::lines_per_cell; + ++line_n) + min_line_length = std::min( + min_line_length, + (cell->line(line_n)->vertex(0) - cell->line(line_n)->vertex(1)) + .norm()); + return min_line_length; + }; + + // the cylindrical triangulation might not match the Cartesian grid: as a + // result the vertices might not be lined up. Get around this by deleting + // duplicated vertices with a very low numerical tolerance. + const double tolerance = + std::min(minimal_line_length(result_1), + minimal_line_length(cylinder_triangulation)) / + 2.0; + + + Triangulation<2> result_2; + GridGenerator::merge_triangulations( + result_1, cylinder_triangulation, result_2, tolerance); + + const types::manifold_id tfi_id = 1; + const types::manifold_id polar_id = 0; + for (const auto cell : result_2.active_cell_iterators()) + { + // set all non-boundary manifold ids to the new TFI manifold id. + if (cell->material_id() == 2) + { + cell->set_manifold_id(tfi_id); + for (unsigned int face_n = 0; + face_n < GeometryInfo<2>::faces_per_cell; + ++face_n) + { + if (cell->face(face_n)->at_boundary()) + cell->face(face_n)->set_manifold_id(polar_id); + else + cell->face(face_n)->set_manifold_id(tfi_id); + } + } + } + + PolarManifold<2> polar_manifold(Point<2>(0.2, 0.2)); + result_2.set_manifold(polar_id, polar_manifold); + TransfiniteInterpolationManifold<2> inner_manifold; + inner_manifold.initialize(result_2); + result_2.set_manifold(tfi_id, inner_manifold); + + std::vector *> inner_pointers; + for (const auto cell : result_2.active_cell_iterators()) + { + for (unsigned int face_n = 0; face_n < GeometryInfo<2>::faces_per_cell; + ++face_n) + { + if (cell->face(face_n)->manifold_id() == polar_id) + { + inner_pointers.push_back(&cell->face(face_n)->vertex(0)); + inner_pointers.push_back(&cell->face(face_n)->vertex(1)); + } + } + } + // de-duplicate + std::sort(inner_pointers.begin(), inner_pointers.end()); + inner_pointers.erase( + std::unique(inner_pointers.begin(), inner_pointers.end()), + inner_pointers.end()); + + // find the current center... + Point<2> center; + for (const Point<2> *const ptr : inner_pointers) + center += *ptr / double(inner_pointers.size()); + + // and recenter at (0.2, 0.2) + for (Point<2> *const ptr : inner_pointers) + *ptr += Point<2>(0.2, 0.2) - center; + + result_2.refine_global(1); + for (const auto cell : result_2.active_cell_iterators()) + { + if (cell->at_boundary()) + { + for (unsigned int face_n = 0; + face_n < GeometryInfo<2>::faces_per_cell; + ++face_n) + { + auto face = cell->face(face_n); + if (face->at_boundary()) + { + deallog << "boundary face with center " << face->center() + << std::endl; + } + } + } + } + deallog << "OK" << std::endl; +} diff --git a/tests/grid/merge_triangulations_04.output b/tests/grid/merge_triangulations_04.output new file mode 100644 index 0000000000..ebb884e9cd --- /dev/null +++ b/tests/grid/merge_triangulations_04.output @@ -0,0 +1,122 @@ + +DEAL::boundary face with center 0.00000 0.0768750 +DEAL::boundary face with center 0.0750000 0.00000 +DEAL::boundary face with center 0.00000 0.0256250 +DEAL::boundary face with center 0.0250000 0.00000 +DEAL::boundary face with center 0.175000 0.00000 +DEAL::boundary face with center 0.125000 0.00000 +DEAL::boundary face with center 0.275000 0.00000 +DEAL::boundary face with center 0.225000 0.00000 +DEAL::boundary face with center 0.375000 0.00000 +DEAL::boundary face with center 0.325000 0.00000 +DEAL::boundary face with center 0.475000 0.00000 +DEAL::boundary face with center 0.425000 0.00000 +DEAL::boundary face with center 0.575000 0.00000 +DEAL::boundary face with center 0.525000 0.00000 +DEAL::boundary face with center 0.675000 0.00000 +DEAL::boundary face with center 0.625000 0.00000 +DEAL::boundary face with center 0.775000 0.00000 +DEAL::boundary face with center 0.725000 0.00000 +DEAL::boundary face with center 0.875000 0.00000 +DEAL::boundary face with center 0.825000 0.00000 +DEAL::boundary face with center 0.975000 0.00000 +DEAL::boundary face with center 0.925000 0.00000 +DEAL::boundary face with center 1.07500 0.00000 +DEAL::boundary face with center 1.02500 0.00000 +DEAL::boundary face with center 1.17500 0.00000 +DEAL::boundary face with center 1.12500 0.00000 +DEAL::boundary face with center 1.27500 0.00000 +DEAL::boundary face with center 1.22500 0.00000 +DEAL::boundary face with center 1.37500 0.00000 +DEAL::boundary face with center 1.32500 0.00000 +DEAL::boundary face with center 1.47500 0.00000 +DEAL::boundary face with center 1.42500 0.00000 +DEAL::boundary face with center 1.57500 0.00000 +DEAL::boundary face with center 1.52500 0.00000 +DEAL::boundary face with center 1.67500 0.00000 +DEAL::boundary face with center 1.62500 0.00000 +DEAL::boundary face with center 1.77500 0.00000 +DEAL::boundary face with center 1.72500 0.00000 +DEAL::boundary face with center 1.87500 0.00000 +DEAL::boundary face with center 1.82500 0.00000 +DEAL::boundary face with center 1.97500 0.00000 +DEAL::boundary face with center 1.92500 0.00000 +DEAL::boundary face with center 2.07500 0.00000 +DEAL::boundary face with center 2.02500 0.00000 +DEAL::boundary face with center 2.20000 0.0768750 +DEAL::boundary face with center 2.20000 0.0256250 +DEAL::boundary face with center 2.17500 0.00000 +DEAL::boundary face with center 2.12500 0.00000 +DEAL::boundary face with center 0.00000 0.179375 +DEAL::boundary face with center 0.00000 0.128125 +DEAL::boundary face with center 2.20000 0.179375 +DEAL::boundary face with center 2.20000 0.128125 +DEAL::boundary face with center 0.00000 0.281875 +DEAL::boundary face with center 0.00000 0.230625 +DEAL::boundary face with center 2.20000 0.281875 +DEAL::boundary face with center 2.20000 0.230625 +DEAL::boundary face with center 0.0750000 0.410000 +DEAL::boundary face with center 0.00000 0.384375 +DEAL::boundary face with center 0.0250000 0.410000 +DEAL::boundary face with center 0.00000 0.333125 +DEAL::boundary face with center 0.175000 0.410000 +DEAL::boundary face with center 0.125000 0.410000 +DEAL::boundary face with center 0.275000 0.410000 +DEAL::boundary face with center 0.225000 0.410000 +DEAL::boundary face with center 0.375000 0.410000 +DEAL::boundary face with center 0.325000 0.410000 +DEAL::boundary face with center 0.475000 0.410000 +DEAL::boundary face with center 0.425000 0.410000 +DEAL::boundary face with center 0.575000 0.410000 +DEAL::boundary face with center 0.525000 0.410000 +DEAL::boundary face with center 0.675000 0.410000 +DEAL::boundary face with center 0.625000 0.410000 +DEAL::boundary face with center 0.775000 0.410000 +DEAL::boundary face with center 0.725000 0.410000 +DEAL::boundary face with center 0.875000 0.410000 +DEAL::boundary face with center 0.825000 0.410000 +DEAL::boundary face with center 0.975000 0.410000 +DEAL::boundary face with center 0.925000 0.410000 +DEAL::boundary face with center 1.07500 0.410000 +DEAL::boundary face with center 1.02500 0.410000 +DEAL::boundary face with center 1.17500 0.410000 +DEAL::boundary face with center 1.12500 0.410000 +DEAL::boundary face with center 1.27500 0.410000 +DEAL::boundary face with center 1.22500 0.410000 +DEAL::boundary face with center 1.37500 0.410000 +DEAL::boundary face with center 1.32500 0.410000 +DEAL::boundary face with center 1.47500 0.410000 +DEAL::boundary face with center 1.42500 0.410000 +DEAL::boundary face with center 1.57500 0.410000 +DEAL::boundary face with center 1.52500 0.410000 +DEAL::boundary face with center 1.67500 0.410000 +DEAL::boundary face with center 1.62500 0.410000 +DEAL::boundary face with center 1.77500 0.410000 +DEAL::boundary face with center 1.72500 0.410000 +DEAL::boundary face with center 1.87500 0.410000 +DEAL::boundary face with center 1.82500 0.410000 +DEAL::boundary face with center 1.97500 0.410000 +DEAL::boundary face with center 1.92500 0.410000 +DEAL::boundary face with center 2.07500 0.410000 +DEAL::boundary face with center 2.02500 0.410000 +DEAL::boundary face with center 2.20000 0.384375 +DEAL::boundary face with center 2.17500 0.410000 +DEAL::boundary face with center 2.12500 0.410000 +DEAL::boundary face with center 2.20000 0.333125 +DEAL::boundary face with center 0.240775 0.227245 +DEAL::boundary face with center 0.248097 0.209567 +DEAL::boundary face with center 0.227245 0.240775 +DEAL::boundary face with center 0.209567 0.248097 +DEAL::boundary face with center 0.190433 0.248097 +DEAL::boundary face with center 0.172755 0.240775 +DEAL::boundary face with center 0.159225 0.227245 +DEAL::boundary face with center 0.151903 0.209567 +DEAL::boundary face with center 0.151903 0.190433 +DEAL::boundary face with center 0.159225 0.172755 +DEAL::boundary face with center 0.190433 0.151903 +DEAL::boundary face with center 0.172755 0.159225 +DEAL::boundary face with center 0.227245 0.159225 +DEAL::boundary face with center 0.209567 0.151903 +DEAL::boundary face with center 0.248097 0.190433 +DEAL::boundary face with center 0.240775 0.172755 +DEAL::OK -- 2.39.5