namespace VectorTools
{
- template <int dim,
- int spacedim,
- typename VectorType,
- template <int, int> class DoFHandlerType>
- void interpolate(
- const Mapping<dim, spacedim> &mapping,
- const DoFHandlerType<dim, spacedim> &dof_handler,
- const Function<spacedim, typename VectorType::value_type> &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<dim, spacedim>::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 <int dim, int spacedim, typename number, typename T1, typename T2, typename T3>
+ void transform(const typename FiniteElementData<dim>::Conformity conformity,
+ const T1 &cell,
+ const unsigned int offset,
+ T2 &fe_values_jacobians,
+ T3 &function_values)
+ {
+ switch (conformity)
+ {
+ case FiniteElementData<dim>::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<number> 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<dim>::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<number> 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<dim>::H1:
+ DEAL_II_FALLTHROUGH;
+ case FiniteElementData<dim>::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<dim, spacedim> fe(dof_handler.get_fe());
- std::vector<types::global_dof_index> 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 <int dim, int spacedim, typename number, typename T1, typename T2, typename T3>
+ void apply_transform(const FiniteElement<dim, spacedim> &fe,
+ const T1 &cell,
+ unsigned int &offset, /* modifies offset */
+ T2 &fe_values_jacobians,
+ T3 &function_values)
+ {
+ if (const auto *system =
+ dynamic_cast<const FESystem<dim, spacedim> *>(&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<dim, spacedim, number>(base_fe,
+ cell,
+ offset,
+ fe_values_jacobians,
+ function_values);
+ }
+ }
+ }
+ else
+ {
+ transform<dim, spacedim, number>(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<std::vector<Vector<number> > > fe_function_values(fe.size());
- std::vector<std::vector<number> > fe_dof_values(fe.size());
+ // Internal implementation of interpolate that takes a generic functor
+ // function such that function(cell) is of type
+ // Function<spacedim, typename VectorType::value_type>
+ template <int dim,
+ int spacedim,
+ typename VectorType,
+ template <int, int> class DoFHandlerType,
+ typename T>
+ void interpolate(
+ const Mapping<dim, spacedim> &mapping,
+ const DoFHandlerType<dim, spacedim> &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<dim> support_quadrature;
- for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
- {
- const auto &points = fe[fe_index].get_generalized_support_points();
- support_quadrature.push_back(Quadrature<dim>(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<dim, spacedim> mapping_collection(mapping);
+ typedef typename VectorType::value_type number;
- hp::FEValues<dim, spacedim> 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<dim, spacedim> fe_values_jacobians(
- mapping_collection,
- fe,
- support_quadrature,
- update_jacobians | update_inverse_jacobians);
+ const hp::FECollection<dim, spacedim> fe(dof_handler.get_fe());
- //
- // Now loop over all locally owned, active cells.
- //
+ std::vector<types::global_dof_index> 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<std::vector<Vector<number> > > fe_function_values(fe.size());
+ std::vector<std::vector<number> > 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<Point<spacedim> > &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<number>(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<dim, spacedim>::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<dim>::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<dim> support_quadrature;
+ for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
{
- switch (conformity)
- {
- case FiniteElementData<dim>::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<number> 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<dim>(points));
+ }
- case FiniteElementData<dim>::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<dim, spacedim> 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<number> 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<dim, spacedim> fe_values(
+ mapping_collection,
+ fe,
+ support_quadrature,
+ update_quadrature_points);
- case FiniteElementData<dim>::H1:
- DEAL_II_FALLTHROUGH;
- case FiniteElementData<dim>::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<dim, spacedim> 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<Point<spacedim> > &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<number>(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<dim, spacedim, number>(
+ 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<dim, spacedim> &fe,
- unsigned int &offset) -> void
- {
- if (const auto *system =
- dynamic_cast<const FESystem<dim, spacedim> *>(&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<VectorType>::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<VectorType>::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<VectorType>::get(
+ vec, dofs_on_cell[i]);
+ ::dealii::internal::ElementAccess<VectorType>::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<VectorType>::get(interpolation, i);
+ const auto weight =
+ ::dealii::internal::ElementAccess<VectorType>::get(weights, i);
+ ::dealii::internal::ElementAccess<VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::get(
- vec, dofs_on_cell[i]);
- ::dealii::internal::ElementAccess<VectorType>::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<VectorType>::get(interpolation, i);
- const auto weight =
- ::dealii::internal::ElementAccess<VectorType>::get(weights, i);
- ::dealii::internal::ElementAccess<VectorType>::set(
- value / weight, i, vec);
- }
- vec.compress(VectorOperation::insert);
+
+ template <int dim,
+ int spacedim,
+ typename VectorType,
+ template <int, int> class DoFHandlerType>
+ void interpolate(
+ const Mapping<dim, spacedim> &mapping,
+ const DoFHandlerType<dim, spacedim> &dof_handler,
+ const Function<spacedim, typename VectorType::value_type> &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<dim, spacedim>::active_cell_iterator &)
+ -> const Function<spacedim, typename VectorType::value_type> &
+ {
+ return function;
+ };
+
+ interpolate(mapping, dof_handler, function_map, vec, component_mask);
}