if (tria.n_cells() != 3)
AssertThrow(false, ExcNotImplemented());
- double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
- double eps = 1e-3 * middle;
- Triangulation<3>::cell_iterator cell = tria.begin();
+ const double middle_radius =
+ (outer_radius - inner_radius) / 2e0 + inner_radius;
+ const double eps = 1e-3 * middle_radius;
- for (; cell != tria.end(); ++cell)
- for (const unsigned int f : GeometryInfo<3>::face_indices())
+ for (const auto &cell : tria.cell_iterators())
+ for (const unsigned int f : cell->face_indices())
{
- if (!cell->face(f)->at_boundary())
+ const auto face = cell->face(f);
+ if (!face->at_boundary())
continue;
- double radius = cell->face(f)->center().norm() - center.norm();
- if (std::fabs(cell->face(f)->center()[0]) <
- eps) // x = 0 set boundary 2
+ const double face_center_radius =
+ (face->center(true) - center).norm();
+
+ if (std::fabs(face->center()[0]) < eps) // x = 0 set boundary 2
{
- cell->face(f)->set_boundary_id(2);
- for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
- ++j)
- if (cell->face(f)->line(j)->at_boundary())
- if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
- cell->face(f)->line(j)->vertex(1).norm()) >
- eps)
- cell->face(f)->line(j)->set_boundary_id(2);
+ face->set_boundary_id(2);
+ for (const unsigned int line_no : face->line_indices())
+ if (face->line(line_no)->at_boundary())
+ if (std::fabs(face->line(line_no)->vertex(0).norm() -
+ face->line(line_no)->vertex(1).norm()) > eps)
+ face->line(line_no)->set_boundary_id(2);
}
- else if (std::fabs(cell->face(f)->center()[1]) <
- eps) // y = 0 set boundary 3
+ else if (std::fabs(face->center()[1]) < eps) // y = 0 set boundary 3
{
- cell->face(f)->set_boundary_id(3);
- for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
- ++j)
- if (cell->face(f)->line(j)->at_boundary())
- if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
- cell->face(f)->line(j)->vertex(1).norm()) >
- eps)
- cell->face(f)->line(j)->set_boundary_id(3);
+ face->set_boundary_id(3);
+ for (const unsigned int line_no : face->line_indices())
+ if (face->line(line_no)->at_boundary())
+ if (std::fabs(face->line(line_no)->vertex(0).norm() -
+ face->line(line_no)->vertex(1).norm()) > eps)
+ face->line(line_no)->set_boundary_id(3);
}
- else if (std::fabs(cell->face(f)->center()[2]) <
- eps) // z = 0 set boundary 4
+ else if (std::fabs(face->center()[2]) < eps) // z = 0 set boundary 4
{
- cell->face(f)->set_boundary_id(4);
- for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
- ++j)
- if (cell->face(f)->line(j)->at_boundary())
- if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
- cell->face(f)->line(j)->vertex(1).norm()) >
- eps)
- cell->face(f)->line(j)->set_boundary_id(4);
+ face->set_boundary_id(4);
+ for (const unsigned int line_no : face->line_indices())
+ if (face->line(line_no)->at_boundary())
+ if (std::fabs(face->line(line_no)->vertex(0).norm() -
+ face->line(line_no)->vertex(1).norm()) > eps)
+ face->line(line_no)->set_boundary_id(4);
}
- else if (radius < middle) // inner radius set boundary 0
+ else if (face_center_radius <
+ middle_radius) // inner radius set boundary 0
{
- cell->face(f)->set_boundary_id(0);
- for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
- ++j)
- if (cell->face(f)->line(j)->at_boundary())
- if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
- cell->face(f)->line(j)->vertex(1).norm()) <
- eps)
- cell->face(f)->line(j)->set_boundary_id(0);
+ face->set_boundary_id(0);
+ for (const unsigned int line_no : face->line_indices())
+ if (face->line(line_no)->at_boundary())
+ if (std::fabs(face->line(line_no)->vertex(0).norm() -
+ face->line(line_no)->vertex(1).norm()) < eps)
+ face->line(line_no)->set_boundary_id(0);
}
- else if (radius > middle) // outer radius set boundary 1
+ else if (face_center_radius >
+ middle_radius) // outer radius set boundary 1
{
- cell->face(f)->set_boundary_id(1);
- for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
- ++j)
- if (cell->face(f)->line(j)->at_boundary())
- if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
- cell->face(f)->line(j)->vertex(1).norm()) <
- eps)
- cell->face(f)->line(j)->set_boundary_id(1);
+ face->set_boundary_id(1);
+ for (const unsigned int line_no : face->line_indices())
+ if (face->line(line_no)->at_boundary())
+ if (std::fabs(face->line(line_no)->vertex(0).norm() -
+ face->line(line_no)->vertex(1).norm()) < eps)
+ face->line(line_no)->set_boundary_id(1);
}
else
DEAL_II_ASSERT_UNREACHABLE();
AssertThrow(false, ExcNotImplemented());
}
- if (colorize)
- colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
-
tria.set_all_manifold_ids(0);
tria.set_manifold(0, SphericalManifold<3>(center));
+
+ if (colorize)
+ colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
}
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2013 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// Test if colorizing a thin quarter_hyper_shell yields correct
+// boundary indicators. The function used to not colorize
+// boundaries correctly if the hyper shell was too thin.
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+void
+test(const Point<dim> origin = Point<dim>())
+{
+ SphericalManifold<dim> manifold(origin);
+ {
+ deallog << "quarter_hyper_shell with origin: " << origin << std::endl;
+ Triangulation<dim> tr;
+ GridGenerator::quarter_hyper_shell(tr, origin, 1460800, 1560800, 0, true);
+ tr.set_manifold(0, manifold);
+
+ bool found_top_boundary = false;
+ const types::boundary_id top_id = 1;
+ deallog << std::endl << "Boundary ids of boundary faces:";
+
+ for (const auto &cell : tr.active_cell_iterators())
+ for (const unsigned int face_no : cell->face_indices())
+ {
+ const typename Triangulation<dim>::face_iterator face =
+ cell->face(face_no);
+
+ if (face->at_boundary())
+ {
+ const types::boundary_id boundary_id = face->boundary_id();
+ deallog << " " << boundary_id;
+
+ if (boundary_id == top_id)
+ {
+ found_top_boundary = true;
+ }
+ }
+ }
+
+ if (dim > 2)
+ {
+ deallog << std::endl << std::endl << "Boundary ids of boundary lines:";
+
+ for (const auto &cell : tr.active_cell_iterators())
+ for (const unsigned int face_no : cell->face_indices())
+ {
+ const typename Triangulation<dim>::face_iterator face =
+ cell->face(face_no);
+
+ if (face->at_boundary())
+ {
+ for (const unsigned l : face->line_indices())
+ deallog << " " << face->line(l)->boundary_id();
+ }
+ }
+ }
+
+ deallog << std::endl;
+
+ if (found_top_boundary)
+ deallog << "Ok. Found top boundary." << std::endl;
+ else
+ deallog << "Error. Failed to find top boundary." << std::endl;
+ }
+}
+
+
+int
+main()
+{
+ initlog();
+
+ deallog.get_file_stream() << std::setprecision(9);
+
+ deallog.push("2d");
+ test<2>();
+ test<2>(Point<2>(1.0, 0.0));
+ deallog.pop();
+
+ deallog.push("3d");
+ test<3>();
+ test<3>(Point<3>(1.0, 0.0, 0.0));
+
+ deallog.pop();
+}