From: Rene Gassmoeller Date: Tue, 24 Oct 2017 17:52:20 +0000 (-0600) Subject: Fix VectorTools::interpolate X-Git-Tag: v9.0.0-rc1~876^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5315%2Fhead;p=dealii.git Fix VectorTools::interpolate --- diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index bc81d25255..16cf85f6da 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -331,6 +331,10 @@ namespace VectorTools interpolation.reinit(vec); weights.reinit(vec); + // Store locally owned dofs, so that we can skip all non-local dofs, + // if they do not need to be interpolated. + const IndexSet locally_owned_dofs = vec.locally_owned_elements(); + // 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 @@ -424,11 +428,6 @@ namespace VectorTools for (unsigned int i=0; i < n_dofs; ++i) { - ::dealii::internal::ElementAccess::add( - typename VectorType::value_type(1.0), - dofs_on_cell[i], - weights); - const auto &nonzero_components = fe[fe_index].get_nonzero_components(i); @@ -443,18 +442,29 @@ namespace VectorTools if (selected) { // Add local values to the global vectors - ::dealii::internal::ElementAccess::add( - dof_values[i], dofs_on_cell[i], interpolation); + ::dealii::internal::ElementAccess::add(dof_values[i], + 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, simply copy all dof values - // from the vector "vec": - const auto value = - ::dealii::internal::ElementAccess::get( - vec, dofs_on_cell[i]); - ::dealii::internal::ElementAccess::add( - value, dofs_on_cell[i], interpolation); + // 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); + } } } } /* loop over dof_handler.active_cell_iterators() */ diff --git a/tests/mpi/interpolate_05.cc b/tests/mpi/interpolate_05.cc new file mode 100644 index 0000000000..7cedd689b1 --- /dev/null +++ b/tests/mpi/interpolate_05.cc @@ -0,0 +1,159 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// This is a modified copy from interpolate_01.cc. +// VectorTools::interpolate used to accidentally read from non-owned +// dofs when called with a component mask that excluded some components. +// This triggered an assertion. Test that the solution, +// namely skipping all non-local dofs that are not selected to be interpolated, works. + +#include "../tests.h" +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + + +template +void test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + // create a distributed mesh + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + GridGenerator::hyper_cube (tr); + tr.refine_global(2); + + // Create a system with more than one component + FESystem fe(FE_Q(1),2); + DoFHandler dofh(tr); + dofh.distribute_dofs (fe); + + IndexSet owned_set = dofh.locally_owned_dofs(); + TrilinosWrappers::MPI::Vector x; + + x.reinit(owned_set, MPI_COMM_WORLD); + + // Select the first of the two components + std::vector components(2); + components[0] = true; + + // This failed before this test was introduced. + // Interpolate onto the first component only. + VectorTools::interpolate(dofh, + Functions::ConstantFunction(1,2), + x, + ComponentMask(components)); + + // Integrate the difference in the first component, if everything went + // well, this should be zero. + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dofh, relevant_set); + TrilinosWrappers::MPI::Vector x_rel(relevant_set, MPI_COMM_WORLD); + x_rel = x; + Vector error (tr.n_active_cells()); + ComponentSelectFunction right_component_select(0,2); + ComponentSelectFunction wrong_component_select(1,2); + + VectorTools::integrate_difference (dofh, + x_rel, + Functions::ConstantFunction(1,2), + error, + QMidpoint(), + VectorTools::L2_norm, + &right_component_select); + + double norm = VectorTools::compute_global_error(tr, error, VectorTools::L2_norm); + + if (myid == 0) + deallog << dofh.n_locally_owned_dofs() << ' ' << dofh.n_dofs() + << std::endl + << "Error of interpolated component: " << norm + << std::endl; + + // Integrate the difference in the second component. Since we did not interpolate + // the function into the finite element space, this should be equal to the + // integral of the ConstantFunction (=1) over the domain (unit square/cube). + // Thus the integral should be one. + VectorTools::integrate_difference (dofh, + x_rel, + Functions::ConstantFunction(1,2), + error, + QMidpoint(), + VectorTools::L2_norm, + &wrong_component_select); + + norm = VectorTools::compute_global_error(tr, error, VectorTools::L2_norm); + + if (myid == 0) + deallog << "Error of not interpolated component: " << norm + << std::endl; + + +} + + +int main(int argc, char *argv[]) +{ +#ifdef DEAL_II_WITH_MPI + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); +#else + (void)argc; + (void)argv; + compile_time_error; +#endif + + unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + + + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + initlog(); + + deallog.push("2d"); + test<2>(); + deallog.pop(); + deallog.push("3d"); + test<3>(); + deallog.pop(); + } + else + { + test<2>(); + test<3>(); + } +} diff --git a/tests/mpi/interpolate_05.with_trilinos=true.mpirun=2.output b/tests/mpi/interpolate_05.with_trilinos=true.mpirun=2.output new file mode 100644 index 0000000000..02b0f1ab98 --- /dev/null +++ b/tests/mpi/interpolate_05.with_trilinos=true.mpirun=2.output @@ -0,0 +1,7 @@ + +DEAL:0:2d::30 50 +DEAL:0:2d::Error of interpolated component: 0.00000 +DEAL:0:2d::Error of not interpolated component: 1.00000 +DEAL:0:3d::150 250 +DEAL:0:3d::Error of interpolated component: 0.00000 +DEAL:0:3d::Error of not interpolated component: 1.00000