* 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.
*
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
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."));
// 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
&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(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff