From 5ea338d0cae46296ce61e0e6462308c2bf3c289a Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Wed, 27 Sep 2017 16:00:01 -0500 Subject: [PATCH] Reimplement VectorTools::interpolate_based_on_material_id --- .../deal.II/numerics/vector_tools.templates.h | 190 ++++-------------- 1 file changed, 36 insertions(+), 154 deletions(-) diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index ada229938a..e15e17fb70 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -248,8 +248,10 @@ namespace VectorTools // Internal implementation of interpolate that takes a generic functor - // function such that function(cell) is of type - // Function + // function such that function(cell).second is of type + // Function* and + // function(cell).first is a boolean indicating whether the current + // cell should be processed template > &generalized_support_points = @@ -390,11 +396,11 @@ namespace VectorTools dof_values.resize(n_dofs); // Get all function values: - Assert(n_components == function(cell).n_components, + Assert(n_components == function(cell).second->n_components, ExcDimensionMismatch(dof_handler.get_fe().n_components(), - function(cell).n_components)); - function(cell).vector_value_list(generalized_support_points, - function_values); + function(cell).second->n_components)); + function(cell).second->vector_value_list(generalized_support_points, + function_values); { // Before we can average, we have to transform all function values @@ -500,9 +506,9 @@ namespace VectorTools // internal implementation const auto function_map = [&function]( const typename DoFHandlerType::active_cell_iterator &) - -> const Function & + -> std::pair *> { - return function; + return std::make_pair(true, &function); }; interpolate(mapping, dof_handler, function_map, vec, component_mask); @@ -591,158 +597,34 @@ namespace VectorTools void interpolate_based_on_material_id (const Mapping &mapping, - const DoFHandlerType &dof, - const std::map *> &function_map, - VectorType &dst, + const DoFHandlerType &dof_handler, + const std::map *> &functions, + 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.") ); - - if ( function_map.size() == 0 ) - return; - - Assert( function_map.find(numbers::invalid_material_id) == function_map.end(), - ExcMessage("You cannot specify the invalid material indicator " - "in your function map.")); - - for (typename std::map* > - ::const_iterator - iter = function_map.begin(); - iter != function_map.end(); - ++iter ) - { - Assert( dof.get_fe(0).n_components() == iter->second->n_components, - ExcDimensionMismatch(dof.get_fe(0).n_components(), iter->second->n_components) ); - } - - 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(); - - std::vector< std::vector< Point > > unit_support_points(fe.size()); - for (unsigned int fe_index = 0; fe_index < fe.size(); ++fe_index) - { - unit_support_points[fe_index] = fe[fe_index].get_unit_support_points(); - Assert( unit_support_points[fe_index].size() != 0, - ExcNonInterpolatingFE() ); - } - - std::vector< std::vector > dofs_of_rep_points(fe.size()); - std::vector< 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 < fe.size(); ++fe_index) - { - for (unsigned int i = 0; i < fe[fe_index].dofs_per_cell; ++i) - { - bool representative = true; - - for (unsigned int j = dofs_of_rep_points[fe_index].size(); j > 0; --j) - if ( unit_support_points[fe_index][i] - == unit_support_points[fe_index][dofs_of_rep_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(dofs_of_rep_points[fe_index].size()); - dofs_of_rep_points[fe_index].push_back(i); - ++n_rep_points[fe_index]; - } - } - - Assert( dofs_of_rep_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() ); - } - - const unsigned int max_rep_points = *std::max_element(n_rep_points.begin(), - n_rep_points.end()); - std::vector< types::global_dof_index> dofs_on_cell(fe.max_dofs_per_cell()); - std::vector< Point > rep_points(max_rep_points); - - std::vector< std::vector > function_values_scalar(fe.size()); - std::vector< std::vector< Vector > > function_values_system(fe.size()); - - hp::QCollection support_quadrature; - for (unsigned int fe_index = 0; fe_index < fe.size(); ++fe_index) - support_quadrature.push_back( Quadrature(unit_support_points[fe_index]) ); - - 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() ) - if ( function_map.find(cell->material_id()) != function_map.end() ) - { - const unsigned int fe_index = cell->active_fe_index(); - - fe_values.reinit(cell); - - const std::vector< Point > &support_points - = fe_values.get_present_fe_values().get_quadrature_points(); - - rep_points.resize( dofs_of_rep_points[fe_index].size() ); - for (unsigned int i = 0; i < dofs_of_rep_points[fe_index].size(); ++i) - rep_points[i] = support_points[dofs_of_rep_points[fe_index][i]]; - - dofs_on_cell.resize( fe[fe_index].dofs_per_cell ); - cell->get_dof_indices(dofs_on_cell); - - if (fe_is_system) - { - function_values_system[fe_index].resize( n_rep_points[fe_index], - Vector(fe[fe_index].n_components()) ); - - function_map.find(cell->material_id())->second->vector_value_list(rep_points, - function_values_system[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; + 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.")); - if ( component_mask[component] ) - { - const unsigned int rep_dof = dof_to_rep_index_table[fe_index][i]; - ::dealii::internal::ElementAccess::set( - function_values_system[fe_index][rep_dof](component), - dofs_on_cell[i], dst); - } - } - } - else - { - function_values_scalar[fe_index].resize(n_rep_points[fe_index]); + Assert (vec.size() == dof_handler.n_dofs(), + ExcDimensionMismatch (vec.size(), dof_handler.n_dofs())); - function_map.find(cell->material_id())->second->value_list - (rep_points, - function_values_scalar[fe_index], - 0); + Assert (component_mask.n_selected_components(dof_handler.get_fe().n_components()) > 0, + ComponentMask::ExcNoComponentSelected()); - for (unsigned int i = 0; i < fe[fe_index].dofs_per_cell; ++i) - ::dealii::internal::ElementAccess::set( - function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]], - dofs_on_cell[i], dst); - } - } + // Create a small lambda capture wrapping the function map and call the + // internal implementation + const auto function_map = [&functions]( + const typename DoFHandlerType::active_cell_iterator &cell) + -> std::pair *> + { + const auto function = functions.find(cell->material_id()); + bool update = function != functions.end(); + return std::make_pair(update, function->second); + }; - dst.compress (VectorOperation::insert); + interpolate(mapping, dof_handler, function_map, vec, component_mask); } -- 2.39.5