From 559b632571c28c7b26211f3f70ee023a83647bcc Mon Sep 17 00:00:00 2001
From: Matthias Maier <tamiko@43-1.org>
Date: Sun, 27 May 2018 19:44:37 -0500
Subject: [PATCH] MatrixFree: Allow for non-matching number types in constraint
 handling

In order to facilitate an AffineConstraints<float> in the matrix free
framework:

 - templatify internal ConstraintValues::insert_entries to convert from
   compatible number2 types to Number types (aka exclusively 'double').
---
 .../deal.II/matrix_free/dof_info.templates.h  | 18 +++++++++------
 source/matrix_free/matrix_free.inst.in        | 22 +++++++++++++++++--
 2 files changed, 31 insertions(+), 9 deletions(-)

diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h
index 93d9de3a23..2acb55d565 100644
--- a/include/deal.II/matrix_free/dof_info.templates.h
+++ b/include/deal.II/matrix_free/dof_info.templates.h
@@ -61,9 +61,11 @@ namespace internal
        * new_indices, and returns the storage position the double array for
        * access later on.
        */
+      template <typename number2>
       unsigned short
       insert_entries(
-        const std::vector<std::pair<types::global_dof_index, double>> &entries);
+        const std::vector<std::pair<types::global_dof_index, number2>>
+          &entries);
 
       std::vector<std::pair<types::global_dof_index, double>>
                                            constraint_entries;
@@ -83,15 +85,18 @@ namespace internal
     {}
 
     template <typename Number>
+    template <typename number2>
     unsigned short
     ConstraintValues<Number>::insert_entries(
-      const std::vector<std::pair<types::global_dof_index, double>> &entries)
+      const std::vector<std::pair<types::global_dof_index, number2>> &entries)
     {
       next_constraint.first.resize(entries.size());
       if (entries.size() > 0)
         {
           constraint_indices.resize(entries.size());
-          constraint_entries = entries;
+          // Use assign so that values for nonmatching Number / number2 are
+          // converted:
+          constraint_entries.assign(entries.begin(), entries.end());
           std::sort(constraint_entries.begin(),
                     constraint_entries.end(),
                     ConstraintComparator());
@@ -208,8 +213,8 @@ namespace internal
             {
               types::global_dof_index current_dof =
                 local_indices[lexicographic_inv[i]];
-              const std::vector<std::pair<types::global_dof_index, double>>
-                *entries_ptr = constraints.get_constraint_entries(current_dof);
+              const auto *entries_ptr =
+                constraints.get_constraint_entries(current_dof);
 
               // dof is constrained
               if (entries_ptr != nullptr)
@@ -226,8 +231,7 @@ namespace internal
                   // check whether this dof is identity constrained to another
                   // dof. then we can simply insert that dof and there is no
                   // need to actually resolve the constraint entries
-                  const std::vector<std::pair<types::global_dof_index, double>>
-                    &                           entries   = *entries_ptr;
+                  const auto &                  entries   = *entries_ptr;
                   const types::global_dof_index n_entries = entries.size();
                   if (n_entries == 1 &&
                       std::abs(entries[0].second - 1.) <
diff --git a/source/matrix_free/matrix_free.inst.in b/source/matrix_free/matrix_free.inst.in
index 694b7904df..b7cdcee976 100644
--- a/source/matrix_free/matrix_free.inst.in
+++ b/source/matrix_free/matrix_free.inst.in
@@ -34,7 +34,8 @@ for (deal_II_dimension : DIMENSIONS)
     template void MatrixFree<deal_II_dimension, float>::
       print_memory_consumption<ConditionalOStream>(ConditionalOStream &) const;
 
-    template void MatrixFree<deal_II_dimension, double>::internal_reinit<double>(
+    template void
+    MatrixFree<deal_II_dimension, double>::internal_reinit<double>(
       const Mapping<deal_II_dimension> &,
       const std::vector<const DoFHandler<deal_II_dimension> *> &,
       const std::vector<const AffineConstraints<double> *> &,
@@ -42,7 +43,8 @@ for (deal_II_dimension : DIMENSIONS)
       const std::vector<hp::QCollection<1>> &,
       const AdditionalData &);
 
-    template void MatrixFree<deal_II_dimension, double>::internal_reinit<double>(
+    template void
+    MatrixFree<deal_II_dimension, double>::internal_reinit<double>(
       const Mapping<deal_II_dimension> &,
       const std::vector<const hp::DoFHandler<deal_II_dimension> *> &,
       const std::vector<const AffineConstraints<double> *> &,
@@ -66,6 +68,22 @@ for (deal_II_dimension : DIMENSIONS)
       const std::vector<hp::QCollection<1>> &,
       const AdditionalData &);
 
+    template void MatrixFree<deal_II_dimension, float>::internal_reinit<float>(
+      const Mapping<deal_II_dimension> &,
+      const std::vector<const DoFHandler<deal_II_dimension> *> &,
+      const std::vector<const AffineConstraints<float> *> &,
+      const std::vector<IndexSet> &,
+      const std::vector<hp::QCollection<1>> &,
+      const AdditionalData &);
+
+    template void MatrixFree<deal_II_dimension, float>::internal_reinit<float>(
+      const Mapping<deal_II_dimension> &,
+      const std::vector<const hp::DoFHandler<deal_II_dimension> *> &,
+      const std::vector<const AffineConstraints<float> *> &,
+      const std::vector<IndexSet> &,
+      const std::vector<hp::QCollection<1>> &,
+      const AdditionalData &);
+
 #ifndef DEAL_II_MSVC
     template void
     internal::MatrixFreeFunctions::ShapeInfo<double>::reinit<deal_II_dimension>(
-- 
2.39.5