}
+ /**
+ * Older version of the function directly below that does the same thing. This
+ * version also works for FESystems that contain non-interpolatory elements
+ * (as long as their components are not selected in the component mask).
+ * TODO: Remove this function when the other function is fixed to also work
+ * for the mentioned case.
+ */
+ template <int dim, int spacedim, typename VectorType,
+ template <int, int> class DoFHandlerType, typename T>
+ void interpolate_selected_components (const Mapping<dim,spacedim> &mapping,
+ const DoFHandlerType<dim,spacedim> &dof,
+ const T &function,
+ VectorType &vec,
+ const ComponentMask &component_mask)
+ {
+ typedef typename VectorType::value_type number;
+ Assert (component_mask.represents_n_components(dof.get_fe(0).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.n_dofs(),
+ ExcDimensionMismatch (vec.size(), dof.n_dofs()));
+ Assert (component_mask.n_selected_components(dof.get_fe(0).n_components()) > 0,
+ ComponentMask::ExcNoComponentSelected());
+
+ const hp::FECollection<dim,spacedim> &fe = dof.get_fe_collection();
+ const unsigned int n_components = fe.n_components();
+ const bool fe_is_system = (n_components != 1);
+
+ typename DoFHandlerType<dim,spacedim>::active_cell_iterator
+ cell = dof.begin_active(),
+ endc = dof.end();
+
+ // For FESystems many of the
+ // unit_support_points will appear
+ // multiple times, as a point may be
+ // unit_support_point for several of the
+ // components of the system. The following
+ // is rather complicated, but at least
+ // attempts to avoid evaluating the
+ // vectorfunction multiple times at the
+ // same point on a cell.
+ //
+ // First check that the desired components are interpolating.
+ for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
+ {
+ for (unsigned int component_index = 0; component_index < n_components; ++component_index)
+ {
+ if (component_mask[component_index] == true)
+ {
+ Assert ((fe[fe_index].base_element(fe[fe_index].component_to_base_index
+ (component_index).first).has_support_points()) ||
+ (fe[fe_index].base_element(fe[fe_index].component_to_base_index
+ (component_index).first).dofs_per_cell == 0),
+ ExcNonInterpolatingFE());
+ }
+ }
+ }
+
+ // Find the support points on a cell that are mentioned multiple times, and
+ // only add each once. Each multiple point gets to know the dof index of
+ // its representative point by the dof_to_rep_dof_table.
+
+ // the following vector collects all unit support points p[i],
+ // 0<=i<fe.dofs_per_cell, for which unit_support_point(i) is unique
+ // (representative). The position of a support point within this vector is
+ // called the rep index.
+ std::vector<std::vector<Point<dim> > > rep_unit_support_points (fe.size());
+ // the following table converts a dof i
+ // to the rep index.
+ std::vector<std::vector<types::global_dof_index> > dof_to_rep_index_table(fe.size());
+
+ std::vector<unsigned int> n_rep_points (fe.size(), 0);
+
+ for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
+ {
+ for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
+ {
+ const unsigned int component
+ = fe[fe_index].system_to_component_index(i).first;
+ if (component_mask[component] == true)
+ {
+ const Point<dim> dof_support_point = fe[fe_index].unit_support_point(i);
+
+ bool representative=true;
+ // the following loop is looped
+ // the other way round to get
+ // the minimal effort of
+ // O(fe.dofs_per_cell) for multiple
+ // support points that are placed
+ // one after the other.
+ for (unsigned int j=rep_unit_support_points[fe_index].size(); j>0; --j)
+ if (dof_support_point
+ == rep_unit_support_points[fe_index][j-1])
+ {
+ dof_to_rep_index_table[fe_index].push_back(j-1);
+ representative=false;
+ break;
+ }
+
+ if (representative)
+ {
+ dof_to_rep_index_table[fe_index].push_back
+ (rep_unit_support_points[fe_index].size());
+ rep_unit_support_points[fe_index].push_back(dof_support_point);
+ ++n_rep_points[fe_index];
+ }
+ }
+ else
+ {
+ // If correct component not to be interpolated, append invalid index to in table
+ dof_to_rep_index_table[fe_index].push_back(numbers::invalid_dof_index);
+ }
+ }
+
+ Assert(rep_unit_support_points[fe_index].size()==n_rep_points[fe_index],
+ ExcInternalError());
+ Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell,
+ ExcInternalError());
+ }
+
+ std::vector<types::global_dof_index> dofs_on_cell (fe.max_dofs_per_cell());
+
+ // get space for the values of the
+ // function at the rep support points.
+ //
+ // have two versions, one for system fe
+ // and one for scalar ones, to take the
+ // more efficient one respectively
+ std::vector<std::vector<number> > function_values_scalar(fe.size());
+ std::vector<std::vector<Vector<number> > > function_values_system(fe.size());
+
+ // Make a quadrature rule from support points
+ // to feed it into FEValues
+ hp::QCollection<dim> support_quadrature;
+ for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
+ support_quadrature.push_back (Quadrature<dim>(rep_unit_support_points[fe_index]));
+
+ // Transformed support points are computed by
+ // FEValues
+ hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
+
+ hp::FEValues<dim,spacedim> fe_values (mapping_collection,
+ fe, support_quadrature, update_quadrature_points);
+
+ for (; cell!=endc; ++cell)
+ if (cell->is_locally_owned())
+ {
+ const unsigned int fe_index = cell->active_fe_index();
+ if (fe[fe_index].dofs_per_cell != 0)
+ {
+ // for each cell:
+ // get location of finite element
+ // support_points
+ fe_values.reinit(cell);
+ const std::vector<Point<spacedim> > &rep_support_points =
+ fe_values.get_present_fe_values().get_quadrature_points();
+
+ // get indices of the dofs on this cell
+ dofs_on_cell.resize (fe[fe_index].dofs_per_cell);
+ cell->get_dof_indices (dofs_on_cell);
+
+
+ if (fe_is_system)
+ {
+ // get function values at
+ // these points. Here: get
+ // all components
+ function_values_system[fe_index]
+ .resize (n_rep_points[fe_index],
+ Vector<number> (fe[fe_index].n_components()));
+
+ Assert (dof.get_fe(fe_index).n_components() == function(cell)->n_components,
+ ExcDimensionMismatch(dof.get_fe(0).n_components(),
+ function(cell)->n_components));
+
+ function(cell)->vector_value_list (rep_support_points,
+ function_values_system[fe_index]);
+ // distribute the function
+ // values to the global
+ // vector
+ for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
+ {
+ const unsigned int component
+ = fe[fe_index].system_to_component_index(i).first;
+ if (component_mask[component] == true)
+ {
+ const unsigned int rep_dof=dof_to_rep_index_table[fe_index][i];
+ ::dealii::internal::ElementAccess<VectorType>::set(
+ function_values_system[fe_index][rep_dof](component),
+ dofs_on_cell[i], vec);
+ }
+ }
+ }
+ else
+ {
+ // get first component only,
+ // which is the only component
+ // in the function anyway
+ function_values_scalar[fe_index].resize (n_rep_points[fe_index]);
+ function(cell)->value_list (rep_support_points,
+ function_values_scalar[fe_index],
+ 0);
+ // distribute the function
+ // values to the global
+ // vector
+ for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
+ ::dealii::internal::ElementAccess<VectorType>::set(
+ function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]],
+ dofs_on_cell[i], vec);
+ }
+ }
+ }
+ vec.compress(VectorOperation::insert);
+ }
+
+
// Internal implementation of interpolate that takes a generic functor
// function such that function(cell) is of type
// Function<spacedim, typename VectorType::value_type>*
//
// A given cell is skipped if function(cell) == nullptr
+ // TODO: Make this function also work for FESystems that
+ // contain elements without generalized support points, as
+ // long as their components are not selected for interpolation.
template <int dim,
int spacedim,
typename VectorType,
const hp::FECollection<dim, spacedim> &fe(dof_handler.get_fe_collection());
+ // Check if the element has generalized support points
+ bool has_generalized_support_points = true;
+ for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
+ {
+ has_generalized_support_points = has_generalized_support_points && fe[fe_index].has_generalized_support_points();
+ }
+
+ // If the element has no generalized support points the algorithm in this
+ // function does not work. Fall back to an earlier version.
+ if (!has_generalized_support_points)
+ {
+ interpolate_selected_components(mapping,dof_handler,function,vec,component_mask);
+ return;
+ }
+
std::vector<types::global_dof_index> dofs_on_cell(fe.max_dofs_per_cell());
// Temporary storage for cell-wise interpolation operation. We store a
projects / dealii.git / commitdiff