From 1f6d5a6aec2d5f5b325e77cfa4760d3f93c7f799 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 30 Jul 2020 19:57:41 +0200 Subject: [PATCH] Fix confilct with Simplex::GridGenerator --- include/deal.II/base/exceptions.h | 8 + include/deal.II/simplex/grid_generator.h | 254 ++++------------------- source/simplex/CMakeLists.txt | 2 + source/simplex/grid_generator.cc | 215 +++++++++++++++++++ source/simplex/grid_generator.inst.in | 37 ++++ 5 files changed, 306 insertions(+), 210 deletions(-) create mode 100644 source/simplex/grid_generator.cc create mode 100644 source/simplex/grid_generator.inst.in diff --git a/include/deal.II/base/exceptions.h b/include/deal.II/base/exceptions.h index b244dd3b70..171411eb07 100644 --- a/include/deal.II/base/exceptions.h +++ b/include/deal.II/base/exceptions.h @@ -1091,6 +1091,14 @@ namespace StandardExceptions "You are attempting to use functionality that is only available " "if deal.II was configured to use MPI."); + /** + * This function requires simplex support. + */ + DeclExceptionMsg( + ExcNeedsSimplexSupport, + "You are attempting to use functionality that is only available " + "if deal.II was configured with DEAL_II_WITH_SIMPLEX_SUPPORT enabled."); + /** * This function requires support for the FunctionParser library. */ diff --git a/include/deal.II/simplex/grid_generator.h b/include/deal.II/simplex/grid_generator.h index de0cf2d7db..2e68f9e1bd 100644 --- a/include/deal.II/simplex/grid_generator.h +++ b/include/deal.II/simplex/grid_generator.h @@ -27,223 +27,57 @@ DEAL_II_NAMESPACE_OPEN -namespace Simplex +#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT + +namespace GridGenerator { /** - * This namespace provides a collection of functions to generate simplex - * triangulations for some basic geometries. + * 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`. * - * @ingroup simplex + * @ingroup simplex */ - 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 - const 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 - const Point dx((p2[0] - p1[0]) / repetitions[0], - (p2[1] - p1[1]) / repetitions[1], - (p2[2] - p1[2]) / repetitions[2]); - - // 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}; + template + void + subdivided_hyper_rectangle_with_simplices( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize = false); - // 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 + /** + * 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. + * + * @ingroup simplex + */ + template + void + subdivided_hyper_cube_with_simplices(Triangulation &tria, + const unsigned int repetitions, + const double p1 = 0.0, + const double p2 = 1.0, + const bool colorize = false); +} // namespace GridGenerator +#endif DEAL_II_NAMESPACE_CLOSE diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt index 604ce1d1f2..f547691c5a 100644 --- a/source/simplex/CMakeLists.txt +++ b/source/simplex/CMakeLists.txt @@ -17,6 +17,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_unity_include_src fe_lib.cc + grid_generator.cc polynomials.cc quadrature_lib.cc ) @@ -35,6 +36,7 @@ SETUP_SOURCE_LIST("${_unity_include_src}" SET(_inst fe_lib.inst.in + grid_generator.inst.in ) FILE(GLOB _header diff --git a/source/simplex/grid_generator.cc b/source/simplex/grid_generator.cc new file mode 100644 index 0000000000..e17e69e39d --- /dev/null +++ b/source/simplex/grid_generator.cc @@ -0,0 +1,215 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT + +namespace GridGenerator +{ + template + void + subdivided_hyper_rectangle_with_simplices( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize) + { + AssertDimension(dim, spacedim); + + AssertThrow(colorize == false, ExcNotImplemented()); + + std::vector> vertices; + std::vector> cells; + + if (dim == 2) + { + // determine cell sizes + const 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 + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1], + (p2[2] - p1[2]) / repetitions[2]); + + // 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()); + } + + template + void + subdivided_hyper_cube_with_simplices(Triangulation &tria, + const unsigned int repetitions, + const double p1, + const double p2, + const bool colorize) + { + if (dim == 2) + { + subdivided_hyper_rectangle_with_simplices( + tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize); + } + else if (dim == 3) + { + subdivided_hyper_rectangle_with_simplices( + tria, + {{repetitions, repetitions, repetitions}}, + {p1, p1, p1}, + {p2, p2, p2}, + colorize); + } + else + { + AssertThrow(false, ExcNotImplemented()) + } + } +} // namespace GridGenerator + +#endif + +// explicit instantiations +#include "grid_generator.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/simplex/grid_generator.inst.in b/source/simplex/grid_generator.inst.in new file mode 100644 index 0000000000..3c855b64ae --- /dev/null +++ b/source/simplex/grid_generator.inst.in @@ -0,0 +1,37 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +#ifdef DEAL_II_WITH_SIMPLEX_SUPPORT + +for (deal_II_dimension : DIMENSIONS) + { + template void GridGenerator::subdivided_hyper_rectangle_with_simplices( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize); + + template void GridGenerator::subdivided_hyper_cube_with_simplices( + Triangulation & tria, + const unsigned int repetitions, + const double p1, + const double p2, + const bool colorize); + } + +#endif -- 2.39.5