From eb09f25d08a1c9f92efe21ac6a6651cf1a2b1d2a Mon Sep 17 00:00:00 2001
From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Thu, 6 Jan 2022 20:51:36 +0100
Subject: [PATCH] Add method that returns the default linear mapping from a
 hp::FECollection

---
 include/deal.II/hp/fe_collection.h      | 26 ++++++++++++++++++
 include/deal.II/hp/mapping_collection.h | 34 +++++++++++++++++++-----
 source/hp/fe_collection.cc              | 35 +++++++++++++++++++++++++
 3 files changed, 89 insertions(+), 6 deletions(-)

diff --git a/include/deal.II/hp/fe_collection.h b/include/deal.II/hp/fe_collection.h
index abc10fb004..6fc19edd91 100644
--- a/include/deal.II/hp/fe_collection.h
+++ b/include/deal.II/hp/fe_collection.h
@@ -28,6 +28,14 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+// Forward declarations
+namespace hp
+{
+  template <int dim, int spacedim>
+  class MappingCollection;
+}
+
+
 namespace hp
 {
   /**
@@ -268,6 +276,17 @@ namespace hp
     unsigned int
     max_dofs_per_cell() const;
 
+    /**
+     * Return a mapping collection that consists of the default linear mappings
+     * matching the reference cells for each hp index. More details may be found
+     * in the documentation for ReferenceCell::get_default_linear_mapping().
+     *
+     * @note This FECollection object must remain in scope for as long as the
+     * reference cell default linear mapping is in use.
+     */
+    const MappingCollection<dim, spacedim> &
+    get_reference_cell_default_linear_mapping() const;
+
     /**
      * @}
      */
@@ -829,6 +848,13 @@ namespace hp
      */
 
   private:
+    /**
+     * A linear mapping collection for all reference cell types of each index
+     * of this object.
+     */
+    std::shared_ptr<MappingCollection<dim, spacedim>>
+      reference_cell_default_linear_mapping;
+
     /**
      * %Function returning the index of the finite element following the given
      * one in hierarchy.
diff --git a/include/deal.II/hp/mapping_collection.h b/include/deal.II/hp/mapping_collection.h
index 5b309eb401..34006e0104 100644
--- a/include/deal.II/hp/mapping_collection.h
+++ b/include/deal.II/hp/mapping_collection.h
@@ -70,12 +70,6 @@ namespace hp
      */
     explicit MappingCollection(const Mapping<dim, spacedim> &mapping);
 
-    /**
-     * Copy constructor.
-     */
-    MappingCollection(
-      const MappingCollection<dim, spacedim> &mapping_collection);
-
     /**
      * Constructor. This constructor creates a MappingCollection from one or
      * more mapping objects passed to the constructor. For this
@@ -85,6 +79,34 @@ namespace hp
     template <class... MappingTypes>
     explicit MappingCollection(const MappingTypes &...mappings);
 
+    /**
+     * Copy constructor.
+     */
+    MappingCollection(
+      const MappingCollection<dim, spacedim> &mapping_collection);
+
+    /**
+     * Move constructor.
+     *
+     * @note The implementation of standard datatypes may change with different
+     * libraries, so their move members may or may not be flagged non-throwing.
+     * We need to explicitly set the noexcept specifier according to its
+     * member variables to still get the performance benefits (and to satisfy
+     * clang-tidy).
+     */
+    MappingCollection(MappingCollection<dim, spacedim> &&) noexcept(
+      std::is_nothrow_move_constructible<
+        std::vector<std::shared_ptr<const Mapping<dim, spacedim>>>>::value
+        &&std::is_nothrow_move_constructible<std::function<
+          unsigned int(const typename hp::MappingCollection<dim, spacedim> &,
+                       const unsigned int)>>::value) = default;
+
+    /**
+     * Move assignment operator.
+     */
+    MappingCollection<dim, spacedim> &
+    operator=(MappingCollection<dim, spacedim> &&) = default; // NOLINT
+
     /**
      * Add a new mapping to the MappingCollection. Generally, you will
      * want to use the same order for mappings as for the elements of
diff --git a/source/hp/fe_collection.cc b/source/hp/fe_collection.cc
index f8022758fc..b171fa3afd 100644
--- a/source/hp/fe_collection.cc
+++ b/source/hp/fe_collection.cc
@@ -17,6 +17,7 @@
 #include <deal.II/base/memory_consumption.h>
 
 #include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/mapping_collection.h>
 
 #include <set>
 
@@ -76,6 +77,40 @@ namespace hp
 
 
 
+  template <int dim, int spacedim>
+  const MappingCollection<dim, spacedim> &
+  FECollection<dim, spacedim>::get_reference_cell_default_linear_mapping() const
+  {
+    Assert(this->size() > 0, ExcNoFiniteElements());
+
+    // Since we can only add elements to an FECollection, we are safe comparing
+    // the sizes of this object and the MappingCollection. One circumstance that
+    // might lead to their sizes diverging is this:
+    // - An FECollection is created and then this function is called. The shared
+    //   map is now initialized.
+    // - A second FECollection is made as a copy of this one. The shared map is
+    //   not recreated.
+    // - The second FECollection is then resized by adding a new FE. The shared
+    //   map is thus invalid for the second instance.
+    if (!reference_cell_default_linear_mapping ||
+        reference_cell_default_linear_mapping->size() != this->size())
+      {
+        auto &this_nc = const_cast<FECollection<dim, spacedim> &>(*this);
+
+        this_nc.reference_cell_default_linear_mapping =
+          std::make_shared<MappingCollection<dim, spacedim>>();
+
+        for (const auto &fe : *this)
+          this_nc.reference_cell_default_linear_mapping->push_back(
+            fe.reference_cell()
+              .template get_default_linear_mapping<dim, spacedim>());
+      }
+
+    return *reference_cell_default_linear_mapping;
+  }
+
+
+
   template <int dim, int spacedim>
   std::set<unsigned int>
   FECollection<dim, spacedim>::find_common_fes(
-- 
2.39.5