From 414f68206887bc359e12e20e2b709351a418a316 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 24 Nov 2017 10:26:07 -0600 Subject: [PATCH] Bugfix: Allow non-interpolatory elements in FESystem for generalized s.p. In case we encounter a non-interpolatory base element we now simply fill the corresponding degrees of freedom with NaNs. --- source/fe/fe_system.cc | 109 +++++++++++++++++++++++++---------------- 1 file changed, 68 insertions(+), 41 deletions(-) diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 575822c237..b004250988 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -1607,18 +1607,6 @@ 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 @@ -1629,6 +1617,18 @@ void FESystem::initialize (const std::vectorn_base_elements(); ++base) { + // If the current base element does not have generalized support + // points, ignore it. Note that + // * FESystem::convert_generalized_support_point_values_to_dof_values + // will simply skip such non-interpolatory base elements by + // assigning NaN to all dofs. + // * If this routine does not pick up any generalized support + // points the corresponding vector will be empty and + // FiniteElement::has_generalized_support_points will return + // false. + if (!base_element(base).has_generalized_support_points()) + continue; + for (const auto &point : base_element(base).get_generalized_support_points()) { @@ -1658,6 +1658,9 @@ void FESystem::initialize (const std::vectorget_generalized_support_points(); for (unsigned int j = 0; j < points.size(); ++j) @@ -2289,8 +2292,6 @@ convert_generalized_support_point_values_to_dof_values (const std::vectorbase_element(base); const unsigned int multiplicity = this->element_multiplicity(base); const unsigned int n_base_dofs = base_element.dofs_per_cell; - const unsigned int n_base_points = - base_element.get_generalized_support_points().size(); const unsigned int n_base_components = base_element.n_components(); // If the number of base degrees of freedom is zero, there is nothing @@ -2302,44 +2303,70 @@ convert_generalized_support_point_values_to_dof_values (const std::vectordofs_per_cell; ++i) + if (this->system_to_base_index(i).first == + std::make_pair(base, m)) + dof_values[i] = + base_dof_values[this->system_to_base_index(i).second]; + } /*for*/ + } + else + { + // If the base element is non-interpolatory, assign NaN to all + // DoFs associated to it. // To do this, we could really use a base_to_system_index() // function, but that doesn't exist -- just do it by using the // reverse table -- the amount of work done here is not worth // trying to optimizing this. - for (unsigned int i = 0; i < this->dofs_per_cell; ++i) - if (this->system_to_base_index(i).first == std::make_pair(base, m)) - dof_values[i] = - base_dof_values[this->system_to_base_index(i).second]; + for (unsigned int m = 0; m < multiplicity; ++m) + for (unsigned int i = 0; i < this->dofs_per_cell; ++i) + if (this->system_to_base_index(i).first == + std::make_pair(base, m)) + dof_values[i] = std::numeric_limits::signaling_NaN(); + + current_vector_component += multiplicity * n_base_components; } - } + } /*for*/ } -- 2.39.5