From 19898c76c1cf69f1bf9949ce511f773de4300f45 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Mon, 6 Jun 2016 16:44:56 +0100 Subject: [PATCH] fix FEFieldFunction in parallel This fixes ExcPointNotAvailableHere() thrown from several places inside FEFieldFunction, because we assert that the cell is locally owned. This is not true, because we can happen to look at a ghost cell. Of course it is totally fine to evaluate function values on ghost cells if we have a ghosted solution vector. Triggering all code paths in the test is somewhat tricky, because we need find_active_cell_around_point to return a ghost cell. --- .../numerics/fe_field_function.templates.h | 10 +- tests/mpi/fe_field_function_02.cc | 176 ++++++++++++++++++ ...tion_02.with_trilinos=true.mpirun=3.output | 7 + 3 files changed, 188 insertions(+), 5 deletions(-) create mode 100644 tests/mpi/fe_field_function_02.cc create mode 100644 tests/mpi/fe_field_function_02.with_trilinos=true.mpirun=3.output diff --git a/include/deal.II/numerics/fe_field_function.templates.h b/include/deal.II/numerics/fe_field_function.templates.h index aea385d4b6..82d36f4954 100644 --- a/include/deal.II/numerics/fe_field_function.templates.h +++ b/include/deal.II/numerics/fe_field_function.templates.h @@ -74,7 +74,7 @@ namespace Functions { const std::pair::type, Point > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, p); - AssertThrow (my_pair.first->is_locally_owned(), + AssertThrow (!my_pair.first->is_artificial(), VectorTools::ExcPointNotAvailableHere()); cell = my_pair.first; @@ -127,7 +127,7 @@ namespace Functions { const std::pair::type, Point > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, p); - AssertThrow (my_pair.first->is_locally_owned(), + AssertThrow (!my_pair.first->is_artificial(), VectorTools::ExcPointNotAvailableHere()); cell = my_pair.first; @@ -195,7 +195,7 @@ namespace Functions { const std::pair::type, Point > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, p); - AssertThrow (my_pair.first->is_locally_owned(), + AssertThrow (!my_pair.first->is_artificial(), VectorTools::ExcPointNotAvailableHere()); cell = my_pair.first; @@ -470,7 +470,7 @@ namespace Functions const std::pair::type, Point > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, points[0]); - AssertThrow (my_pair.first->is_locally_owned(), + AssertThrow (!my_pair.first->is_artificial(), VectorTools::ExcPointNotAvailableHere()); cell = my_pair.first; @@ -534,7 +534,7 @@ namespace Functions { const std::pair::type, Point > my_pair = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]); - AssertThrow (my_pair.first->is_locally_owned(), + AssertThrow (!my_pair.first->is_artificial(), VectorTools::ExcPointNotAvailableHere()); cells.push_back(my_pair.first); diff --git a/tests/mpi/fe_field_function_02.cc b/tests/mpi/fe_field_function_02.cc new file mode 100644 index 0000000000..aa6a476dbe --- /dev/null +++ b/tests/mpi/fe_field_function_02.cc @@ -0,0 +1,176 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2015 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. +// +// --------------------------------------------------------------------- + + + +// Test FEFieldFunction for parallel computations. Specifically, we should be able +// to evaluate it in owned or ghost cells. + +// This used to crash with VectorTools::ExcPointNotAvailableHere() because we +// didn't expect getting a ghost cell. + +#include "../tests.h" +#include "coarse_grid_common.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + + +template +class LinearFunction : public Function +{ +public: + double value (const Point &p, + const unsigned int) const + { + return p[0] + 2; + } +}; + + +template +void test() +{ + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + // the mesh needs to be irregular enough to throw off the "cell hint" and + // GridTools::find_active_cell_around_point() to return ghost cells + + GridGenerator::quarter_hyper_shell (tr, + Point(), + 0.5, + 1.0, + 0, + true); + static const SphericalManifold spherical_manifold; + tr.set_manifold (99, spherical_manifold); + + for (typename Triangulation::active_cell_iterator + cell = tr.begin_active(); + cell != tr.end(); ++cell) + cell->set_all_manifold_ids (99); + for (typename Triangulation::active_cell_iterator + cell = tr.begin_active(); + cell != tr.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f)) + cell->face(f)->set_all_manifold_ids (numbers::invalid_manifold_id); + static const HyperShellBoundary boundary_shell; + tr.set_boundary (0, boundary_shell); + tr.set_boundary (1, boundary_shell); + + tr.refine_global ((dim==2)?3:1); + + const FE_Q fe(1); + DoFHandler dofh(tr); + dofh.distribute_dofs (fe); + + TrilinosWrappers::MPI::Vector interpolated(dofh.locally_owned_dofs(), + MPI_COMM_WORLD); + VectorTools::interpolate (dofh, + LinearFunction(), + interpolated); + + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs (dofh, relevant_set); + TrilinosWrappers::MPI::Vector x_rel(relevant_set, MPI_COMM_WORLD); + x_rel = interpolated; + + Functions::FEFieldFunction,TrilinosWrappers::MPI::Vector> + field_function (dofh, x_rel); + + std::vector > points; + + for (typename Triangulation::active_cell_iterator cell=tr.begin_active(); + cell!=tr.end(); ++cell) + { + if (cell->is_artificial()) + continue; + + points.push_back(cell->center()); + + // throw off the current cell hint + field_function.set_active_cell(dofh.end()); + + { + Point p = cell->center(); + Assert(std::abs(field_function.value (p)-(p[0]+2.))<1e-10, ExcInternalError()); + + Tensor<1,dim> gradient = field_function.gradient (p); + Tensor<1,dim> exact_gradient; + exact_gradient[0]=1.0; + Assert((gradient-exact_gradient).norm()<1e-10, ExcInternalError()); + + Assert(std::abs(field_function.laplacian(p))<1e-10, ExcInternalError()); + } + + if (cell->is_locally_owned()) + { + // more evil: use a corner so we can end up in a different cell + Point p = cell->vertex(0); + Assert(std::abs(field_function.value (p)-(p[0]+2.))<1e-10, ExcInternalError()); + + Tensor<1,dim> gradient = field_function.gradient (p); + Tensor<1,dim> exact_gradient; + exact_gradient[0]=1.0; + Assert((gradient-exact_gradient).norm()<1e-10, ExcInternalError()); + + Assert(std::abs(field_function.laplacian(p))<1e-10, ExcInternalError()); + } + } + + // also hit the code path for list version (no need to cover gradients/laplace here): + std::vector values(points.size()); + field_function.value_list(points, values, 0); + + if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) + deallog << "OK" << std::endl; +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + + deallog.push("2d"); + test<2>(); + deallog.pop(); + + deallog.push("3d"); + test<3>(); + deallog.pop(); +} diff --git a/tests/mpi/fe_field_function_02.with_trilinos=true.mpirun=3.output b/tests/mpi/fe_field_function_02.with_trilinos=true.mpirun=3.output new file mode 100644 index 0000000000..9b4ff1aa57 --- /dev/null +++ b/tests/mpi/fe_field_function_02.with_trilinos=true.mpirun=3.output @@ -0,0 +1,7 @@ + +DEAL:0:2d::OK +DEAL:0:3d::OK + + + + -- 2.39.5