From 0b4e04b78a816308d18ca394fbc5640cd62f1377 Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Thu, 26 Oct 2017 09:35:14 -0600 Subject: [PATCH] Reinsert old function --- .../deal.II/numerics/vector_tools.templates.h | 236 ++++++++++++++++++ 1 file changed, 236 insertions(+) diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 16cf85f6da2b..aaf40160efd5 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -247,11 +247,232 @@ namespace VectorTools } + /** + * 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 class DoFHandlerType, typename T> + void interpolate_selected_components (const Mapping &mapping, + const DoFHandlerType &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 &fe = dof.get_fe_collection(); + const unsigned int n_components = fe.n_components(); + const bool fe_is_system = (n_components != 1); + + typename DoFHandlerType::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 > > rep_unit_support_points (fe.size()); + // the following table converts a dof i + // to the rep index. + std::vector > dof_to_rep_index_table(fe.size()); + + std::vector n_rep_points (fe.size(), 0); + + for (unsigned int fe_index=0; fe_index 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 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 > function_values_scalar(fe.size()); + std::vector > > function_values_system(fe.size()); + + // Make a quadrature rule from support points + // to feed it into FEValues + hp::QCollection support_quadrature; + for (unsigned int fe_index=0; fe_index(rep_unit_support_points[fe_index])); + + // Transformed support points are computed by + // FEValues + hp::MappingCollection mapping_collection (mapping); + + hp::FEValues 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 > &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 (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::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::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* // // 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 &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 dofs_on_cell(fe.max_dofs_per_cell()); // Temporary storage for cell-wise interpolation operation. We store a