From: Peter Munch Date: Fri, 26 Jul 2024 06:44:49 +0000 (+0200) Subject: VT::interpolate() for multigrid levels X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=673f8c0f75605e1f69a6c80a772ef79d4317eb74;p=dealii.git VT::interpolate() for multigrid levels --- diff --git a/doc/news/changes/minor/20240726Munch b/doc/news/changes/minor/20240726Munch new file mode 100644 index 0000000000..05403034cc --- /dev/null +++ b/doc/news/changes/minor/20240726Munch @@ -0,0 +1,3 @@ +New: VectorTools::interpolate() can now also used on multigrid levels. +
+(Peter Munch, Richard Schussnig, 2024/07/26) diff --git a/include/deal.II/numerics/vector_tools_interpolate.h b/include/deal.II/numerics/vector_tools_interpolate.h index b4006dfbfa..684fd2b2d0 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.h +++ b/include/deal.II/numerics/vector_tools_interpolate.h @@ -77,7 +77,8 @@ namespace VectorTools const DoFHandler &dof, const Function &function, VectorType &vec, - const ComponentMask &component_mask = {}); + const ComponentMask &component_mask = {}, + const unsigned int level = numbers::invalid_unsigned_int); /** * Same as above but in an hp-context. @@ -91,7 +92,8 @@ namespace VectorTools const DoFHandler &dof, const Function &function, VectorType &vec, - const ComponentMask &component_mask = {}); + const ComponentMask &component_mask = {}, + const unsigned int level = numbers::invalid_unsigned_int); /** @@ -106,7 +108,8 @@ namespace VectorTools const DoFHandler &dof, const Function &function, VectorType &vec, - const ComponentMask &component_mask = {}); + const ComponentMask &component_mask = {}, + const unsigned int level = numbers::invalid_unsigned_int); /** * Interpolate different finite element spaces. The interpolation of vector diff --git a/include/deal.II/numerics/vector_tools_interpolate.templates.h b/include/deal.II/numerics/vector_tools_interpolate.templates.h index 96a617522e..f5fafb1488 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.templates.h +++ b/include/deal.II/numerics/vector_tools_interpolate.templates.h @@ -231,7 +231,8 @@ namespace VectorTools const DoFHandler &dof_handler, T &function, VectorType &vec, - const ComponentMask &component_mask) + const ComponentMask &component_mask, + const unsigned int level = numbers::invalid_unsigned_int) { Assert(component_mask.represents_n_components( dof_handler.get_fe_collection().n_components()), @@ -240,7 +241,10 @@ namespace VectorTools "zero or equal to the number of components in the finite " "element.")); - AssertDimension(vec.size(), dof_handler.n_dofs()); + if (level == numbers::invalid_unsigned_int) + AssertDimension(vec.size(), dof_handler.n_dofs()); + else + AssertDimension(vec.size(), dof_handler.n_dofs(level)); Assert(component_mask.n_selected_components( dof_handler.get_fe_collection().n_components()) > 0, @@ -354,143 +358,152 @@ namespace VectorTools // // Now loop over all locally owned, active cells. // + const auto runner = [&](const auto &cell) { + const unsigned int fe_index = cell->active_fe_index(); + + // Do nothing if there are no local degrees of freedom. + if (fe[fe_index].n_dofs_per_cell() == 0) + return; + + // Skip processing of the current cell if the function object is + // invalid. This is used by interpolate_by_material_id to skip + // interpolating over cells with unknown material id. + if (!function(cell)) + return; + + // Get transformed, generalized support points + fe_values.reinit(cell); + const std::vector> &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].n_dofs_per_cell(); + dofs_on_cell.resize(n_dofs); + cell->get_active_or_mg_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(); + // Only resize (and create sample entry) if sizes do not match + if (function_values.size() != generalized_support_points.size()) + function_values.resize(generalized_support_points.size(), + Vector(n_components)); + + // Get all function values: + AssertDimension(n_components, function(cell)->n_components); + function(cell)->vector_value_list(generalized_support_points, + function_values); + + // For the simple case with elements with support points, we will + // simply use the interpolated DoF values in the access loop further + // down. Otherwise, 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 (needs_expensive_algorithm[fe_index]) + { + dof_values.resize(n_dofs); + const unsigned int offset = + apply_transform(fe[fe_index], + /* starting_offset = */ 0, + fe_values, + function_values); + (void)offset; + Assert(offset == n_components, ExcInternalError()); + + FETools::convert_generalized_support_point_values_to_dof_values( + fe[fe_index], function_values, dof_values); + } - for (const 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].n_dofs_per_cell() == 0) - continue; - - // Skip processing of the current cell if the function object is - // invalid. This is used by interpolate_by_material_id to skip - // interpolating over cells with unknown material id. - if (!function(cell)) - continue; - - // Get transformed, generalized support points - fe_values.reinit(cell); - const std::vector> &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].n_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(); - // Only resize (and create sample entry) if sizes do not match - if (function_values.size() != generalized_support_points.size()) - function_values.resize(generalized_support_points.size(), - Vector(n_components)); - - // Get all function values: - AssertDimension(n_components, function(cell)->n_components); - function(cell)->vector_value_list(generalized_support_points, - function_values); - - // For the simple case with elements with support points, we will - // simply use the interpolated DoF values in the access loop further - // down. Otherwise, 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 (needs_expensive_algorithm[fe_index]) - { - dof_values.resize(n_dofs); - const unsigned int offset = - apply_transform(fe[fe_index], - /* starting_offset = */ 0, - fe_values, - function_values); - (void)offset; - Assert(offset == n_components, ExcInternalError()); - - 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) - { - 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 c = 0; c < nonzero_components.size(); ++c) - selected = - selected || (nonzero_components[c] && component_mask[c]); - - if (selected) - { + for (unsigned int i = 0; i < n_dofs; ++i) + { + 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 c = 0; c < nonzero_components.size(); ++c) + selected = + selected || (nonzero_components[c] && component_mask[c]); + + if (selected) + { #ifdef DEBUG - // make sure that all selected base elements are indeed - // interpolatory - - if (const auto fe_system = - dynamic_cast *>(&fe[fe_index])) - { - const auto index = - fe_system->system_to_base_index(i).first.first; - Assert(fe_system->base_element(index) - .has_generalized_support_points(), - ExcMessage("The component mask supplied to " - "VectorTools::interpolate selects a " - "non-interpolatory element.")); - } + // make sure that all selected base elements are indeed + // interpolatory + + if (const auto fe_system = + dynamic_cast *>(&fe[fe_index])) + { + const auto index = + fe_system->system_to_base_index(i).first.first; + Assert(fe_system->base_element(index) + .has_generalized_support_points(), + ExcMessage("The component mask supplied to " + "VectorTools::interpolate selects a " + "non-interpolatory element.")); + } #endif - // Add local values to the global vectors - if (needs_expensive_algorithm[fe_index]) - ::dealii::internal::ElementAccess::add( - dof_values[i], dofs_on_cell[i], interpolation); - else - { - const auto base_index = - fe[fe_index].system_to_base_index(i); - ::dealii::internal::ElementAccess::add( - function_values[base_index.second] - [base_index.first.second], - dofs_on_cell[i], - interpolation); - } + // Add local values to the global vectors + if (needs_expensive_algorithm[fe_index]) ::dealii::internal::ElementAccess::add( - typename VectorType::value_type(1.0), - dofs_on_cell[i], - weights); - } - else - { - // If a component is ignored, copy the dof values - // from the vector "vec", but only if they are locally - // available - if (locally_owned_dofs.is_element(dofs_on_cell[i])) - { - const auto value = - ::dealii::internal::ElementAccess::get( - vec, dofs_on_cell[i]); - ::dealii::internal::ElementAccess::add( - value, dofs_on_cell[i], interpolation); - ::dealii::internal::ElementAccess::add( - typename VectorType::value_type(1.0), - dofs_on_cell[i], - weights); - } - } + dof_values[i], dofs_on_cell[i], interpolation); + else + { + const auto base_index = + fe[fe_index].system_to_base_index(i); + ::dealii::internal::ElementAccess::add( + function_values[base_index.second] + [base_index.first.second], + dofs_on_cell[i], + interpolation); + } + ::dealii::internal::ElementAccess::add( + typename VectorType::value_type(1.0), + dofs_on_cell[i], + weights); + } + else + { + // If a component is ignored, copy the dof values + // from the vector "vec", but only if they are locally + // available + if (locally_owned_dofs.is_element(dofs_on_cell[i])) + { + const auto value = + ::dealii::internal::ElementAccess::get( + vec, dofs_on_cell[i]); + ::dealii::internal::ElementAccess::add( + value, dofs_on_cell[i], interpolation); + ::dealii::internal::ElementAccess::add( + typename VectorType::value_type(1.0), + dofs_on_cell[i], + weights); + } + } + } + }; + + if (level == numbers::invalid_unsigned_int) + { + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned()) + runner(cell); } - } /* loop over dof_handler.active_cell_iterators() */ + } + else + { + for (const auto &cell : dof_handler.mg_cell_iterators_on_level(level)) + if (cell->is_locally_owned_on_level()) + runner(cell); + } interpolation.compress(VectorOperation::add); weights.compress(VectorOperation::add); @@ -527,22 +540,21 @@ namespace VectorTools const DoFHandler &dof_handler, const Function &function, VectorType &vec, - const ComponentMask &component_mask) + const ComponentMask &component_mask, + const unsigned int level) { AssertDimension(dof_handler.get_fe_collection().n_components(), function.n_components); // Create a small lambda capture wrapping function and call the // internal implementation - const auto function_map = - [&function]( - const typename DoFHandler::active_cell_iterator &) + const auto function_map = [&function](const auto &) -> const Function * { return &function; }; internal::interpolate( - mapping, dof_handler, function_map, vec, component_mask); + mapping, dof_handler, function_map, vec, component_mask, level); } @@ -554,13 +566,15 @@ namespace VectorTools const DoFHandler &dof_handler, const Function &function, VectorType &vec, - const ComponentMask &component_mask) + const ComponentMask &component_mask, + const unsigned int level) { interpolate(hp::MappingCollection(mapping), dof_handler, function, vec, - component_mask); + component_mask, + level); } @@ -571,7 +585,8 @@ namespace VectorTools const DoFHandler &dof, const Function &function, VectorType &vec, - const ComponentMask &component_mask) + const ComponentMask &component_mask, + const unsigned int level) { AssertDimension(dof.get_fe_collection().n_components(), function.n_components); @@ -579,7 +594,8 @@ namespace VectorTools dof, function, vec, - component_mask); + component_mask, + level); } diff --git a/source/numerics/vector_tools_interpolate.inst.in b/source/numerics/vector_tools_interpolate.inst.in index d0da45e0e8..b79665303f 100644 --- a/source/numerics/vector_tools_interpolate.inst.in +++ b/source/numerics/vector_tools_interpolate.inst.in @@ -28,7 +28,8 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const DoFHandler &, const Function &, VEC &, - const ComponentMask &); + const ComponentMask &, + const unsigned int); template void interpolate( @@ -36,14 +37,16 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; const DoFHandler &, const Function &, VEC &, - const ComponentMask &); + const ComponentMask &, + const unsigned int); template void interpolate( const DoFHandler &, const Function &, VEC &, - const ComponentMask &); + const ComponentMask &, + const unsigned int); template void interpolate( diff --git a/tests/multigrid/mg_data_out_04.cc b/tests/multigrid/mg_data_out_04.cc index 3b31dc3ba7..a73eff6766 100644 --- a/tests/multigrid/mg_data_out_04.cc +++ b/tests/multigrid/mg_data_out_04.cc @@ -36,7 +36,7 @@ template void -do_test() +do_test(const bool interpolate_with_vector_tools) { parallel::distributed::Triangulation triangulation( MPI_COMM_WORLD, @@ -87,6 +87,19 @@ do_test() transfer.build(dof_handler); transfer.interpolate_to_mg(dof_handler, dof_vector, global_dof_vector); + if (interpolate_with_vector_tools) + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + { + dof_vector[level] = 0.0; + + VectorTools::interpolate(dof_handler, + Functions::SquareFunction(), + dof_vector[level], + {}, + level); + } + for (unsigned int level = 0; level < triangulation.n_global_levels(); ++level) { @@ -115,7 +128,9 @@ main(int argc, char **argv) Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); MPILogInitAll log; - do_test<2>(); + do_test<2>(true); + do_test<2>(true); + // do_test<3>(); return 0; } diff --git a/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output b/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output index 568ac41070..c7bf970230 100644 --- a/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output +++ b/tests/multigrid/mg_data_out_04.with_p4est=true.mpirun=2.output @@ -100,7 +100,167 @@ DEAL:0::* level 2 0.500000 0.500000 0.500000 +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 +0.250000 0.00000 0.0625000 + +0.00000 0.250000 0.0625000 +0.250000 0.250000 0.125000 + + +0.250000 0.00000 0.0625000 +0.500000 0.00000 0.250000 + +0.250000 0.250000 0.125000 +0.500000 0.250000 0.312500 + + +0.00000 0.250000 0.0625000 +0.250000 0.250000 0.125000 + +0.00000 0.500000 0.250000 +0.250000 0.500000 0.312500 + + +0.250000 0.250000 0.125000 +0.500000 0.250000 0.312500 + +0.250000 0.500000 0.312500 +0.500000 0.500000 0.500000 + + +DEAL:0::* level 0 +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 +1.00000 0.00000 1.00000 + +0.00000 1.00000 1.00000 +1.00000 1.00000 2.00000 + + +DEAL:0::* level 1 +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 +0.500000 0.00000 0.250000 + +0.00000 0.500000 0.250000 +0.500000 0.500000 0.500000 + + +DEAL:0::* level 2 +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.00000 0.00000 +0.250000 0.00000 0.0625000 + +0.00000 0.250000 0.0625000 +0.250000 0.250000 0.125000 + + +0.250000 0.00000 0.0625000 +0.500000 0.00000 0.250000 + +0.250000 0.250000 0.125000 +0.500000 0.250000 0.312500 + + +0.00000 0.250000 0.0625000 +0.250000 0.250000 0.125000 + +0.00000 0.500000 0.250000 +0.250000 0.500000 0.312500 + + +0.250000 0.250000 0.125000 +0.500000 0.250000 0.312500 + +0.250000 0.500000 0.312500 +0.500000 0.500000 0.500000 + + + +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.500000 0.00000 0.250000 +1.00000 0.00000 1.00000 + +0.500000 0.500000 0.500000 +1.00000 0.500000 1.25000 + + +0.00000 0.500000 0.250000 +0.500000 0.500000 0.500000 + +0.00000 1.00000 1.00000 +0.500000 1.00000 1.25000 + + +0.500000 0.500000 0.500000 +1.00000 0.500000 1.25000 + +0.500000 1.00000 1.25000 +1.00000 1.00000 2.00000 + + +DEAL:1::* level 0 +DEAL:1::* level 1 +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.500000 0.00000 0.250000 +1.00000 0.00000 1.00000 + +0.500000 0.500000 0.500000 +1.00000 0.500000 1.25000 + + +0.00000 0.500000 0.250000 +0.500000 0.500000 0.500000 + +0.00000 1.00000 1.00000 +0.500000 1.00000 1.25000 + + +0.500000 0.500000 0.500000 +1.00000 0.500000 1.25000 +0.500000 1.00000 1.25000 +1.00000 1.00000 2.00000 + + +DEAL:1::* level 2 # This file was generated by the deal.II library.