From d545d14877b30506bb9007573fcc0a059c22c078 Mon Sep 17 00:00:00 2001
From: Matthias Maier <tamiko@43-1.org>
Date: Tue, 29 Aug 2017 20:41:47 -0500
Subject: [PATCH] FESystem: Fix setup of generalized_support_points

---
 include/deal.II/fe/fe_system.h | 14 ++++++++++
 source/fe/fe.cc                |  6 ++--
 source/fe/fe_system.cc         | 51 ++++++++++++++++++++++++++++++++++
 3 files changed, 67 insertions(+), 4 deletions(-)

diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h
index 7cec10c953..7ff99bbfb4 100644
--- a/include/deal.II/fe/fe_system.h
+++ b/include/deal.II/fe/fe_system.h
@@ -969,6 +969,20 @@ private:
       unsigned int> >
       base_elements;
 
+  /**
+   * An index table that maps generalized support points of a base element
+   * to the vector of generalized support points of the FE System.
+   * It holds true that
+   * @code
+   *   auto n = generalized_support_points_index_table[i][j];
+   *   generalized_support_points[n] ==
+   *           base_elements[i].generalized_support_points[j];
+   * @endcode
+   * for each base element (indexed by i) and each g. s. point of the base
+   * element (index by j).
+   */
+  std::vector<std::vector<std::size_t>> generalized_support_points_index_table;
+
   /**
    * This function is simply singled out of the constructors since there are
    * several of them. It sets up the index table for the system as well as @p
diff --git a/source/fe/fe.cc b/source/fe/fe.cc
index aecf5e86c7..c2e5179af8 100644
--- a/source/fe/fe.cc
+++ b/source/fe/fe.cc
@@ -979,10 +979,8 @@ template <int dim, int spacedim>
 const std::vector<Point<dim> > &
 FiniteElement<dim,spacedim>::get_generalized_support_points () const
 {
-  // a finite element may define
-  // support points, but only if
-  // there are as many as there are
-  // degrees of freedom
+  // If the finite element implements generalized support points, return
+  // those. Otherwise fall back to unit support points.
   return ((generalized_support_points.size() == 0)
           ? unit_support_points
           : generalized_support_points);
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index 40ee23cd78..10d8fd1c95 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -1587,6 +1587,57 @@ void FESystem<dim,spacedim>::initialize (const std::vector<const FiniteElement<d
         }
     });
 
+  // Initialize generalized support points and an (internal) index table
+  init_tasks += Threads::new_task ([&]()
+  {
+    // If one of the base elements has no generalized support points, then
+    // it makes no sense to define generalized support points for the
+    // composed element, so return an empty array to demonstrate that fact.
+    // Note that we ignore FE_Nothing in this logic.
+    for (unsigned int base_el=0; base_el < this->n_base_elements(); ++base_el)
+      if (!base_element(base_el).has_generalized_support_points() &&
+          base_element(base_el).dofs_per_cell != 0)
+        {
+          this->generalized_support_points.resize(0);
+          return;
+        }
+
+    // Iterate over all base elements, extract a representative set of
+    // _unique_ generalized support points and store the information how
+    // generalized support points of base elements are mapped to this list
+    // of representatives. Complexity O(n^2), where n is the number of
+    // generalized support points.
+
+    generalized_support_points_index_table.resize(this->n_base_elements());
+
+    for (unsigned int base_el = 0; base_el < this->n_base_elements(); ++base_el)
+      {
+        for (const auto &point :
+             base_element(base_el).get_generalized_support_points())
+          {
+            // Is point already an element of generalized_support_points?
+            const auto p =
+              std::find(std::begin(this->generalized_support_points),
+                        std::end(this->generalized_support_points),
+                        point);
+
+            if (p == std::end(this->generalized_support_points))
+              {
+                // If no, update the table and add the point to the vector
+                const auto n = this->generalized_support_points.size();
+                generalized_support_points_index_table[base_el].push_back(n);
+                this->generalized_support_points.push_back(point);
+              }
+            else
+              {
+                // If yes, just add the correct index to the table.
+                const auto n = p - std::begin(this->generalized_support_points);
+                generalized_support_points_index_table[base_el].push_back(n);
+              }
+          }
+      }
+  });
+
   // initialize quad dof index permutation in 3d and higher
   if (dim >= 3)
     init_tasks += Threads::new_task ([&]()
-- 
2.39.5