From d457ad9b34172cd3f0b8e89cea13b242eec3cf94 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 14 Nov 2019 10:22:06 +0100
Subject: [PATCH] Add initialization of MappingQCache by a lambda

---
 include/deal.II/fe/mapping_q_cache.h | 27 +++++++++++++++++++++++++++
 source/fe/mapping_q_cache.cc         | 21 ++++++++++++++++++++-
 2 files changed, 47 insertions(+), 1 deletion(-)

diff --git a/include/deal.II/fe/mapping_q_cache.h b/include/deal.II/fe/mapping_q_cache.h
index 6f1a4225e9..4d9bf652d5 100644
--- a/include/deal.II/fe/mapping_q_cache.h
+++ b/include/deal.II/fe/mapping_q_cache.h
@@ -86,6 +86,33 @@ public:
   initialize(const Triangulation<dim, spacedim> &  triangulation,
              const MappingQGeneric<dim, spacedim> &mapping);
 
+  /**
+   * Initialize the data cache by letting the function given as an argument
+   * provide the mapping support points for all cells (on all levels) of the
+   * given triangulation. The function must return a vector of
+   * `Point<spacedim>` whose length is the same as the size of the polynomial
+   * space, $(p+1)^\text{dim}$, where $p$ is the polynomial degree of the
+   * mapping, and it must be in the order the mapping or FE_Q sort their
+   * points, i.e., all $2^\text{dim}$ vertex points first, then the points on
+   * the lines, quads, and hexes according to the usual hierarchical
+   * numbering. No attempt is made to validate these points internally, except
+   * for the number of given points.
+   *
+   * @note If multiple threads are enabled, this function will run in
+   * parallel, invoking the function passed in several times. Thus, in case
+   * MultithreadInfo::n_threads()>1, the user code must make sure that the
+   * function, typically a lambda, does not write into data shared with other
+   * threads.
+   *
+   * @note The cache is invalidated upon the signal
+   * Triangulation::Signals::any_change of the underlying triangulation.
+   */
+  void
+  initialize(const Triangulation<dim, spacedim> &triangulation,
+             const std::function<std::vector<Point<spacedim>>(
+               const typename Triangulation<dim, spacedim>::cell_iterator &)>
+               &compute_points_on_cell);
+
   /**
    * Return the memory consumption (in bytes) of the cache.
    */
diff --git a/source/fe/mapping_q_cache.cc b/source/fe/mapping_q_cache.cc
index 9c7cae10e6..9f95efb1c6 100644
--- a/source/fe/mapping_q_cache.cc
+++ b/source/fe/mapping_q_cache.cc
@@ -77,7 +77,23 @@ MappingQCache<dim, spacedim>::initialize(
   const MappingQGeneric<dim, spacedim> &mapping)
 {
   AssertDimension(this->get_degree(), mapping.get_degree());
+  initialize(triangulation,
+             [&mapping](const typename Triangulation<dim>::cell_iterator &cell)
+               -> std::vector<Point<spacedim>> {
+               return mapping.compute_mapping_support_points(cell);
+             });
+}
+
 
+
+template <int dim, int spacedim>
+void
+MappingQCache<dim, spacedim>::initialize(
+  const Triangulation<dim, spacedim> &triangulation,
+  const std::function<std::vector<Point<spacedim>>(
+    const typename Triangulation<dim, spacedim>::cell_iterator &)>
+    &compute_points_on_cell)
+{
   clear_signal = triangulation.signals.any_change.connect(
     [&]() -> void { this->support_point_cache.reset(); });
 
@@ -94,7 +110,10 @@ MappingQCache<dim, spacedim>::initialize(
         void *,
         void *) {
       (*support_point_cache)[cell->level()][cell->index()] =
-        mapping.compute_mapping_support_points(cell);
+        compute_points_on_cell(cell);
+      AssertDimension(
+        (*support_point_cache)[cell->level()][cell->index()].size(),
+        Utilities::pow(this->get_degree() + 1, dim));
     },
     /* copier */ std::function<void(void *)>(),
     /* scratch_data */ nullptr,
-- 
2.39.5