]> https://gitweb.dealii.org/ - dealii.git/commitdiff
FESystem: Fix setup of generalized_support_points
authorMatthias Maier <tamiko@43-1.org>
Wed, 30 Aug 2017 01:41:47 +0000 (20:41 -0500)
committerMatthias Maier <tamiko@43-1.org>
Thu, 31 Aug 2017 15:57:38 +0000 (10:57 -0500)
include/deal.II/fe/fe_system.h
source/fe/fe.cc
source/fe/fe_system.cc

index 7cec10c95303e7d8958f1a4ac5bde4cf6fa54a5c..7ff99bbfb4ef8d6d19d27ed43434df7b60720472 100644 (file)
@@ -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
index aecf5e86c79f95755624a31c3c05a8f32af029be..c2e5179af89a153aa0a52f42f7b5950765ad850a 100644 (file)
@@ -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);
index 40ee23cd78c586f992dab1172d0843a9b44adf03..10d8fd1c95245966723872b29f5a794a2e9dbf54 100644 (file)
@@ -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 ([&]()

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.