From bd838d8f9e9aa5f9c3529040b4b4607e80bf74c5 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 27 Dec 2020 15:11:57 +0100 Subject: [PATCH] Add utility functions to ReferenceCell namespace --- include/deal.II/grid/reference_cell.h | 37 ++++++ source/grid/CMakeLists.txt | 2 + source/grid/reference_cell.cc | 168 ++++++++++++++++++++++++++ source/grid/reference_cell.inst.in | 33 +++++ 4 files changed, 240 insertions(+) create mode 100644 source/grid/reference_cell.cc create mode 100644 source/grid/reference_cell.inst.in diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 2bfa0943b1..595889907b 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -23,6 +23,15 @@ DEAL_II_NAMESPACE_OPEN +// Forward declarations +#ifndef DOXYGEN +template +class Triangulation; +template +class Mapping; +template +class Quadrature; +#endif /** * A namespace for reference cells. @@ -479,6 +488,34 @@ namespace ReferenceCell return {}; } + /* + * Create a (coarse) grid with a single cell of the shape of the provided + * reference cell. + */ + template + void + make_triangulation(const Type & reference_cell, + Triangulation &tria); + + /** + * Return a default linear mapping matching the given reference cell + * (MappingQ1 for hypercube cells and MappingFE else). + */ + template + const Mapping & + get_default_linear_mapping(const Type &reference_cell); + + /** + * Return a Gauss-type quadrature matching the given reference cell(QGauss, + * Simplex::QGauss, Simplex::QGaussPyramid, Simplex::QGaussWedge) and + * @p n_points_1D the number of quadrature points in each direction (QGuass) + * or the indication of what polynomial degree to be integrated exactly. + */ + template + Quadrature + get_gauss_type_quadrature(const Type & reference_cell, + const unsigned n_points_1D); + namespace internal { /** diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index d50cbdec9b..8d0bf1bdae 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -23,6 +23,7 @@ SET(_unity_include_src manifold.cc manifold_lib.cc persistent_tria.cc + reference_cell.cc tria_accessor.cc tria_description.cc tria_faces.cc @@ -64,6 +65,7 @@ SET(_inst intergrid_map.inst.in manifold.inst.in manifold_lib.inst.in + reference_cell.inst.in tria_accessor.inst.in tria_description.inst.in tria.inst.in diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc new file mode 100644 index 0000000000..31d9acacab --- /dev/null +++ b/source/grid/reference_cell.cc @@ -0,0 +1,168 @@ +// --------------------------------------------------------------------- +// +// 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 +#include + +#include +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace ReferenceCell +{ + template + void + make_triangulation(const Type & reference_cell, + Triangulation &tria) + { + AssertDimension(dim, get_dimension(reference_cell)); + + if (reference_cell == get_hypercube(dim)) + { + GridGenerator::hyper_cube(tria, 0, 1); + } + else if (reference_cell == Type::Tri) + { + static const std::vector> vertices = { + {{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}}}; + + std::vector> cells(1); + cells[0].vertices = {0, 1, 2}; + + tria.create_triangulation(vertices, cells, {}); + } + else if (reference_cell == Type::Tet) + { + AssertDimension(spacedim, 3); + + static const std::vector> vertices = { + {{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}}; + + std::vector> cells(1); + cells[0].vertices = {0, 1, 2, 3}; + + tria.create_triangulation(vertices, cells, {}); + } + else if (reference_cell == Type::Pyramid) + { + AssertDimension(spacedim, 3); + + static const std::vector> vertices = { + {{-1.0, -1.0, 0.0}, + {+1.0, -1.0, 0.0}, + {-1.0, +1.0, 0.0}, + {+1.0, +1.0, 0.0}, + {+0.0, +0.0, 1.0}}}; + + std::vector> cells(1); + cells[0].vertices = {0, 1, 2, 3, 4}; + + tria.create_triangulation(vertices, cells, {}); + } + else if (reference_cell == Type::Wedge) + { + AssertDimension(spacedim, 3); + + static const std::vector> vertices = { + {{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 0.0}, + {1.0, 0.0, 1.0}, + {0.0, 1.0, 1.0}, + {0.0, 0.0, 1.0}}}; + + std::vector> cells(1); + cells[0].vertices = {0, 1, 2, 3, 4, 5}; + + tria.create_triangulation(vertices, cells, {}); + } + else + { + Assert(false, ExcNotImplemented()); + } + } + + + template + const Mapping & + get_default_linear_mapping(const Type &reference_cell) + { + AssertDimension(dim, get_dimension(reference_cell)); + + if (reference_cell == get_hypercube(dim)) + { + return StaticMappingQ1::mapping; + } + else if (reference_cell == Type::Tri || reference_cell == Type::Tet) + { + static const MappingFE mapping( + Simplex::FE_P(1)); + return mapping; + } + else if (reference_cell == Type::Pyramid) + { + static const MappingFE mapping( + Simplex::FE_PyramidP(1)); + return mapping; + } + else if (reference_cell == Type::Wedge) + { + static const MappingFE mapping( + Simplex::FE_WedgeP(1)); + return mapping; + } + else + { + Assert(false, ExcNotImplemented()); + } + + return StaticMappingQ1::mapping; // never reached + } + + + + template + Quadrature + get_gauss_type_quadrature(const Type & reference_cell, + const unsigned n_points_1D) + { + AssertDimension(dim, get_dimension(reference_cell)); + + if (reference_cell == get_hypercube(dim)) + return QGauss(n_points_1D); + else if (reference_cell == Type::Tri || reference_cell == Type::Tet) + return Simplex::QGauss(n_points_1D); + else if (reference_cell == Type::Pyramid) + return Simplex::QGaussPyramid(n_points_1D); + else if (reference_cell == Type::Wedge) + return Simplex::QGaussWedge(n_points_1D); + else + Assert(false, ExcNotImplemented()); + + return Quadrature(); // never reached + } + +#include "reference_cell.inst" + +} // namespace ReferenceCell + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/grid/reference_cell.inst.in b/source/grid/reference_cell.inst.in new file mode 100644 index 0000000000..35a6a2ff05 --- /dev/null +++ b/source/grid/reference_cell.inst.in @@ -0,0 +1,33 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template void make_triangulation( + const Type & reference_cell, + Triangulation &tria); + + template const Mapping + &get_default_linear_mapping(const Type &reference_cell); + +#endif + } + +for (deal_II_dimension : DIMENSIONS) + { + template Quadrature get_gauss_type_quadrature( + const Type &reference_cell, const unsigned n_points_1D); + } -- 2.39.5