Add utility functions to ReferenceCell namespace

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 2bfa0943b121fe729a3fccaa587c60a57efb49c8..595889907b8e5c7eb044cc3fe0002822abd6266d 100644 (file)
--- 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 <int dim, int spacedim>
+class Triangulation;
+template <int dim, int spacedim>
+class Mapping;
+template <int dim>
+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 <int dim, int spacedim>
+  void
+  make_triangulation(const Type &                  reference_cell,
+                     Triangulation<dim, spacedim> &tria);
+
+  /**
+   * Return a default linear mapping matching the given reference cell
+   * (MappingQ1 for hypercube cells and MappingFE else).
+   */
+  template <int dim, int spacedim>
+  const Mapping<dim, spacedim> &
+  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 <int dim>
+  Quadrature<dim>
+  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 d50cbdec9b389d2088f0e38118447e952002c570..8d0bf1bdaed8e73f1653c7026537d70505cfc962 100644 (file)
--- 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 (file)
index 0000000..31d9aca
--- /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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/reference_cell.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace ReferenceCell
+{
+  template <int dim, int spacedim>
+  void
+  make_triangulation(const Type &                  reference_cell,
+                     Triangulation<dim, spacedim> &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<Point<spacedim>> vertices = {
+          {{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}}};
+
+        std::vector<CellData<dim>> 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<Point<spacedim>> 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<CellData<dim>> 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<Point<spacedim>> 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<CellData<dim>> 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<Point<spacedim>> 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<CellData<dim>> cells(1);
+        cells[0].vertices = {0, 1, 2, 3, 4, 5};
+
+        tria.create_triangulation(vertices, cells, {});
+      }
+    else
+      {
+        Assert(false, ExcNotImplemented());
+      }
+  }
+
+
+  template <int dim, int spacedim>
+  const Mapping<dim, spacedim> &
+  get_default_linear_mapping(const Type &reference_cell)
+  {
+    AssertDimension(dim, get_dimension(reference_cell));
+
+    if (reference_cell == get_hypercube(dim))
+      {
+        return StaticMappingQ1<dim, spacedim>::mapping;
+      }
+    else if (reference_cell == Type::Tri || reference_cell == Type::Tet)
+      {
+        static const MappingFE<dim, spacedim> mapping(
+          Simplex::FE_P<dim, spacedim>(1));
+        return mapping;
+      }
+    else if (reference_cell == Type::Pyramid)
+      {
+        static const MappingFE<dim, spacedim> mapping(
+          Simplex::FE_PyramidP<dim, spacedim>(1));
+        return mapping;
+      }
+    else if (reference_cell == Type::Wedge)
+      {
+        static const MappingFE<dim, spacedim> mapping(
+          Simplex::FE_WedgeP<dim, spacedim>(1));
+        return mapping;
+      }
+    else
+      {
+        Assert(false, ExcNotImplemented());
+      }
+
+    return StaticMappingQ1<dim, spacedim>::mapping; // never reached
+  }
+
+
+
+  template <int dim>
+  Quadrature<dim>
+  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<dim>(n_points_1D);
+    else if (reference_cell == Type::Tri || reference_cell == Type::Tet)
+      return Simplex::QGauss<dim>(n_points_1D);
+    else if (reference_cell == Type::Pyramid)
+      return Simplex::QGaussPyramid<dim>(n_points_1D);
+    else if (reference_cell == Type::Wedge)
+      return Simplex::QGaussWedge<dim>(n_points_1D);
+    else
+      Assert(false, ExcNotImplemented());
+
+    return Quadrature<dim>(); // 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 (file)
index 0000000..35a6a2f
--- /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<deal_II_dimension, deal_II_space_dimension> &tria);
+
+    template const Mapping<deal_II_dimension, deal_II_space_dimension>
+      &get_default_linear_mapping(const Type &reference_cell);
+
+#endif
+  }
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+    template Quadrature<deal_II_dimension> get_gauss_type_quadrature(
+      const Type &reference_cell, const unsigned n_points_1D);
+  }