From: David Wells <drwells@email.unc.edu>
Date: Mon, 1 Jul 2024 14:18:24 +0000 (-0400)
Subject: GridTools::Cache: set up the default Mapping in another constructor.
X-Git-Tag: v9.6.0-rc1~142^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e09819225c9ec65da5ce994a95edf227522e4614;p=dealii.git

GridTools::Cache: set up the default Mapping in another constructor.

This preserves the current behavior and makes the Triangulation-only version
work correctly with simplices.
---

diff --git a/include/deal.II/grid/grid_tools_cache.h b/include/deal.II/grid/grid_tools_cache.h
index d7f110e63e..f48710d010 100644
--- a/include/deal.II/grid/grid_tools_cache.h
+++ b/include/deal.II/grid/grid_tools_cache.h
@@ -72,18 +72,19 @@ namespace GridTools
      * If you provide the optional `mapping` argument, then this is used
      * whenever a mapping is required.
      *
-     * @param tria The triangulation for which to store information
-     * @param mapping The mapping to use when computing cached objects
+     * @param tria The triangulation for which to store information.
+     * @param mapping The Mapping to use when computing cached objects.
      */
     Cache(const Triangulation<dim, spacedim> &tria,
-          const Mapping<dim, spacedim>       &mapping =
-            (ReferenceCells::get_hypercube<dim>()
-#ifndef _MSC_VER
-               .template get_default_linear_mapping<dim, spacedim>()
-#else
-               .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
-#endif
-               ));
+          const Mapping<dim, spacedim>       &mapping);
+
+    /**
+     * Constructor. Uses the default linear or multilinear mapping for the
+     * underlying ReferenceCell of @p tria.
+     *
+     * @param tria The triangulation for which to store information.
+     */
+    Cache(const Triangulation<dim, spacedim> &tria);
 
     /**
      * Destructor.
@@ -327,9 +328,14 @@ namespace GridTools
     mutable std::mutex vertices_with_ghost_neighbors_mutex;
 
     /**
-     * Storage for the status of the triangulation signal.
+     * Storage for the status of the triangulation change signal.
+     */
+    boost::signals2::connection tria_change_signal;
+
+    /**
+     * Storage for the status of the triangulation creation signal.
      */
-    boost::signals2::connection tria_signal;
+    boost::signals2::connection tria_create_signal;
   };
 
 
@@ -348,6 +354,7 @@ namespace GridTools
   inline const Mapping<dim, spacedim> &
   Cache<dim, spacedim>::get_mapping() const
   {
+    Assert(mapping, ExcNotInitialized());
     return *mapping;
   }
 } // namespace GridTools
diff --git a/source/grid/grid_tools_cache.cc b/source/grid/grid_tools_cache.cc
index 8f831a3c25..ebb42bc3d8 100644
--- a/source/grid/grid_tools_cache.cc
+++ b/source/grid/grid_tools_cache.cc
@@ -30,17 +30,47 @@ namespace GridTools
     , tria(&tria)
     , mapping(&mapping)
   {
-    tria_signal =
+    tria_change_signal =
       tria.signals.any_change.connect([&]() { mark_for_update(update_all); });
   }
 
+
+
+  template <int dim, int spacedim>
+  Cache<dim, spacedim>::Cache(const Triangulation<dim, spacedim> &tria)
+    : update_flags(update_all)
+    , tria(&tria)
+  {
+    tria_change_signal =
+      tria.signals.any_change.connect([&]() { mark_for_update(update_all); });
+
+    // Allow users to set this class up with an empty Triangulation and no
+    // Mapping argument by deferring Mapping assignment until after the
+    // Triangulation exists.
+    if (tria.get_reference_cells().size() == 0)
+      {
+        tria_create_signal = tria.signals.create.connect([&]() {
+          Assert(tria.get_reference_cells().size() > 0, ExcInternalError());
+          Assert(tria.get_reference_cells().size() == 1, ExcNotImplemented());
+          mapping = &tria.get_reference_cells()[0]
+                       .template get_default_linear_mapping<dim, spacedim>();
+        });
+      }
+    else
+      {
+        Assert(tria.get_reference_cells().size() == 1, ExcNotImplemented());
+        mapping = &tria.get_reference_cells()[0]
+                     .template get_default_linear_mapping<dim, spacedim>();
+      }
+  }
+
   template <int dim, int spacedim>
   Cache<dim, spacedim>::~Cache()
   {
-    // Make sure that the signal that was attached to the triangulation
-    // is removed here.
-    if (tria_signal.connected())
-      tria_signal.disconnect();
+    if (tria_change_signal.connected())
+      tria_change_signal.disconnect();
+    if (tria_create_signal.connected())
+      tria_create_signal.disconnect();
   }
 
 
diff --git a/tests/simplex/compute_point_locations_01.cc b/tests/simplex/compute_point_locations_01.cc
index a171f9308c..33cd75ab5b 100644
--- a/tests/simplex/compute_point_locations_01.cc
+++ b/tests/simplex/compute_point_locations_01.cc
@@ -12,14 +12,6 @@
 
 #include <deal.II/base/quadrature_lib.h>
 
-#include <deal.II/fe/fe_pyramid_p.h>
-#include <deal.II/fe/fe_simplex_p.h>
-#include <deal.II/fe/fe_simplex_p_bubbles.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe_wedge_p.h>
-#include <deal.II/fe/mapping_fe.h>
-
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/grid_tools_cache.h>
@@ -27,7 +19,7 @@
 #include "../tests.h"
 
 // Create a simplex mesh in the unit cube. Check points distribute to each cell
-// via GridTools::compute_point_locations
+// via GridTools::compute_point_locations()
 
 
 template <int dim>
@@ -37,13 +29,9 @@ test_in_unit_cube(const std::vector<Point<dim>> &points)
   Triangulation<dim> tria;
   GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
 
-  MappingFE<dim> mapping(FE_SimplexP<dim>(1));
-
-  const auto tria_cache =
-    std::make_unique<GridTools::Cache<dim>>(tria, mapping);
-
-  const auto point_locations =
-    GridTools::compute_point_locations(*tria_cache, points);
+  GridTools::Cache<dim> cache(tria);
+  const auto            point_locations =
+    GridTools::compute_point_locations(cache, points);
 
   const auto cells   = std::get<0>(point_locations);
   const auto qpoints = std::get<1>(point_locations);
@@ -61,8 +49,8 @@ test_in_unit_cube(const std::vector<Point<dim>> &points)
           deallog << "        physical position   : (" << points[indices[i][j]]
                   << ')' << std::endl;
           deallog << "        FE mapping position : ("
-                  << mapping.transform_unit_to_real_cell(cells[i],
-                                                         qpoints[i][j])
+                  << cache.get_mapping().transform_unit_to_real_cell(
+                       cells[i], qpoints[i][j])
                   << ')' << std::endl;
         }
     }