From: Matthias Maier Date: Sun, 17 Sep 2017 21:50:26 +0000 (-0500) Subject: VectorTools::interpolate - support transformations with FESystems X-Git-Tag: v9.0.0-rc1~1002^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6aecdc77d89b07756d6017834afaa3ff0c61c193;p=dealii.git VectorTools::interpolate - support transformations with FESystems --- diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index a553901b77..04966b4bff 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -20,8 +20,8 @@ #include #include #include -#include #include +#include #include #include #include @@ -71,6 +71,8 @@ #include #include +#include + #include #include #include @@ -234,74 +236,124 @@ namespace VectorTools function.vector_value_list(generalized_support_points, function_values); - // FIXME: In case of an FESystem we have to apply this - // transformation according to the conformity of each base 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. + // + // 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 - + // fe.conforming_space is not the right call (thing about BDM). + const auto transform = [&](const typename FiniteElementData::Conformity conformity, + unsigned int offset) + { + 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; - // 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. - switch (fe[fe_index].conforming_space) - { - 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(); + 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 - // value[m] <- sum jacobian_transpose[m][n] * old_value[n]: - const auto old_value = function_values[i]; - TensorAccessors::contract<1, 2, 1, dim>( - function_values[i], jacobians[i].transpose(), old_value); - } - break; + 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; - 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 + 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; - 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(); - - // value[m] <- sum inverse_jacobians[m][n] * old_value[n]: - const auto old_value = function_values[i]; - TensorAccessors::contract<1, 2, 1, dim>( - function_values[i], inverse_jacobians[i], old_value); - - for (unsigned int j = 0; j < n_components; ++j) - function_values[i][j] *= jacobians[i].determinant(); - } - break; + default: + // In case we deal with an unknown conformity, just assume we + // deal with a Lagrange element and do nothing. + break; - 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; + } /*switch*/ + }; - default: - // In case we deal with an unknown conformity, just assume we - // deal with a Lagrange element and do nothing. - break; + // 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. + if (const auto *system = + dynamic_cast *>(&fe[fe_index])) + { + // In case of an FESystem transform every (vector) component + // separately: - } /*switch*/ + unsigned int offset = 0; + for (unsigned int i = 0; i < system->n_base_elements(); ++i) + { + const auto &fe = system->base_element(i); + const auto multiplicity = system->element_multiplicity(i); + for (unsigned int m = 0; m < multiplicity; ++m) + { + transform(fe.conforming_space, offset); + offset += fe.n_components(); + } + } + } + else + { + transform(fe[fe_index].conforming_space, 0); + } FETools::convert_generalized_support_point_values_to_dof_values( fe[fe_index], function_values, dof_values);