From: Matthias Maier Date: Wed, 30 Aug 2017 01:41:47 +0000 (-0500) Subject: FESystem: Fix setup of generalized_support_points X-Git-Tag: v9.0.0-rc1~1140^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d545d14877b30506bb9007573fcc0a059c22c078;p=dealii.git FESystem: Fix setup of generalized_support_points --- 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> 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 const std::vector > & FiniteElement::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::initialize (const std::vectorn_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 ([&]()