Add parallel test for find_active_cell_around_point

[dangars]

The following test succeeds in serial (mpi = 1), but fails in parallel (mpi = 7). I think that it is related to #10478, but I'm not sure.

I get the following error with find_active_cell_around_point():

The point <0 0> could not be found inside any of the subcells of a coarse grid cell.

But the point <0 0> is in the coarse grid.

This is a follow-up from #10737

[peterrum]

The function find_active_cell_around_point only succeeds on locally owned and ghost cells (i.e. active cells). My guess is that the problematic point is in neither type of cells on one of the processes. The fact that it part of the coarse grid is not relevant. @dangars Could you check this?

FYI @nfehn

[dangars]

@peterrum That is right. The point is not in any of the active cells. But I do not understand why this worked before. In step-62 I used this function to find out if a point is in the locally owned cells.

According to the documentation, this function will raise an exception only if the point does not lie inside any of the cells of the mesh.

If the point requested does not lie in any of the cells of the mesh given,
then this function throws an exception of type GridTools::ExcPointNotFound.

I'm going to try to find out why the old version didn't raise an exception when the point does not lie in an active cell.

@peterrum Is find_active_cell_around_point() the right function to find out if a point is in a locally owned cell?

[peterrum]

We have changed recently the behavior of the function(s) because it turned out that working internally with non-active cells was not a good idea.

There is still a discussion how to deal with this changed behavior (see #10535). The documentation will be adjusted once we found a common ground.

Is find_active_cell_around_point() the right function to find out if a point is in a locally owned cell?

You can do that. But what you need to do is to catch the exception to be on the safe side.

[dangars]

Is find_active_cell_around_point() the right function to find out if a point is in a locally owned cell?

You can do that. But what you need to do is to catch the exception to be on the safe side.

Thanks @peterrum, I understand now :) I think that exceptions are quite inefficient for this. I will leave step-62 in a broken state until it has been decided what to do with this function. Then I will fix it.

[bangerth]

What's the status here?

[dangars]

I'll update this once #11276 is in master :)

[peterrum]

@dangars Could you rebase this PR!? Thanks!

[masterleinad]

/rebuild

[dangars]

I rebased the test. I didn't have the time update the test. I'll try to find some time tomorrow to do it.

[bangerth]

@dangars Today would be the day to get it in before the release :-)

[dangars]

Hi @bangerth, I didn't have time to finish the test, I guess that this will make it in deal.II 9.4

I pushed my last version of the test, something is not working. It shouldn't be difficult to find the problem, but I need a little bit more time.

step-62 shows that find_active_cell_around_point() works well in parallel with MPI.

I'll let you know when it is ready for review!

FYI @luca-heltai

[marcfehling]

I pushed my last version of the test, something is not working. It shouldn't be difficult to find the problem, but I need a little bit more time.

Can you be more specific?

In the meantime, I will check your version against the current master branch. If it runs smoothly, I will merge the test as is.

[marcfehling]

I get errors of the following type:

37: /dealii/tests/bits/point_value_03.cc: In instantiation of ‘void check() [with int dim = 2]’:
2437: /dealii/tests/bits/point_value_03.cc:208:12:   required from here
2437: /dealii/tests/bits/point_value_03.cc:160:34: error: ‘dealii::GridTools::Cache<2, 2> cache’ has incomplete type
2437:    GridTools::Cache<dim, dim>     cache(tria, mapping);
2437:                                   ^~~~~

[marcfehling]

You missed to include a header. I will push directly to your branch and rebase on the upstream master.

The serial version of the test passes, but the parallel run with 7 processes still fails:

2439: --------------------------------------------------------
2439: An error occurred in line <108> of file <dealii/include/deal.II/numerics/vector_tools_point_value.templates.h> in function
2439:     void dealii::VectorTools::point_value(const dealii::Mapping<dim, spacedim>&, const dealii::DoFHandler<dim, spacedim>&, const VectorType&, const dealii::Point<spacedim>&, dealii::Vector<typename VectorType::value_type>&) [with int dim = 2; VectorType = dealii::PETScWrappers::MPI::Vector; int spacedim = 2; typename VectorType::value_type = double]
2439: The violated condition was: 
2439:     cell_point.first.state() == IteratorState::valid && cell_point.first->is_locally_owned()
2439: Additional information: 
2439:     The given point is inside a cell of a
2439:     parallel::distributed::Triangulation that is not locally owned by this
2439:     processor.
2439: --------------------------------------------------------

Can you have a look?

[masterleinad]

I pushed a working version.