]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a tolerance for merge_triangulations. 6674/head
authorDavid Wells <wellsd2@rpi.edu>
Thu, 24 May 2018 18:36:03 +0000 (14:36 -0400)
committerDavid Wells <drwells@vt.edu>
Mon, 28 May 2018 13:46:55 +0000 (09:46 -0400)
The default setting is not robust with respect to roundoff errors (or
other minor numerical errors).

include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc
source/grid/grid_generator.inst.in
tests/grid/merge_triangulations_04.cc [new file with mode: 0644]
tests/grid/merge_triangulations_04.output [new file with mode: 0644]

index 284d63e49b90a89125bda3871ba7dc1c3c75a43e..02e1139e47783038fa4574eea9816c37b9c4525c 100644 (file)
@@ -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<dim> &tria) -> double
+   * {
+   *   double length = std::numeric_limits<double>::max();
+   *   for (const auto cell : tria.active_cell_iterators())
+   *     for (unsigned int n = 0; n < GeometryInfo<dim>::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<dim, spacedim> &triangulation_1,
                        const Triangulation<dim, spacedim> &triangulation_2,
-                       Triangulation<dim, spacedim> &      result);
+                       Triangulation<dim, spacedim> &      result,
+                       const double duplicated_vertex_tolerance = 1.0e-12);
 
   /**
    * Given the two triangulations specified as the first two arguments, create
index f7fa27710c0930b17160b5d8dd35b48301eaff8b..29df390576778f231b27f464636068234ac97401 100644 (file)
@@ -3784,7 +3784,8 @@ namespace GridGenerator
   void
   merge_triangulations(const Triangulation<dim, spacedim> &triangulation_1,
                        const Triangulation<dim, spacedim> &triangulation_2,
-                       Triangulation<dim, spacedim> &      result)
+                       Triangulation<dim, spacedim> &      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<unsigned int> 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
index 4fd6e2958ff64c58e01e3150ccc2010b5b107f36..136480afb96480c2a972c35f1e6454f82d95371e 100644 (file)
@@ -75,7 +75,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
           &triangulation_1,
         const Triangulation<deal_II_dimension, deal_II_space_dimension>
           &triangulation_2,
-        Triangulation<deal_II_dimension, deal_II_space_dimension> &result);
+        Triangulation<deal_II_dimension, deal_II_space_dimension> &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 (file)
index 0000000..04cc6ef
--- /dev/null
@@ -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 <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#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<Triangulation<2>::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<double>::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<Point<2> *> 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 (file)
index 0000000..ebb884e
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.