From 12909d1617d7d0caeb4264e3498c93ddbc14ec75 Mon Sep 17 00:00:00 2001 From: Bruno Date: Thu, 1 Apr 2021 14:44:23 -0400 Subject: [PATCH] Add a subdivided_cylinder grid generator and test --- include/deal.II/grid/grid_generator.h | 40 +++++ source/grid/grid_generator.cc | 150 ++++++++++++++++ source/grid/grid_generator_from_name.cc | 4 + .../grid_generator_subdivided_cylinder.cc | 46 +++++ .../grid_generator_subdivided_cylinder.output | 164 ++++++++++++++++++ 5 files changed, 404 insertions(+) create mode 100644 tests/grid/grid_generator_subdivided_cylinder.cc create mode 100644 tests/grid/grid_generator_subdivided_cylinder.output diff --git a/include/deal.II/grid/grid_generator.h b/include/deal.II/grid/grid_generator.h index 8e641c20ea..f5b8e90b62 100644 --- a/include/deal.II/grid/grid_generator.h +++ b/include/deal.II/grid/grid_generator.h @@ -948,6 +948,46 @@ namespace GridGenerator const double radius = 1., const double half_length = 1.); + + /** + * Create a @p dim dimensional cylinder where the $x$-axis serves as + * the axis of the cylinder. For the purposes of this function, a + * cylinder is defined as a (@p dim - 1) dimensional disk of given + * @p radius, extruded along the axis of the cylinder (which is the + * first coordinate direction). Consequently, in three dimensions, + * the cylinder extends from `x=-half_length` to `x=+half_length` + * and its projection into the @p yz-plane is a circle of radius @p + * radius. In two dimensions, the cylinder is a rectangle from + * `x=-half_length` to `x=+half_length` and from `y=-radius` to + * `y=radius`. + * + * The boundaries are colored according to the following scheme: 0 for the + * hull of the cylinder, 1 for the left hand face and 2 for the right hand + * face (see + * @ref GlossColorization "the glossary entry on colorization"). + * + * The manifold id for the hull of the cylinder is set to zero, and a + * CylindricalManifold is attached to it. + * + * @pre The triangulation passed as argument needs to be empty when calling + * this function. + * + * @param x_subdivisions A positive integer denoting the number + * of cells to generate in the x direction. The default cylinder has + * x_repetitions=2. + * + * @param radius The radius of the circle in the yz-plane used to extrude the cylinder. + * + * @param half_length The half-length of the cylinder in the x direction. + */ + template + void + subdivided_cylinder(Triangulation &tria, + const unsigned int x_subdivisions, + const double radius = 1., + const double half_length = 1.); + + /** * Create a cut cone around the x-axis. The cone extends from * x=-half_length to x=half_length and its projection into diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index c845f2d454..3f2fb6382e 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -3952,6 +3952,16 @@ namespace GridGenerator } + template <> + void subdivided_cylinder(Triangulation<1> &, + const unsigned int, + const double, + const double) + { + Assert(false, ExcNotImplemented()); + } + + template <> void @@ -4470,6 +4480,15 @@ namespace GridGenerator } } + template <> + void subdivided_cylinder(Triangulation<2> &, + const unsigned int, + const double, + const double) + { + Assert(false, ExcNotImplemented()); + } + // Implementation for 2D only @@ -5560,6 +5579,137 @@ namespace GridGenerator tria.set_manifold(0, CylindricalManifold<3>()); } + // Implementation for 3D only + template <> + void subdivided_cylinder(Triangulation<3> & tria, + const unsigned int x_subdivisions, + const double radius, + const double half_length) + { + // Copy the base from hyper_ball<3> + // and transform it to yz + const double d = radius / std::sqrt(2.0); + const double a = d / (1 + std::sqrt(2.0)); + + // const unsigned int npts_per_plane = 8; + std::vector> vertices; + const double initial_height = -half_length; + const double height_increment = 2. * half_length / x_subdivisions; + + for (unsigned int rep = 0; rep < (x_subdivisions + 1); ++rep) + { + const double height = initial_height + height_increment * rep; + + vertices.emplace_back(Point<3>(-d, height, -d)); + vertices.emplace_back(Point<3>(d, height, -d)); + vertices.emplace_back(Point<3>(-a, height, -a)); + vertices.emplace_back(Point<3>(a, height, -a)); + vertices.emplace_back(Point<3>(-a, height, a)); + vertices.emplace_back(Point<3>(a, height, a)); + vertices.emplace_back(Point<3>(-d, height, d)); + vertices.emplace_back(Point<3>(d, height, d)); + } + + // Turn cylinder such that y->x + for (auto &vertex : vertices) + { + const double h = vertex(1); + vertex(1) = -vertex(0); + vertex(0) = h; + } + + std::vector> cell_vertices; + cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11}); + cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12}); + cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13}); + cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13}); + cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13}); + + for (unsigned int rep = 1; rep < x_subdivisions; ++rep) + { + for (unsigned int i = 0; i < 5; ++i) + { + std::vector new_cell_vertices(8); + for (unsigned int j = 0; j < 8; ++j) + new_cell_vertices[j] = cell_vertices[i][j] + 8 * rep; + cell_vertices.push_back(new_cell_vertices); + } + } + + unsigned int n_cells = x_subdivisions * 5; + + std::vector> cells(n_cells, CellData<3>()); + + for (unsigned int i = 0; i < n_cells; ++i) + { + for (unsigned int j = 0; j < 8; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + tria.create_triangulation(std::vector>(std::begin(vertices), + std::end(vertices)), + cells, + SubCellData()); // no boundary information + + // set boundary indicators for the + // faces at the ends to 1 and 2, + // respectively. note that we also + // have to deal with those lines + // that are purely in the interior + // of the ends. we determine whether + // an edge is purely in the + // interior if one of its vertices + // is at coordinates '+-a' as set + // above + Triangulation<3>::cell_iterator cell = tria.begin(); + Triangulation<3>::cell_iterator end = tria.end(); + + tria.set_all_manifold_ids_on_boundary(0); + + for (; cell != end; ++cell) + for (unsigned int i : GeometryInfo<3>::face_indices()) + if (cell->at_boundary(i)) + { + if (cell->face(i)->center()(0) > half_length - 1.e-5) + { + cell->face(i)->set_boundary_id(2); + cell->face(i)->set_manifold_id(numbers::flat_manifold_id); + + for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face; + ++e) + if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) + { + cell->face(i)->line(e)->set_boundary_id(2); + cell->face(i)->line(e)->set_manifold_id( + numbers::flat_manifold_id); + } + } + else if (cell->face(i)->center()(0) < -half_length + 1.e-5) + { + cell->face(i)->set_boundary_id(1); + cell->face(i)->set_manifold_id(numbers::flat_manifold_id); + + for (unsigned int e = 0; e < GeometryInfo<3>::lines_per_face; + ++e) + if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) + { + cell->face(i)->line(e)->set_boundary_id(1); + cell->face(i)->line(e)->set_manifold_id( + numbers::flat_manifold_id); + } + } + } + tria.set_manifold(0, CylindricalManifold<3>()); + } + + template <> void quarter_hyper_ball(Triangulation<3> &tria, diff --git a/source/grid/grid_generator_from_name.cc b/source/grid/grid_generator_from_name.cc index d482394086..4fcb174f2d 100644 --- a/source/grid/grid_generator_from_name.cc +++ b/source/grid/grid_generator_from_name.cc @@ -142,6 +142,10 @@ namespace GridGenerator else if (name == "cylinder") parse_and_create(cylinder, arguments, tria); + else if (name == "subdivided_cylinder") + parse_and_create( + subdivided_cylinder, arguments, tria); + else if (name == "truncated_cone") parse_and_create(truncated_cone, arguments, diff --git a/tests/grid/grid_generator_subdivided_cylinder.cc b/tests/grid/grid_generator_subdivided_cylinder.cc new file mode 100644 index 0000000000..da9999beba --- /dev/null +++ b/tests/grid/grid_generator_subdivided_cylinder.cc @@ -0,0 +1,46 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2020 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Test GridReordering<3>::reorder_cells + +#include +#include +#include + +#include "../tests.h" + + + +void +dim_3(std::ostream &os) +{ + const unsigned int d = 3; + Triangulation tr; + + GridGenerator::subdivided_cylinder(tr, 4, 1, 1); + + GridOut gout; + gout.write_vtk(tr, os); +} + +int +main() +{ + initlog(true); + std::ostream &logfile = deallog.get_file_stream(); + dim_3(logfile); +} diff --git a/tests/grid/grid_generator_subdivided_cylinder.output b/tests/grid/grid_generator_subdivided_cylinder.output new file mode 100644 index 0000000000..662964b4eb --- /dev/null +++ b/tests/grid/grid_generator_subdivided_cylinder.output @@ -0,0 +1,164 @@ + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 40 double +-1.00000 0.707107 -0.707107 +-1.00000 -0.707107 -0.707107 +-1.00000 0.292893 -0.292893 +-1.00000 -0.292893 -0.292893 +-1.00000 0.292893 0.292893 +-1.00000 -0.292893 0.292893 +-1.00000 0.707107 0.707107 +-1.00000 -0.707107 0.707107 +-0.500000 0.707107 -0.707107 +-0.500000 -0.707107 -0.707107 +-0.500000 0.292893 -0.292893 +-0.500000 -0.292893 -0.292893 +-0.500000 0.292893 0.292893 +-0.500000 -0.292893 0.292893 +-0.500000 0.707107 0.707107 +-0.500000 -0.707107 0.707107 +0.00000 0.707107 -0.707107 +0.00000 -0.707107 -0.707107 +0.00000 0.292893 -0.292893 +0.00000 -0.292893 -0.292893 +0.00000 0.292893 0.292893 +0.00000 -0.292893 0.292893 +0.00000 0.707107 0.707107 +0.00000 -0.707107 0.707107 +0.500000 0.707107 -0.707107 +0.500000 -0.707107 -0.707107 +0.500000 0.292893 -0.292893 +0.500000 -0.292893 -0.292893 +0.500000 0.292893 0.292893 +0.500000 -0.292893 0.292893 +0.500000 0.707107 0.707107 +0.500000 -0.707107 0.707107 +1.00000 0.707107 -0.707107 +1.00000 -0.707107 -0.707107 +1.00000 0.292893 -0.292893 +1.00000 -0.292893 -0.292893 +1.00000 0.292893 0.292893 +1.00000 -0.292893 0.292893 +1.00000 0.707107 0.707107 +1.00000 -0.707107 0.707107 + +CELLS 98 466 +8 0 1 9 8 2 3 11 10 +8 0 2 10 8 6 4 12 14 +8 2 3 11 10 4 5 13 12 +8 1 7 15 9 3 5 13 11 +8 6 4 12 14 7 5 13 15 +8 8 9 17 16 10 11 19 18 +8 8 10 18 16 14 12 20 22 +8 10 11 19 18 12 13 21 20 +8 9 15 23 17 11 13 21 19 +8 14 12 20 22 15 13 21 23 +8 16 17 25 24 18 19 27 26 +8 16 18 26 24 22 20 28 30 +8 18 19 27 26 20 21 29 28 +8 17 23 31 25 19 21 29 27 +8 22 20 28 30 23 21 29 31 +8 24 25 33 32 26 27 35 34 +8 24 26 34 32 30 28 36 38 +8 26 27 35 34 28 29 37 36 +8 25 31 39 33 27 29 37 35 +8 30 28 36 38 31 29 37 39 +4 0 2 3 1 +4 0 6 4 2 +4 0 8 14 6 +4 0 1 9 8 +4 1 3 5 7 +4 1 7 15 9 +4 2 4 5 3 +4 6 7 5 4 +4 6 14 15 7 +4 8 16 22 14 +4 8 9 17 16 +4 9 15 23 17 +4 14 22 23 15 +4 16 24 30 22 +4 16 17 25 24 +4 17 23 31 25 +4 22 30 31 23 +4 24 32 38 30 +4 24 25 33 32 +4 25 31 39 33 +4 30 38 39 31 +4 32 34 35 33 +4 32 38 36 34 +4 33 35 37 39 +4 34 36 37 35 +4 38 39 37 36 +2 0 1 +2 2 3 +2 0 2 +2 1 3 +2 6 4 +2 0 6 +2 2 4 +2 0 8 +2 8 14 +2 6 14 +2 1 9 +2 8 9 +2 1 7 +2 3 5 +2 7 5 +2 7 15 +2 9 15 +2 4 5 +2 6 7 +2 14 15 +2 8 16 +2 16 22 +2 14 22 +2 9 17 +2 16 17 +2 15 23 +2 17 23 +2 22 23 +2 16 24 +2 24 30 +2 22 30 +2 17 25 +2 24 25 +2 23 31 +2 25 31 +2 30 31 +2 32 34 +2 24 32 +2 32 38 +2 30 38 +2 25 33 +2 32 33 +2 33 35 +2 31 39 +2 33 39 +2 34 36 +2 34 35 +2 35 37 +2 36 37 +2 38 39 +2 38 36 +2 39 37 + +CELL_TYPES 98 +12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 +9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 +3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + +CELL_DATA 98 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 1 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 +0 1 1 1 1 0 1 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 2 0 0 2 2 2 2 0 2 2 + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 +-1 -1 0 0 -1 0 -1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 -1 -1 -1 -1 +0 -1 -1 -1 -1 0 -1 0 0 0 0 0 0 -1 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 -1 0 0 -1 -1 -1 -1 0 -1 -1 -- 2.39.5