From: Timo Heister Date: Sat, 18 Mar 2017 17:54:19 +0000 (-0400) Subject: FEFieldFunction behavior for nonlocal cells X-Git-Tag: v8.5.0-rc1~24^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8ca19a847eb44cb95643ba0557bb717e0f928378;p=dealii.git FEFieldFunction behavior for nonlocal cells This fixes a bug in FEFieldFunction in several functions were it tries to evaluate the function on an artificial cell and failing instead of throwing ExcPointNotAvailableHere like it should. Fix this. --- diff --git a/doc/news/changes/minor/20170318TimoHeister b/doc/news/changes/minor/20170318TimoHeister new file mode 100644 index 0000000000..cc9dd24b40 --- /dev/null +++ b/doc/news/changes/minor/20170318TimoHeister @@ -0,0 +1,3 @@ +Fixed: FEFieldFunction now correctly throws ExcPointNotAvailableHere if used in parallel with points not available locally. +
+(Timo Heister, 2017/03/18) diff --git a/include/deal.II/numerics/fe_field_function.templates.h b/include/deal.II/numerics/fe_field_function.templates.h index e30ff8df30..266bb74dab 100644 --- a/include/deal.II/numerics/fe_field_function.templates.h +++ b/include/deal.II/numerics/fe_field_function.templates.h @@ -85,6 +85,10 @@ namespace Functions cell_hint.get() = cell; + // check that the current cell is available: + AssertThrow (!cell->is_artificial(), + VectorTools::ExcPointNotAvailableHere()); + // Now we can find out about the point Quadrature quad(qp.get()); FEValues fe_v(mapping, cell->get_fe(), quad, @@ -136,6 +140,10 @@ namespace Functions qp = my_pair.second; } + // check that the current cell is available: + AssertThrow (!cell->is_artificial(), + VectorTools::ExcPointNotAvailableHere()); + cell_hint.get() = cell; // Now we can find out about the point @@ -204,6 +212,10 @@ namespace Functions qp = my_pair.second; } + // check that the current cell is available: + AssertThrow (!cell->is_artificial(), + VectorTools::ExcPointNotAvailableHere()); + cell_hint.get() = cell; // Now we can find out about the point @@ -480,6 +492,10 @@ namespace Functions point_flags[0] = true; } + // check that the cell is available: + AssertThrow (!cell->is_artificial(), + VectorTools::ExcPointNotAvailableHere()); + // Put in the first point. cells.push_back(cell); qpoints.push_back(std::vector >(1, qp.get())); diff --git a/tests/mpi/fe_field_function_03.cc b/tests/mpi/fe_field_function_03.cc new file mode 100644 index 0000000000..e4f5cf8e3c --- /dev/null +++ b/tests/mpi/fe_field_function_03.cc @@ -0,0 +1,210 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 2016 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 regression in FEFieldFunction that doesn't throw VectorTools::ExcPointNotAvailableHere() + +/* +An error occurred in line <3405> of file in function + void dealii::DoFCellAccessor, false>::get_dof_values(const InputVector &, ForwardIterator, ForwardIterator) const [DoFHandlerType = dealii::DoFHandler<2, 2>, lda = false, InputVector = dealii::TrilinosWrappers::MPI::Vector, ForwardIterator = double *] +The violated condition was: + this->is_artificial() == false +Additional information: + Can't ask for DoF indices on artificial cells. + */ + +#include "../tests.h" +#include +#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); + + GridGenerator::hyper_cube_with_cylindrical_hole (tr, + 0.05 /* cylinder radius */, + 0.4/2.0 /* box radius */); + tr.refine_global (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; + + typename Functions::FEFieldFunction, TrilinosWrappers::MPI::Vector> + field_func(dofh, x_rel); + + Point<2> p (0.1, 0.0); + std::vector > points; + points.push_back(p); + + for (int i=0;; ++i) + { + try + { + if (i==0) + { + deallog << "vector_value:" << std::endl; + Vector out(1); + field_func.vector_value(p, out); + deallog << " OK: " << out[0] << std::endl; + } + else if (i==1) + { + deallog << "value:" << std::endl; + double out = field_func.value(p); + deallog << " OK: " << out << std::endl; + } + else if (i==2) + { + deallog << "value_list:" << std::endl; + std::vector out(1, -42.0f); + field_func.value_list(points, out); + deallog << " OK: " << out[0] << std::endl; + } + else if (i==3) + { + deallog << "vector_value_list:" << std::endl; + std::vector > out(1, Vector(1)); + field_func.vector_value_list(points, out); + deallog << " OK: " << out[0][0] << std::endl; + } + else if (i==4) + { + deallog << "vector_gradient:" << std::endl; + std::vector > out(1, Tensor<1,2>()); + field_func.vector_gradient(p, out); + deallog << " OK: " << out[0] << std::endl; + } + else if (i==5) + { + deallog << "gradient:" << std::endl; + Tensor<1,2> out = field_func.gradient(p); + deallog << " OK: " << out[0] << std::endl; + } + else if (i==6) + { + deallog << "vector_gradient_list:" << std::endl; + std::vector > > out(1, std::vector >(1, Tensor<1,2>())); + field_func.vector_gradient_list(points, out); + deallog << " OK: " << out[0][0] << std::endl; + } + else if (i==7) + { + deallog << "gradient_list:" << std::endl; + std::vector > out(1, Tensor<1,2>()); + field_func.gradient_list(points, out); + deallog << " OK: " << out[0] << std::endl; + } + else if (i==8) + { + deallog << "laplacian:" << std::endl; + double out = field_func.laplacian(p); + deallog << " OK: " << out << std::endl; + } + else if (i==9) + { + deallog << "vector_laplacian:" << std::endl; + Vector out(1); + field_func.vector_laplacian(p, out); + deallog << " OK: " << out[0] << std::endl; + } + else if (i==10) + { + deallog << "laplacian_list:" << std::endl; + std::vector out(1); + field_func.laplacian_list(points, out); + deallog << " OK: " << out[0] << std::endl; + } + else if (i==11) + { + deallog << "vector_laplacian_list:" << std::endl; + std::vector > out(1, Vector(1)); + field_func.vector_laplacian_list(points, out); + deallog << " OK: " << out[0][0] << std::endl; + } + else + break; + } + catch (const typename Functions::FEFieldFunction,TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &) + { + deallog << " ExcPointNotAvailableHere" << std::endl; + } + catch (...) + { + deallog << " Oh no! Some other error that we shouldn't get." << std::endl; + } + } + + deallog << "OK" << std::endl; +} + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + deal_II_exceptions::disable_abort_on_exception(); + + test<2>(); +} diff --git a/tests/mpi/fe_field_function_03.with_trilinos=true.mpirun=3.output b/tests/mpi/fe_field_function_03.with_trilinos=true.mpirun=3.output new file mode 100644 index 0000000000..fb805ec861 --- /dev/null +++ b/tests/mpi/fe_field_function_03.with_trilinos=true.mpirun=3.output @@ -0,0 +1,80 @@ + +DEAL:0::vector_value: +DEAL:0:: OK: 2.10000 +DEAL:0::value: +DEAL:0:: OK: 2.10000 +DEAL:0::value_list: +DEAL:0:: OK: 2.10000 +DEAL:0::vector_value_list: +DEAL:0:: OK: 2.10000 +DEAL:0::vector_gradient: +DEAL:0:: OK: 1.00000 -3.55271e-15 +DEAL:0::gradient: +DEAL:0:: OK: 1.00000 +DEAL:0::vector_gradient_list: +DEAL:0:: OK: 1.00000 -3.55271e-15 +DEAL:0::gradient_list: +DEAL:0:: OK: 1.00000 -3.55271e-15 +DEAL:0::laplacian: +DEAL:0:: OK: 0.00000 +DEAL:0::vector_laplacian: +DEAL:0:: OK: 0.00000 +DEAL:0::laplacian_list: +DEAL:0:: OK: 0.00000 +DEAL:0::vector_laplacian_list: +DEAL:0:: OK: 0.00000 +DEAL:0::OK + +DEAL:1::vector_value: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::value: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::value_list: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::vector_value_list: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::vector_gradient: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::gradient: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::vector_gradient_list: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::gradient_list: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::laplacian: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::vector_laplacian: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::laplacian_list: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::vector_laplacian_list: +DEAL:1:: ExcPointNotAvailableHere +DEAL:1::OK + + +DEAL:2::vector_value: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::value: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::value_list: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::vector_value_list: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::vector_gradient: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::gradient: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::vector_gradient_list: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::gradient_list: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::laplacian: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::vector_laplacian: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::laplacian_list: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::vector_laplacian_list: +DEAL:2:: ExcPointNotAvailableHere +DEAL:2::OK +