From: Peter Munch Date: Fri, 10 Jul 2020 09:19:14 +0000 (+0200) Subject: Add Simplex::GridGenerator X-Git-Tag: v9.3.0-rc1~1302^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10688%2Fhead;p=dealii.git Add Simplex::GridGenerator --- diff --git a/include/deal.II/simplex/grid_generator.h b/include/deal.II/simplex/grid_generator.h new file mode 100644 index 0000000000..f35f98ebf7 --- /dev/null +++ b/include/deal.II/simplex/grid_generator.h @@ -0,0 +1,248 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_simplex_grid_generator_h +#define dealii_simplex_grid_generator_h + + +#include + +#include + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + /** + * This namespace provides a collection of functions to generate simplex + * triangulations for some basic geometries. + */ + namespace GridGenerator + { + /** + * Create a coordinate-parallel brick from the two diagonally opposite + * corner points @p p1 and @p p2. The number of vertices in coordinate + * direction @p i is given by repetitions[i]+1. + * + * @note This function connects internally 4/8 vertices to quadrilateral/ + * hexahedral cells and subdivides these into 2/5 triangular/ + * tetrahedral cells. + * + * @note Currently, this function only works for `dim==spacedim`. + */ + template + void + subdivided_hyper_rectangle(Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize = false) + { + AssertDimension(dim, spacedim); + + AssertThrow(colorize == false, ExcNotImplemented()); + + std::vector> vertices; + std::vector> cells; + + if (dim == 2) + { + // determine cell sizes + Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1]); + + // create vertices + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back( + Point(p1[0] + dx[0] * i, p1[1] + dx[1] * j)); + + // create cells + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference QUAD cell + std::array quad{ + (j + 0) * (repetitions[0] + 1) + i + 0, // + (j + 0) * (repetitions[0] + 1) + i + 1, // + (j + 1) * (repetitions[0] + 1) + i + 0, // + (j + 1) * (repetitions[0] + 1) + i + 1 // + }; // + + // TRI cell 0 + { + CellData tri; + tri.vertices = {quad[0], quad[1], quad[2]}; + cells.push_back(tri); + } + + // TRI cell 1 + { + CellData tri; + tri.vertices = {quad[3], quad[2], quad[1]}; + cells.push_back(tri); + } + } + } + else if (dim == 3) + { + // determine cell sizes + Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1], + (p2[2] - p1[2]) / repetitions[1]); + + // create vertices + for (unsigned int k = 0; k <= repetitions[2]; ++k) + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back(Point(p1[0] + dx[0] * i, + p1[1] + dx[1] * j, + p1[2] + dx[2] * k)); + + // create cells + for (unsigned int k = 0; k < repetitions[2]; ++k) + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference HEX cell + std::array quad{ + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 0) * (repetitions[0] + 1) + i + 1, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 0, + (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) + + (j + 1) * (repetitions[0] + 1) + i + 1}; + + // TET cell 0 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {quad[0], quad[1], quad[2], quad[4]}; + else + cell.vertices = {quad[0], quad[1], quad[3], quad[5]}; + + cells.push_back(cell); + } + + // TET cell 1 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {quad[2], quad[1], quad[3], quad[7]}; + else + cell.vertices = {quad[0], quad[3], quad[2], quad[6]}; + cells.push_back(cell); + } + + // TET cell 2 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {quad[1], quad[4], quad[5], quad[7]}; + else + cell.vertices = {quad[0], quad[4], quad[5], quad[6]}; + cells.push_back(cell); + } + + // TET cell 3 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {quad[2], quad[4], quad[7], quad[6]}; + else + cell.vertices = {quad[3], quad[5], quad[7], quad[6]}; + cells.push_back(cell); + } + + // TET cell 4 + { + CellData cell; + if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0) + cell.vertices = {quad[1], quad[2], quad[4], quad[7]}; + else + cell.vertices = {quad[0], quad[3], quad[6], quad[5]}; + cells.push_back(cell); + } + } + } + else + { + AssertThrow(colorize == false, ExcNotImplemented()); + } + + // actually create triangulation + tria.create_triangulation(vertices, cells, SubCellData()); + } + + /** + * Initialize the given triangulation with a hypercube (square in 2D and + * cube in 3D) consisting of @p repetitions cells in each direction. + * The hypercube volume is the tensor product interval + * $[left,right]^{\text{dim}}$ in the present number of dimensions, where + * the limits are given as arguments. They default to zero and unity, then + * producing the unit hypercube. + * + * @note This function connects internally 4/8 vertices to quadrilateral/ + * hexahedral cells and subdivides these into 2/5 triangular/ + * tetrahedral cells. + */ + template + void + subdivided_hyper_cube(Triangulation &tria, + const unsigned int repetitions, + const double p1 = 0.0, + const double p2 = 1.0, + const bool colorize = false) + { + if (dim == 2) + { + subdivided_hyper_rectangle( + tria, {repetitions, repetitions}, {p1, p1}, {p2, p2}, colorize); + } + else if (dim == 3) + { + subdivided_hyper_rectangle(tria, + {repetitions, repetitions, repetitions}, + {p1, p1, p1}, + {p2, p2, p2}, + colorize); + } + else + { + AssertThrow(false, ExcNotImplemented()) + } + } + } // namespace GridGenerator +} // namespace Simplex + + + +DEAL_II_NAMESPACE_CLOSE + +#endif