From: Matthias Maier Date: Wed, 27 Sep 2017 17:12:59 +0000 (-0500) Subject: refactor code X-Git-Tag: v9.0.0-rc1~1002^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ca0105805b8b835d11b624dfe1cd73cf08d623db;p=dealii.git refactor code --- diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index b1d7eff7e6..2e68509b63 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -87,344 +87,427 @@ DEAL_II_NAMESPACE_OPEN namespace VectorTools { - template class DoFHandlerType> - void interpolate( - const Mapping &mapping, - const DoFHandlerType &dof_handler, - const Function &function, - VectorType &vec, - const ComponentMask &component_mask) + + // This anonymous namespace contains the actual implementation called + // by VectorTools::interpolate and variants (such as + // VectorTools::interpolate_by_material_id). + namespace { - Assert(component_mask.represents_n_components(dof_handler.get_fe().n_components()), - ExcMessage("The number of components in the mask has to be either " - "zero or equal to the number of components in the finite " - "element.")); + // A small helper function to transform a component range starting + // at offset from the real to the unit cell according to the + // supplied conformity. The function_values vector is transformed + // in place. + // + // FIXME: This should be refactored into the mapping (i.e. + // implement the inverse function of Mapping::transform). + // Further, the finite element should make the information about + // the correct mapping directly accessible (i.e. which MappingType + // should be used). Using fe.conforming_space might be a bit of a + // problem because we only support doing nothing, Hcurl, and Hdiv + // conforming mappings. + // + // Input: + // conformity: conformity of the finite element, used to select + // appropriate type of transformation + // fe_values_jacobians, cell: used to reinitialize an fe_values object + // if values of jacbians (and inverses of + // jacobians) are needed + // function_values, offset: function_values is manipulated in place + // starting at position offset + template + void transform(const typename FiniteElementData::Conformity conformity, + const T1 &cell, + const unsigned int offset, + T2 &fe_values_jacobians, + T3 &function_values) + { + switch (conformity) + { + case FiniteElementData::Hcurl: + // See Monk, Finite Element Methods for Maxwell's Equations, + // p. 77ff, formula (3.76) and Corollary 3.58. + // For given mapping F_K: \hat K \to K, we have to transform + // \hat u = (dF_K)^T u\circ F_K + + fe_values_jacobians.reinit(cell); + for (unsigned int i = 0; i < function_values.size(); ++i) + { + const auto &jacobians = + fe_values_jacobians.get_present_fe_values().get_jacobians(); + + auto shifted_view = boost::make_iterator_range( + std::begin(function_values[i]) + offset, + std::begin(function_values[i]) + offset + dim); + std::vector old_value; + std::copy(std::begin(shifted_view), + std::end(shifted_view), + std::back_inserter(old_value)); + + // value[m] <- sum jacobian_transpose[m][n] * old_value[n]: + TensorAccessors::contract<1, 2, 1, dim>( + shifted_view, jacobians[i].transpose(), old_value); + } + break; - Assert (vec.size() == dof_handler.n_dofs(), - ExcDimensionMismatch (vec.size(), dof_handler.n_dofs())); + case FiniteElementData::Hdiv: + // See Monk, Finite Element Methods for Maxwell's Equations, + // p. 79ff, formula (3.77) and Lemma 3.59. + // For given mapping F_K: \hat K \to K, we have to transform + // \hat w = det(dF_K) (dF_K)^{-1} w\circ F_K - Assert (dof_handler.get_fe().n_components() == function.n_components, - ExcDimensionMismatch(dof_handler.get_fe().n_components(), function.n_components)); + fe_values_jacobians.reinit(cell); + for (unsigned int i = 0; i < function_values.size(); ++i) + { + const auto &jacobians = + fe_values_jacobians.get_present_fe_values().get_jacobians(); + const auto &inverse_jacobians = + fe_values_jacobians.get_present_fe_values() + .get_inverse_jacobians(); + + auto shifted_view = boost::make_iterator_range( + std::begin(function_values[i]) + offset, + std::begin(function_values[i]) + offset + dim); + std::vector old_value; + std::copy(std::begin(shifted_view), + std::end(shifted_view), + std::back_inserter(old_value)); + + // value[m] <- sum inverse_jacobians[m][n] * old_value[n]: + TensorAccessors::contract<1, 2, 1, dim>( + shifted_view, inverse_jacobians[i], old_value); + + for (unsigned int j = 0; j < dim; ++j) + shifted_view[j] *= jacobians[i].determinant(); + } + break; - Assert (component_mask.n_selected_components(dof_handler.get_fe().n_components()) > 0, - ComponentMask::ExcNoComponentSelected()); + case FiniteElementData::H1: + DEAL_II_FALLTHROUGH; + case FiniteElementData::L2: + // See Monk, Finite Element Methods for Maxwell's Equations, + // p. 77ff, formula (3.74). + // For given mapping F_K: \hat K \to K, we have to transform + // \hat p = p\circ F_K + // i.e., do nothing. + break; - // - // Computing the generalized interpolant isn't quite as straightforward - // as for classical Lagrange elements. A major complication is the fact - // it generally doesn't hold true that a function evaluates to the same - // dof coefficient on different cells. This means *setting* the value - // of a (global) degree of freedom computed on one cell doesn't - // necessarily lead to the same result when computed on a neighboring - // cell (that shares the same global degree of freedom). - // - // We thus, do the following operation: - // - // On each cell: - // - // - We first determine all function values u(x_i) in generalized - // support points - // - // - We transform these function values back to the unit cell - // according to the conformity of the component (scalar, Hdiv, or - // Hcurl conforming); see [Monk, Finite Element Methods for Maxwell's - // Equations, p.77ff Section 3.9] for details. This results in - // \hat u(\hat x_i) - // - // - We convert these generalized support point values to nodal values - // - // - For every global dof we take the average 1 / n_K \sum_{K} dof_K - // where n_K is the number of cells sharing the global dof and dof_K - // is the computed value on the cell K. - // - // For every degree of freedom that is shared by k cells, we compute - // its value on all k cells and take the weighted average with respect - // to the JxW values. - // + default: + // In case we deal with an unknown conformity, just assume we + // deal with a Lagrange element and do nothing. + break; - typedef typename VectorType::value_type number; + } /*switch*/ + } - const hp::FECollection fe(dof_handler.get_fe()); - std::vector dofs_on_cell(fe.max_dofs_per_cell()); + // A small helper function that iteratively applies above transform + // function to a vector function_values recursing over a given finite + // element decomposing it into base elements: + // + // Input + // fe: the full finite element corresponding to function_values + // [ rest see above] + template + void apply_transform(const FiniteElement &fe, + const T1 &cell, + unsigned int &offset, /* modifies offset */ + T2 &fe_values_jacobians, + T3 &function_values) + { + if (const auto *system = + dynamic_cast *>(&fe)) + { + // In case of an FESystem transform every (vector) component + // separately: + for (unsigned int i = 0; i < system->n_base_elements(); ++i) + { + const auto &base_fe = system->base_element(i); + const auto multiplicity = system->element_multiplicity(i); + for (unsigned int m = 0; m < multiplicity; ++m) + { + // recursively call apply_transform to make sure to + // correctly handle nested fe systems. + apply_transform(base_fe, + cell, + offset, + fe_values_jacobians, + function_values); + } + } + } + else + { + transform(fe.conforming_space, + cell, + offset, + fe_values_jacobians, + function_values); + offset += fe.n_components(); + } + }; - // Temporary storage for cell-wise interpolation operation. We store a - // variant for every fe we encounter to speed up resizing operations. - // The first vector is used for local function evaluation. The vector - // dof_values is used to store intermediate cell-wise interpolation - // results (see the detailed explanation in the for loop further down - // below). - std::vector > > fe_function_values(fe.size()); - std::vector > fe_dof_values(fe.size()); + // Internal implementation of interpolate that takes a generic functor + // function such that function(cell) is of type + // Function + template class DoFHandlerType, + typename T> + void interpolate( + const Mapping &mapping, + const DoFHandlerType &dof_handler, + T &function, + VectorType &vec, + const ComponentMask &component_mask) + { + Assert(component_mask.represents_n_components(dof_handler.get_fe().n_components()), + ExcMessage("The number of components in the mask has to be either " + "zero or equal to the number of components in the finite " + "element.")); - // We will need two temporary global vectors that store the new values - // and weights. - VectorType interpolation; - VectorType weights; - interpolation.reinit(vec); - weights.reinit(vec); + Assert (vec.size() == dof_handler.n_dofs(), + ExcDimensionMismatch (vec.size(), dof_handler.n_dofs())); - // We use an FEValues object to transform all generalized support - // points from the unit cell to the real cell coordinates. Thus, - // initialize a quadrature with all generalized support points and - // create an FEValues object with it. + Assert (component_mask.n_selected_components(dof_handler.get_fe().n_components()) > 0, + ComponentMask::ExcNoComponentSelected()); - hp::QCollection support_quadrature; - for (unsigned int fe_index=0; fe_index(points)); - } + // + // Computing the generalized interpolant isn't quite as straightforward + // as for classical Lagrange elements. A major complication is the fact + // it generally doesn't hold true that a function evaluates to the same + // dof coefficient on different cells. This means *setting* the value + // of a (global) degree of freedom computed on one cell doesn't + // necessarily lead to the same result when computed on a neighboring + // cell (that shares the same global degree of freedom). + // + // We thus, do the following operation: + // + // On each cell: + // + // - We first determine all function values u(x_i) in generalized + // support points + // + // - We transform these function values back to the unit cell + // according to the conformity of the component (scalar, Hdiv, or + // Hcurl conforming); see [Monk, Finite Element Methods for Maxwell's + // Equations, p.77ff Section 3.9] for details. This results in + // \hat u(\hat x_i) + // + // - We convert these generalized support point values to nodal values + // + // - For every global dof we take the average 1 / n_K \sum_{K} dof_K + // where n_K is the number of cells sharing the global dof and dof_K + // is the computed value on the cell K. + // + // For every degree of freedom that is shared by k cells, we compute + // its value on all k cells and take the weighted average with respect + // to the JxW values. + // - hp::MappingCollection mapping_collection(mapping); + typedef typename VectorType::value_type number; - hp::FEValues fe_values( - mapping_collection, - fe, - support_quadrature, - update_quadrature_points); - - // An extra FEValues object to compute jacobians. - // Only re-initialized in case of Hcurl or Hdiv conforming elements, - // i.e. if we really need the information. - hp::FEValues fe_values_jacobians( - mapping_collection, - fe, - support_quadrature, - update_jacobians | update_inverse_jacobians); + const hp::FECollection fe(dof_handler.get_fe()); - // - // Now loop over all locally owned, active cells. - // + std::vector dofs_on_cell(fe.max_dofs_per_cell()); - for (auto cell : dof_handler.active_cell_iterators()) - { - // If this cell is not locally owned, do nothing. - if (! cell->is_locally_owned()) - continue; + // Temporary storage for cell-wise interpolation operation. We store a + // variant for every fe we encounter to speed up resizing operations. + // The first vector is used for local function evaluation. The vector + // dof_values is used to store intermediate cell-wise interpolation + // results (see the detailed explanation in the for loop further down + // below). - const unsigned int fe_index = cell->active_fe_index(); + std::vector > > fe_function_values(fe.size()); + std::vector > fe_dof_values(fe.size()); - // Do nothing if there are no local degrees of freedom. - if (fe[fe_index].dofs_per_cell == 0) - continue; + // We will need two temporary global vectors that store the new values + // and weights. + VectorType interpolation; + VectorType weights; + interpolation.reinit(vec); + weights.reinit(vec); - // Get transformed, generalized support points - fe_values.reinit(cell); - const std::vector > &generalized_support_points = - fe_values.get_present_fe_values().get_quadrature_points(); - - // Get indices of the dofs on this cell - const auto n_dofs = fe[fe_index].dofs_per_cell; - dofs_on_cell.resize (n_dofs); - cell->get_dof_indices (dofs_on_cell); - - // Prepare temporary storage - auto &function_values = fe_function_values[fe_index]; - auto &dof_values = fe_dof_values[fe_index]; - - const auto n_components = fe[fe_index].n_components(); - function_values.resize(generalized_support_points.size(), - Vector(n_components)); - dof_values.resize(n_dofs); - - // Get all function values: - function.vector_value_list(generalized_support_points, - function_values); - - // A small helper function to transform a component range starting - // at offset from the real to the unit cell according to the - // supplied conformity. The function_values vector is transformed - // in place. - // - // FIXME: This should be refactored into the mapping (i.e. - // implement the inverse function of Mapping::transform). - // Further, the finite element should make the information about - // the correct mapping directly accessible (i.e. which MappingType - // should be used). Using fe.conforming_space might be a bit of a - // problem because we only support doing nothing, Hcurl, and Hdiv - // conforming mappings. - - const auto transform = [&function_values, &fe_values_jacobians, &cell]( - const typename FiniteElementData::Conformity conformity, - const unsigned int offset) + // We use an FEValues object to transform all generalized support + // points from the unit cell to the real cell coordinates. Thus, + // initialize a quadrature with all generalized support points and + // create an FEValues object with it. + + hp::QCollection support_quadrature; + for (unsigned int fe_index=0; fe_index::Hcurl: - // See Monk, Finite Element Methods for Maxwell's Equations, - // p. 77ff, formula (3.76) and Corollary 3.58. - // For given mapping F_K: \hat K \to K, we have to transform - // \hat u = (dF_K)^T u\circ F_K - - fe_values_jacobians.reinit(cell); - for (unsigned int i = 0; i < function_values.size(); ++i) - { - const auto &jacobians = - fe_values_jacobians.get_present_fe_values() - .get_jacobians(); - - auto shifted_view = boost::make_iterator_range( - std::begin(function_values[i]) + offset, - std::begin(function_values[i]) + offset + dim); - std::vector old_value; - std::copy(std::begin(shifted_view), - std::end(shifted_view), - std::back_inserter(old_value)); - - // value[m] <- sum jacobian_transpose[m][n] * old_value[n]: - TensorAccessors::contract<1, 2, 1, dim>( - shifted_view, jacobians[i].transpose(), old_value); - } - break; + const auto &points = fe[fe_index].get_generalized_support_points(); + support_quadrature.push_back(Quadrature(points)); + } - case FiniteElementData::Hdiv: - // See Monk, Finite Element Methods for Maxwell's Equations, - // p. 79ff, formula (3.77) and Lemma 3.59. - // For given mapping F_K: \hat K \to K, we have to transform - // \hat w = det(dF_K) (dF_K)^{-1} w\circ F_K + hp::MappingCollection mapping_collection(mapping); - fe_values_jacobians.reinit(cell); - for (unsigned int i = 0; i < function_values.size(); ++i) - { - const auto &jacobians = - fe_values_jacobians.get_present_fe_values().get_jacobians(); - const auto &inverse_jacobians = - fe_values_jacobians.get_present_fe_values().get_inverse_jacobians(); - - auto shifted_view = boost::make_iterator_range( - std::begin(function_values[i]) + offset, - std::begin(function_values[i]) + offset + dim); - std::vector old_value; - std::copy(std::begin(shifted_view), - std::end(shifted_view), - std::back_inserter(old_value)); - - // value[m] <- sum inverse_jacobians[m][n] * old_value[n]: - TensorAccessors::contract<1, 2, 1, dim>( - shifted_view, inverse_jacobians[i], old_value); - - for (unsigned int j = 0; j < dim; ++j) - shifted_view[j] *= jacobians[i].determinant(); - } - break; + hp::FEValues fe_values( + mapping_collection, + fe, + support_quadrature, + update_quadrature_points); - case FiniteElementData::H1: - DEAL_II_FALLTHROUGH; - case FiniteElementData::L2: - // See Monk, Finite Element Methods for Maxwell's Equations, - // p. 77ff, formula (3.74). - // For given mapping F_K: \hat K \to K, we have to transform - // \hat p = p\circ F_K - // i.e., do nothing. - break; + // An extra FEValues object to compute jacobians. + // Only re-initialized in case of Hcurl or Hdiv conforming elements, + // i.e. if we really need the information. + hp::FEValues fe_values_jacobians( + mapping_collection, + fe, + support_quadrature, + update_jacobians | update_inverse_jacobians); - default: - // In case we deal with an unknown conformity, just assume we - // deal with a Lagrange element and do nothing. - break; + // + // Now loop over all locally owned, active cells. + // - } /*switch*/ - }; /* lambda function transform */ + for (auto cell : dof_handler.active_cell_iterators()) + { + // If this cell is not locally owned, do nothing. + if (! cell->is_locally_owned()) + continue; + + const unsigned int fe_index = cell->active_fe_index(); + + // Do nothing if there are no local degrees of freedom. + if (fe[fe_index].dofs_per_cell == 0) + continue; + + // Get transformed, generalized support points + fe_values.reinit(cell); + const std::vector > &generalized_support_points = + fe_values.get_present_fe_values().get_quadrature_points(); + + // Get indices of the dofs on this cell + const auto n_dofs = fe[fe_index].dofs_per_cell; + dofs_on_cell.resize (n_dofs); + cell->get_dof_indices (dofs_on_cell); + + // Prepare temporary storage + auto &function_values = fe_function_values[fe_index]; + auto &dof_values = fe_dof_values[fe_index]; + + const auto n_components = fe[fe_index].n_components(); + function_values.resize(generalized_support_points.size(), + Vector(n_components)); + dof_values.resize(n_dofs); + + // Get all function values: + Assert(n_components == function(cell).n_components, + ExcDimensionMismatch(dof_handler.get_fe().n_components(), + function(cell).n_components)); + function(cell).vector_value_list(generalized_support_points, + function_values); - // Before we can average, we have to transform all function values - // from the real cell back to the unit cell. We query the finite - // element for the correct transformation. Matters get a bit more - // complicated because we have to apply said transformation for - // every base element. + { + // Before we can average, we have to transform all function values + // from the real cell back to the unit cell. We query the finite + // element for the correct transformation. Matters get a bit more + // complicated because we have to apply said transformation for + // every base element. + + unsigned int offset = 0; + apply_transform( + fe[fe_index], cell, offset, fe_values_jacobians, function_values); + Assert(offset == n_components, ExcInternalError()); + } - // modifies offset - const auto apply_transformation = [&](auto &&self, - const FiniteElement &fe, - unsigned int &offset) -> void - { - if (const auto *system = - dynamic_cast *>(&fe)) + FETools::convert_generalized_support_point_values_to_dof_values( + fe[fe_index], function_values, dof_values); + + for (unsigned int i=0; i < n_dofs; ++i) { - // In case of an FESystem transform every (vector) component - // separately: - for (unsigned int i = 0; i < system->n_base_elements(); ++i) + ::dealii::internal::ElementAccess::add( + typename VectorType::value_type(1.0), + dofs_on_cell[i], + weights); + + const auto &nonzero_components = + fe[fe_index].get_nonzero_components(i); + + // Figure out whether the component mask applies. We assume + // that we are allowed to set degrees of freedom if at least + // one of the components (of the dof) is selected. + bool selected = false; + for (unsigned int i = 0; i < nonzero_components.size(); ++i) + selected = + selected || (nonzero_components[i] && component_mask[i]); + + if (selected) { - const auto &base_fe = system->base_element(i); - const auto multiplicity = system->element_multiplicity(i); - for (unsigned int m = 0; m < multiplicity; ++m) - { - // recursively call apply_transform to make sure to - // correctly handle nested fe systems. - self(self, base_fe, offset); - } + // Add local values to the global vectors + ::dealii::internal::ElementAccess::add( + dof_values[i], dofs_on_cell[i], interpolation); + } + else + { + // If a component is ignored, simply copy all dof values + // from the vector "vec": + const auto value = + ::dealii::internal::ElementAccess::get( + vec, dofs_on_cell[i]); + ::dealii::internal::ElementAccess::add( + value, dofs_on_cell[i], interpolation); } } - else - { - transform(fe.conforming_space, offset); - offset += fe.n_components(); - } - }; + } /* loop over dof_handler.active_cell_iterators() */ + interpolation.compress(VectorOperation::add); + weights.compress(VectorOperation::add); + + for (const auto i : interpolation.locally_owned_elements()) { - unsigned int offset = 0; - apply_transformation(apply_transformation, fe[fe_index], offset); - Assert(offset == fe[fe_index].n_components(), ExcInternalError()); + const auto value = + ::dealii::internal::ElementAccess::get(interpolation, i); + const auto weight = + ::dealii::internal::ElementAccess::get(weights, i); + ::dealii::internal::ElementAccess::set( + value / weight, i, vec); } + vec.compress(VectorOperation::insert); + } - FETools::convert_generalized_support_point_values_to_dof_values( - fe[fe_index], function_values, dof_values); + } /* internal namespace */ - for (unsigned int i=0; i < n_dofs; ++i) - { - ::dealii::internal::ElementAccess::add( - typename VectorType::value_type(1.0), - dofs_on_cell[i], - weights); - - const auto &nonzero_components = - fe[fe_index].get_nonzero_components(i); - - // Figure out whether the component mask applies. We assume - // that we are allowed to set degrees of freedom if at least - // one of the components (of the dof) is selected. - bool selected = false; - for (unsigned int i = 0; i < nonzero_components.size(); ++i) - selected = - selected || (nonzero_components[i] && component_mask[i]); - - if (selected) - { - // Add local values to the global vectors - ::dealii::internal::ElementAccess::add( - dof_values[i], dofs_on_cell[i], interpolation); - } - else - { - // If a component is ignored, simply copy all dof values - // from the vector "vec": - const auto value = - ::dealii::internal::ElementAccess::get( - vec, dofs_on_cell[i]); - ::dealii::internal::ElementAccess::add( - value, dofs_on_cell[i], interpolation); - } - } - } /* loop over dof_handler.active_cell_iterators() */ - interpolation.compress(VectorOperation::add); - weights.compress(VectorOperation::add); - for (const auto i : interpolation.locally_owned_elements()) - { - const auto value = - ::dealii::internal::ElementAccess::get(interpolation, i); - const auto weight = - ::dealii::internal::ElementAccess::get(weights, i); - ::dealii::internal::ElementAccess::set( - value / weight, i, vec); - } - vec.compress(VectorOperation::insert); + + template class DoFHandlerType> + void interpolate( + const Mapping &mapping, + const DoFHandlerType &dof_handler, + const Function &function, + VectorType &vec, + const ComponentMask &component_mask) + { + Assert(component_mask.represents_n_components(dof_handler.get_fe().n_components()), + ExcMessage("The number of components in the mask has to be either " + "zero or equal to the number of components in the finite " + "element.")); + + Assert (vec.size() == dof_handler.n_dofs(), + ExcDimensionMismatch (vec.size(), dof_handler.n_dofs())); + + Assert (dof_handler.get_fe().n_components() == function.n_components, + ExcDimensionMismatch(dof_handler.get_fe().n_components(), function.n_components)); + + Assert (component_mask.n_selected_components(dof_handler.get_fe().n_components()) > 0, + ComponentMask::ExcNoComponentSelected()); + + // Create a small lambda capture wrapping function and call the + // internal implementation + const auto function_map = [&function]( + const typename DoFHandlerType::active_cell_iterator &) + -> const Function & + { + return function; + }; + + interpolate(mapping, dof_handler, function_map, vec, component_mask); }