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 (file)
index 0000000..cc9dd24
--- /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.
+<br>
+(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 e30ff8df306bb14dc6a0a79974015567fbc0bbd4..266bb74dab623f95f762ad9ec93df6da40bf9412 100644 (file)
--- 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<dim> quad(qp.get());
     FEValues<dim> 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<Point<dim> >(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 (file)
index 0000000..e4f5cf8
--- /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 </ssd/deal-git/include/deal.II/dofs/dof_accessor.templates.h> in function
+    void dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, 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 <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+template <int dim>
+class LinearFunction : public Function<dim>
+{
+public:
+  double value (const Point<dim> &p,
+                const unsigned int) const
+  {
+    return p[0] + 2;
+  }
+};
+
+
+template<int dim>
+void test()
+{
+  parallel::distributed::Triangulation<dim> 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<dim> fe(1);
+  DoFHandler<dim> dofh(tr);
+  dofh.distribute_dofs (fe);
+
+  TrilinosWrappers::MPI::Vector interpolated(dofh.locally_owned_dofs(),
+                                             MPI_COMM_WORLD);
+  VectorTools::interpolate (dofh,
+                            LinearFunction<dim>(),
+                            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<dim, DoFHandler<dim>, TrilinosWrappers::MPI::Vector>
+  field_func(dofh, x_rel);
+
+  Point<2> p (0.1, 0.0);
+  std::vector<Point<2> > points;
+  points.push_back(p);
+
+  for (int i=0;; ++i)
+    {
+      try
+        {
+          if (i==0)
+            {
+              deallog << "vector_value:" << std::endl;
+              Vector<double> 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<double> 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<Vector<double> > out(1, Vector<double>(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<Tensor<1,2> > 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<std::vector<Tensor<1,2> > > out(1, std::vector<Tensor<1,2> >(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<Tensor<1,2> > 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<double> 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<double> 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<Vector<double> > out(1, Vector<double>(1));
+              field_func.vector_laplacian_list(points, out);
+              deallog << "  OK: " << out[0][0] << std::endl;
+            }
+          else
+            break;
+        }
+      catch (const typename Functions::FEFieldFunction<dim,DoFHandler<dim>,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 (file)
index 0000000..fb805ec
--- /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
+