From e9f77d0444a0564a56596ac75a05c2dc5d445377 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 3 Nov 2020 23:11:34 +0100
Subject: [PATCH] Store reference_types in Triangulation

---
 include/deal.II/distributed/tria_base.h |  6 +++
 include/deal.II/grid/tria.h             | 26 +++++++++++
 source/distributed/tria_base.cc         | 26 +++++++++++
 source/grid/tria.cc                     | 62 ++++++++++++++++++++-----
 4 files changed, 108 insertions(+), 12 deletions(-)

diff --git a/include/deal.II/distributed/tria_base.h b/include/deal.II/distributed/tria_base.h
index 167f9926ad..566f9fc53c 100644
--- a/include/deal.II/distributed/tria_base.h
+++ b/include/deal.II/distributed/tria_base.h
@@ -318,6 +318,12 @@ namespace parallel
     virtual void
     update_number_cache();
 
+    /**
+     * @copydoc dealii::Triangulation::update_reference_cell_types()
+     */
+    void
+    update_reference_cell_types() override;
+
     /**
      * Reset global active cell indices and global level cell indices.
      */
diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h
index 96abef56a3..5a8efeeb31 100644
--- a/include/deal.II/grid/tria.h
+++ b/include/deal.II/grid/tria.h
@@ -3359,6 +3359,20 @@ public:
     std::pair<std::pair<cell_iterator, unsigned int>, std::bitset<3>>> &
   get_periodic_face_map() const;
 
+  /**
+   * Return vector filled with the used reference-cell types of this
+   * triangulation.
+   */
+  const std::vector<ReferenceCell::Type> &
+  get_reference_cell_types() const;
+
+  /**
+   * Indicate if the triangulation only consists of hypercube-like cells, i.e.,
+   * lines, quadrilaterals, or hexahedrons.
+   */
+  bool
+  all_reference_cell_types_are_hyper_cube() const;
+
 #ifdef DOXYGEN
   /**
    * Write and read the data of this object from a stream for the purpose
@@ -3467,6 +3481,12 @@ protected:
    */
   MeshSmoothing smooth_grid;
 
+  /**
+   * Vector caching all reference-cell types of the given triangulation
+   * (also in the distributed case).
+   */
+  std::vector<ReferenceCell::Type> reference_cell_types;
+
   /**
    * Write a bool vector to the given stream, writing a pre- and a postfix
    * magic number. The vector is written in an almost binary format, i.e. the
@@ -3503,6 +3523,12 @@ protected:
   void
   update_periodic_face_map();
 
+  /**
+   * Update the internal reference_cell_types vector.
+   */
+  virtual void
+  update_reference_cell_types();
+
 
 private:
   /**
diff --git a/source/distributed/tria_base.cc b/source/distributed/tria_base.cc
index 2ba8d4f359..97d475ddaf 100644
--- a/source/distributed/tria_base.cc
+++ b/source/distributed/tria_base.cc
@@ -287,6 +287,32 @@ namespace parallel
 
 #endif
 
+  template <int dim, int spacedim>
+  void
+  TriangulationBase<dim, spacedim>::update_reference_cell_types()
+  {
+    // run algorithm for locally-owned cells
+    dealii::Triangulation<dim, spacedim>::update_reference_cell_types();
+
+    // translate ReferenceCell::Type to unsigned int (needed by
+    // Utilities::MPI::compute_set_union)
+    std::vector<unsigned int> reference_cell_types_ui;
+
+    for (const auto &i : this->reference_cell_types)
+      reference_cell_types_ui.push_back(static_cast<unsigned int>(i));
+
+    // create union
+    reference_cell_types_ui =
+      Utilities::MPI::compute_set_union(reference_cell_types_ui,
+                                        this->mpi_communicator);
+
+    // transform back and store result
+    this->reference_cell_types.clear();
+    for (const auto &i : reference_cell_types_ui)
+      this->reference_cell_types.push_back(static_cast<ReferenceCell::Type>(i));
+  }
+
+
   template <int dim, int spacedim>
   types::subdomain_id
   TriangulationBase<dim, spacedim>::locally_owned_subdomain() const
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index a8e053e185..5912648524 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -9402,6 +9402,7 @@ Triangulation<dim, spacedim>::Triangulation(
   Triangulation<dim, spacedim> &&tria) noexcept
   : Subscriptor(std::move(tria))
   , smooth_grid(tria.smooth_grid)
+  , reference_cell_types(std::move(tria.reference_cell_types))
   , periodic_face_pairs_level_0(std::move(tria.periodic_face_pairs_level_0))
   , periodic_face_map(std::move(tria.periodic_face_map))
   , levels(std::move(tria.levels))
@@ -9427,6 +9428,7 @@ operator=(Triangulation<dim, spacedim> &&tria) noexcept
   Subscriptor::operator=(std::move(tria));
 
   smooth_grid                  = tria.smooth_grid;
+  reference_cell_types         = std::move(tria.reference_cell_types);
   periodic_face_pairs_level_0  = std::move(tria.periodic_face_pairs_level_0);
   periodic_face_map            = std::move(tria.periodic_face_map);
   levels                       = std::move(tria.levels);
@@ -9485,6 +9487,7 @@ Triangulation<dim, spacedim>::clear()
   clear_despite_subscriptions();
   periodic_face_pairs_level_0.clear();
   periodic_face_map.clear();
+  reference_cell_types.clear();
 }
 
 
@@ -9714,6 +9717,7 @@ Triangulation<dim, spacedim>::copy_triangulation(
   vertices_used          = other_tria.vertices_used;
   anisotropic_refinement = other_tria.anisotropic_refinement;
   smooth_grid            = other_tria.smooth_grid;
+  reference_cell_types   = other_tria.reference_cell_types;
 
   if (dim > 1)
     faces = std::make_unique<internal::TriangulationImplementation::TriaFaces>(
@@ -9802,20 +9806,16 @@ Triangulation<dim, spacedim>::create_triangulation(
   // because sometimes other objects are already attached to it:
   try
     {
-#ifndef DEAL_II_WITH_SIMPLEX_SUPPORT
-      AssertThrow(
-        std::any_of(cells.begin(),
-                    cells.end(),
-                    [](const auto &cell) {
-                      return cell.vertices.size() !=
-                             GeometryInfo<dim>::vertices_per_cell;
-                    }) == false,
-        ExcMessage(
-          "A cell with invalid number of vertices has been provided."));
-#endif
-
       internal::TriangulationImplementation::Implementation::
         create_triangulation(v, cells, subcelldata, *this);
+
+      this->update_reference_cell_types();
+
+#ifndef DEAL_II_WITH_SIMPLEX_SUPPORT
+      Assert(this->all_reference_cell_types_are_hyper_cube(),
+             ExcMessage(
+               "A cell with invalid number of vertices has been provided."));
+#endif
     }
   catch (...)
     {
@@ -12772,6 +12772,43 @@ Triangulation<dim, spacedim>::update_periodic_face_map()
 
 
 
+template <int dim, int spacedim>
+void
+Triangulation<dim, spacedim>::update_reference_cell_types()
+{
+  std::set<ReferenceCell::Type> reference_cell_types_set;
+  for (auto cell : active_cell_iterators())
+    if (cell->is_locally_owned())
+      reference_cell_types_set.insert(cell->reference_cell_type());
+
+  std::vector<ReferenceCell::Type> reference_cell_types(
+    reference_cell_types_set.begin(), reference_cell_types_set.end());
+
+  this->reference_cell_types = reference_cell_types;
+}
+
+
+
+template <int dim, int spacedim>
+const std::vector<ReferenceCell::Type> &
+Triangulation<dim, spacedim>::get_reference_cell_types() const
+{
+  return this->reference_cell_types;
+}
+
+
+
+template <int dim, int spacedim>
+bool
+Triangulation<dim, spacedim>::all_reference_cell_types_are_hyper_cube() const
+{
+  return (this->reference_cell_types.size() == 0) ||
+         (this->reference_cell_types.size() == 1 &&
+          this->reference_cell_types[0] == ReferenceCell::get_hypercube(dim));
+}
+
+
+
 template <int dim, int spacedim>
 void
 Triangulation<dim, spacedim>::clear_despite_subscriptions()
@@ -12788,6 +12825,7 @@ Triangulation<dim, spacedim>::clear_despite_subscriptions()
 }
 
 
+
 template <int dim, int spacedim>
 typename Triangulation<dim, spacedim>::DistortedCellList
 Triangulation<dim, spacedim>::execute_refinement()
-- 
2.39.5