From 20110d0fac152127691b45bf42b7d99e5f5fff09 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Wed, 13 Jan 2021 16:38:38 -0700
Subject: [PATCH] Introduce ReferenceCell::get_default_linear_mapping() with
 triangulation argument.

---
 include/deal.II/grid/reference_cell.h | 22 ++++++++++++++++++++--
 source/grid/reference_cell.cc         | 20 ++++++++++++++++++++
 source/grid/reference_cell.inst.in    |  4 ++++
 3 files changed, 44 insertions(+), 2 deletions(-)

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 582f3499b9..bf61c466c0 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -340,13 +340,31 @@ namespace ReferenceCell
                      Triangulation<dim, spacedim> &tria);
 
   /**
-   * Return a default linear mapping matching the given reference cell
-   * (MappingQ1 for hypercube cells and MappingFE else).
+   * Return a default linear mapping matching the given reference cell.
+   * If this reference cell is a hypercube, then the returned mapping
+   * is a MappingQ1; otherwise, it is an object of type MappingFE
+   * initialized with Simplex::FE_P (if the reference cell is a triangle and
+   * tetrahedron), with Simplex::FE_PyramidP (if the reference cell is a
+   * pyramid), or with Simplex::FE_WedgeP (if the reference cell is a wedge). In
+   * other words, the term "linear" in the name of the function has to be
+   * understood as $d$-linear (i.e., bilinear or trilinear) for some of the
+   * coordinate directions.
    */
   template <int dim, int spacedim>
   const Mapping<dim, spacedim> &
   get_default_linear_mapping(const Type &reference_cell);
 
+  /**
+   * Return a default linear mapping that works for the given triangulation.
+   * Internally, this function calls the function above for the reference
+   * cell used by the given triangulation, assuming that the triangulation
+   * uses only a single cell type. If the triangulation uses mixed cell
+   * types, then this function will trigger an exception.
+   */
+  template <int dim, int spacedim>
+  const Mapping<dim, spacedim> &
+  get_default_linear_mapping(const Triangulation<dim, spacedim> &triangulation);
+
   /**
    * Return a Gauss-type quadrature matching the given reference cell(QGauss,
    * Simplex::QGauss, Simplex::QGaussPyramid, Simplex::QGaussWedge) and
diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc
index 31d9acacab..816da89dad 100644
--- a/source/grid/reference_cell.cc
+++ b/source/grid/reference_cell.cc
@@ -140,6 +140,26 @@ namespace ReferenceCell
 
 
 
+  template <int dim, int spacedim>
+  const Mapping<dim, spacedim> &
+  get_default_linear_mapping(const Triangulation<dim, spacedim> &triangulation)
+  {
+    const auto &reference_cell_types = triangulation.get_reference_cell_types();
+    Assert(reference_cell_types.size() == 1,
+           ExcMessage(
+             "This function can only work for triangulations that "
+             "use only a single cell type -- for example, only triangles "
+             "or only quadrilaterals. For mixed meshes, there is no "
+             "single linear mapping object that can be used for all "
+             "cells of the triangulation. The triangulation you are "
+             "passing to this function uses multiple cell types."));
+
+    return get_default_linear_mapping<dim, spacedim>(
+      reference_cell_types.front());
+  }
+
+
+
   template <int dim>
   Quadrature<dim>
   get_gauss_type_quadrature(const Type &   reference_cell,
diff --git a/source/grid/reference_cell.inst.in b/source/grid/reference_cell.inst.in
index 35a6a2ff05..198079f766 100644
--- a/source/grid/reference_cell.inst.in
+++ b/source/grid/reference_cell.inst.in
@@ -23,6 +23,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
     template const Mapping<deal_II_dimension, deal_II_space_dimension>
       &get_default_linear_mapping(const Type &reference_cell);
 
+    template const Mapping<deal_II_dimension, deal_II_space_dimension>
+      &get_default_linear_mapping(
+        const Triangulation<deal_II_dimension, deal_II_space_dimension>
+          &triangulation);
 #endif
   }
 
-- 
2.39.5