From 5e6e726dd1d6abeacb5f46213e81c27e718bb795 Mon Sep 17 00:00:00 2001
From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Thu, 2 Jan 2020 15:04:06 -0500
Subject: [PATCH] Template function MatrixFree::get_dof_handler() to work with
 hp::DoFHandler

---
 include/deal.II/matrix_free/matrix_free.h     |  9 ++-
 .../matrix_free/matrix_free.templates.h       | 64 +++++++++++++++----
 source/matrix_free/matrix_free.inst.in        | 13 ++++
 3 files changed, 70 insertions(+), 16 deletions(-)

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index 6539c31885..7e32badbd9 100644
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1520,9 +1520,14 @@ public:
 
   /**
    * Return the DoFHandler with the index as given to the respective
-   * `std::vector` argument in the reinit() function.
+   * `std::vector` argument in the reinit() function. Note that if you want to
+   * call this function with a template parameter different than the default
+   * one, you will need to use the `template` before the function call, i.e.,
+   * you will have something like `matrix_free.template
+   * get_dof_handler<hp::DoFHandler<dim>>()`.
    */
-  const DoFHandler<dim> &
+  template <typename DoFHandlerType = DoFHandler<dim>>
+  const DoFHandlerType &
   get_dof_handler(const unsigned int dof_handler_index = 0) const;
 
   /**
diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index d4de7d1ba1..fd119ddb61 100644
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -53,6 +53,50 @@ DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 DEAL_II_NAMESPACE_OPEN
 
 
+namespace internal
+{
+  template <
+    typename DoFHandlerType,
+    typename std::enable_if<
+      std::is_same<DoFHandlerType,
+                   ::dealii::DoFHandler<DoFHandlerType::dimension>>::value,
+      int>::type = 0>
+  const DoFHandlerType &
+  get_dof_handler(
+    const std::vector<
+      SmartPointer<const ::dealii::DoFHandler<DoFHandlerType::dimension>>>
+      &dof_handler,
+    const std::vector<
+      SmartPointer<const ::dealii::hp::DoFHandler<DoFHandlerType::dimension>>>
+      &,
+    const unsigned int n_dof_handlers,
+    const unsigned int dof_handler_index)
+  {
+    AssertDimension(dof_handler.size(), n_dof_handlers);
+    return *dof_handler[dof_handler_index];
+  }
+
+  template <
+    typename DoFHandlerType,
+    typename std::enable_if<
+      std::is_same<DoFHandlerType,
+                   ::dealii::hp::DoFHandler<DoFHandlerType::dimension>>::value,
+      int>::type = 0>
+  const DoFHandlerType &
+  get_dof_handler(
+    const std::vector<
+      SmartPointer<const ::dealii::DoFHandler<DoFHandlerType::dimension>>> &,
+    const std::vector<
+      SmartPointer<const ::dealii::hp::DoFHandler<DoFHandlerType::dimension>>>
+      &                hp_dof_handler,
+    const unsigned int n_dof_handlers,
+    const unsigned int dof_handler_index)
+  {
+    AssertDimension(hp_dof_handler.size(), n_dof_handlers);
+    return *hp_dof_handler[dof_handler_index];
+  }
+} // namespace internal
+
 // --------------------- MatrixFree -----------------------------------
 
 template <int dim, typename Number, typename VectorizedArrayType>
@@ -132,24 +176,16 @@ MatrixFree<dim, Number, VectorizedArrayType>::renumber_dofs(
 
 
 template <int dim, typename Number, typename VectorizedArrayType>
-const DoFHandler<dim> &
+template <typename DoFHandlerType>
+const DoFHandlerType &
 MatrixFree<dim, Number, VectorizedArrayType>::get_dof_handler(
   const unsigned int dof_handler_index) const
 {
   AssertIndexRange(dof_handler_index, n_components());
-  if (dof_handlers.active_dof_handler == DoFHandlers::usual)
-    {
-      AssertDimension(dof_handlers.dof_handler.size(),
-                      dof_handlers.n_dof_handlers);
-      return *dof_handlers.dof_handler[dof_handler_index];
-    }
-  else
-    {
-      AssertThrow(false, ExcNotImplemented());
-      // put pseudo return argument to avoid compiler error, but trigger a
-      // segfault in case this is only run in optimized mode
-      return *dof_handlers.dof_handler[numbers::invalid_unsigned_int];
-    }
+  return internal::get_dof_handler<DoFHandlerType>(dof_handlers.dof_handler,
+                                                   dof_handlers.hp_dof_handler,
+                                                   dof_handlers.n_dof_handlers,
+                                                   dof_handler_index);
 }
 
 
diff --git a/source/matrix_free/matrix_free.inst.in b/source/matrix_free/matrix_free.inst.in
index 9fe8f0f8fb..3485264e83 100644
--- a/source/matrix_free/matrix_free.inst.in
+++ b/source/matrix_free/matrix_free.inst.in
@@ -66,6 +66,19 @@ for (deal_II_dimension : DIMENSIONS;
         const std::vector<hp::QCollection<1>> &,
         const AdditionalData &);
 
+    template const DoFHandler<deal_II_dimension> &
+    MatrixFree<deal_II_dimension,
+               deal_II_scalar_vectorized::value_type,
+               deal_II_scalar_vectorized>::
+      get_dof_handler<DoFHandler<deal_II_dimension>>(const unsigned int) const;
+
+    template const hp::DoFHandler<deal_II_dimension>
+      &MatrixFree<deal_II_dimension,
+                  deal_II_scalar_vectorized::value_type,
+                  deal_II_scalar_vectorized>::
+        get_dof_handler<hp::DoFHandler<deal_II_dimension>>(const unsigned int)
+          const;
+
     template void internal::MatrixFreeFunctions::
       ShapeInfo<deal_II_scalar_vectorized>::reinit<deal_II_dimension>(
         const Quadrature<1> &,
-- 
2.39.5