From c1d316ad1c75a76e1850ef05d6c44684a481bfe5 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@uni-a.de>
Date: Mon, 30 May 2022 22:34:47 +0200
Subject: [PATCH] Make VectorTools::project compile with complex numbers

---
 .../numerics/vector_tools_project.templates.h | 109 ++++++++++--------
 1 file changed, 60 insertions(+), 49 deletions(-)

diff --git a/include/deal.II/numerics/vector_tools_project.templates.h b/include/deal.II/numerics/vector_tools_project.templates.h
index 45ef53965d..719a20fc12 100644
--- a/include/deal.II/numerics/vector_tools_project.templates.h
+++ b/include/deal.II/numerics/vector_tools_project.templates.h
@@ -445,6 +445,8 @@ namespace VectorTools
       vec_result.compress(VectorOperation::insert);
     }
 
+
+
     /**
      * Return whether the boundary values try to constrain a degree of freedom
      * that is already constrained to something else
@@ -769,19 +771,23 @@ namespace VectorTools
     }
 
     /**
-     * Specialization of project() for the case dim==spacedim.
-     * Check if we can use the MatrixFree implementation or need
-     * to use the matrix based one.
+     * Specialization of project() for the case dim==spacedim and with correct
+     * number types for MatrixFree support.  Check if we actually can use the
+     * MatrixFree implementation or need to use the matrix based one
+     * nonetheless based on the number of components.
      */
-    template <typename VectorType, int dim>
-    void
+    template <typename VectorType, int dim, int spacedim>
+    std::enable_if_t<
+      !numbers::NumberTraits<typename VectorType::value_type>::is_complex &&
+        (dim == spacedim),
+      void>
     project(
-      const Mapping<dim> &                                      mapping,
-      const DoFHandler<dim> &                                   dof,
-      const AffineConstraints<typename VectorType::value_type> &constraints,
-      const Quadrature<dim> &                                   quadrature,
-      const Function<dim, typename VectorType::value_type> &    function,
-      VectorType &                                              vec_result,
+      const Mapping<dim, spacedim> &                             mapping,
+      const DoFHandler<dim, spacedim> &                          dof,
+      const AffineConstraints<typename VectorType::value_type> & constraints,
+      const Quadrature<dim> &                                    quadrature,
+      const Function<spacedim, typename VectorType::value_type> &function,
+      VectorType &                                               vec_result,
       const bool                 enforce_zero_boundary,
       const Quadrature<dim - 1> &q_boundary,
       const bool                 project_to_boundary_first)
@@ -826,6 +832,40 @@ namespace VectorTools
                      project_to_boundary_first);
         }
     }
+
+    /**
+     * Specialization of project() for complex numbers or `dim < spacedim`,
+     * for which we are sure that we cannot use the MatrixFree implementation.
+     */
+    template <typename VectorType, int dim, int spacedim>
+    std::enable_if_t<
+      numbers::NumberTraits<typename VectorType::value_type>::is_complex ||
+        (dim < spacedim),
+      void>
+    project(
+      const Mapping<dim, spacedim> &                             mapping,
+      const DoFHandler<dim, spacedim> &                          dof,
+      const AffineConstraints<typename VectorType::value_type> & constraints,
+      const Quadrature<dim> &                                    quadrature,
+      const Function<spacedim, typename VectorType::value_type> &function,
+      VectorType &                                               vec_result,
+      const bool                 enforce_zero_boundary,
+      const Quadrature<dim - 1> &q_boundary,
+      const bool                 project_to_boundary_first)
+    {
+      Assert((dynamic_cast<const parallel::TriangulationBase<dim> *>(
+                &(dof.get_triangulation())) == nullptr),
+             ExcNotImplemented());
+      do_project(mapping,
+                 dof,
+                 constraints,
+                 quadrature,
+                 function,
+                 vec_result,
+                 enforce_zero_boundary,
+                 q_boundary,
+                 project_to_boundary_first);
+    }
   } // namespace internal
 
   template <int dim, typename VectorType, int spacedim>
@@ -900,44 +940,15 @@ namespace VectorTools
           const Quadrature<dim - 1> &q_boundary,
           const bool                 project_to_boundary_first)
   {
-    if (dim == spacedim)
-      {
-        const Mapping<dim> *const mapping_ptr =
-          dynamic_cast<const Mapping<dim> *>(&mapping);
-        const DoFHandler<dim> *const dof_ptr =
-          dynamic_cast<const DoFHandler<dim> *>(&dof);
-        const Function<dim,
-                       typename VectorType::value_type> *const function_ptr =
-          dynamic_cast<const Function<dim, typename VectorType::value_type> *>(
-            &function);
-        Assert(mapping_ptr != nullptr, ExcInternalError());
-        Assert(dof_ptr != nullptr, ExcInternalError());
-        internal::project<VectorType, dim>(*mapping_ptr,
-                                           *dof_ptr,
-                                           constraints,
-                                           quadrature,
-                                           *function_ptr,
-                                           vec_result,
-                                           enforce_zero_boundary,
-                                           q_boundary,
-                                           project_to_boundary_first);
-      }
-    else
-      {
-        Assert(
-          (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
-             &(dof.get_triangulation())) == nullptr),
-          ExcNotImplemented());
-        internal::do_project(mapping,
-                             dof,
-                             constraints,
-                             quadrature,
-                             function,
-                             vec_result,
-                             enforce_zero_boundary,
-                             q_boundary,
-                             project_to_boundary_first);
-      }
+    internal::project(mapping,
+                      dof,
+                      constraints,
+                      quadrature,
+                      function,
+                      vec_result,
+                      enforce_zero_boundary,
+                      q_boundary,
+                      project_to_boundary_first);
   }
 
 
-- 
2.39.5