Usage of FEFieldFunction.vector_value_list on a parallel::distributed::Triangulation

[Maxi Miller]

I would like to use the function vector_value_list for the interpolation between two different grids (one externally, one internally) instead of having to loop over all points and calling value_list for each point. Now I already know that not all points are within the range of the current MPI-thread, thus both functions will throw ExcPointNotAvailableHere sooner or later. If that happens, I can deal with it when looping over every single point (by simply hopping over it). If that happens when calling vector_value_list, how does the function behave then? Does it throw the error, but all points which can be set have been set, or does it throw the error before all possible points have been set? If the latter, how can I make sure that I got all possible points in that thread?

[Wolfgang Bangerth]

I had to look at the function's implementation and it turns out that it throws
the exception the first time it encounters a cell that is not locally owned.
It would be a reasonable idea to restructure the function in such a way that
it first fills the data for all points for which this is possible and only
then throws the exception. Would that be useful to you? Would you be
interested in writing a patch to this end?

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Maxi Miller]

Hei,

that would definitely be useful for me (else that function is not usable when running multiple MPI threads). I can take a look at the function itself, but am not sure if I can write a patch for it without additional help.

[Wolfgang Bangerth]

On 1/2/19 3:24 AM, 'Maxi Miller' via deal.II User Group wrote:
> that would definitely be useful for me (else that function is not usable when
> running multiple MPI threads). I can take a look at the function itself, but
> am not sure if I can write a patch for it without additional help.

That's where we're happy to help! Take a look at the implementation of the
function here:

https://github.com/dealii/dealii/blob/master/include/deal.II/numerics/fe_field_function.templates.h#L263

The exception is thrown in line 297/298. In essence, what one would need to do
is this:

(i) instead of throwing an exception, just condition the code in lines 300-308
with an `if (cell->is_artificial() == false)` and just ignore all points for
which the cell is artificial;

(ii) add an `else` clause in case where the cell is artificial in which we
fill the output `values` array with invalid values;

(iii) *after* the loop, if one of the cells was artificial, throw the exception.

Want to give this a try? A guide on how to contribute patches is here:
https://github.com/dealii/dealii/wiki/Contributing

[Maxi Miller]

I will take a look, but it could take a while (week, or so?), due to a lot of other things I also have to do...

[Wolfgang Bangerth]

On 1/5/19 4:52 AM, 'Maxi Miller' via deal.II User Group wrote:
> I will take a look, but it could take a while (week, or so?), due to a lot of
> other things I also have to do...

Take your time -- any contribution is welcome!

Best
Wolfgang

[Maxi Miller]

I implemented the suggestion, but get an error at

        const unsigned int n_cells =
                compute_point_locations(points, cells, qpoints, maps);

with the error

An error occurred in line <1816> of file </home/roland/Downloads/git-files/dealii/source/grid/grid_tools.cc> in function
    std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, dealii::Point<dim> > dealii::GridTools::find_active_cell_around_point(const dealii::Mapping<dim, spacedim>&, const MeshType<dim, spacedim>&, const dealii::Point<spacedim>&, const std::vector<std::set<typename MeshType<dim, spacedim>::active_cell_iterator> >&, const std::vector<std::vector<dealii::Tensor<1, spacedim> > >&, const typename MeshType<dim, spacedim>::active_cell_iterator&, const std::vector<bool>&, dealii::RTree<std::pair<dealii::Point<spacedim>, unsigned int> >&) [with int dim = 2; MeshType = dealii::Triangulation; int spacedim = 2; typename MeshType<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >; typename MeshType<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >; typename MeshType<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >; typename MeshType<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::CellAccessor<2, 2> >; dealii::RTree<std::pair<dealii::Point<spacedim>, unsigned int> > = boost::geometry::index::rtree<std::pair<dealii::Point<2, double>, unsigned int>, boost::geometry::index::linear<16>, boost::geometry::index::indexable<std::pair<dealii::Point<2, double>, unsigned int> >, boost::geometry::index::equal_to<std::pair<dealii::Point<2, double>, unsigned int> >, std::allocator<std::pair<dealii::Point<2, double>, unsigned int> > >]
The violated condition was:
    current_cell.state() == IteratorState::valid
Additional information:
    The point <0.5 0.8> could not be found inside any of the subcells of a coarse grid cell.

I assume it tries to find the point, which is not available due to being within a cell located in another thread. I assume I have to modify the cell iterator in that case. Is that correct?

[giovann...@hotmail.it]

Hi everyone!

  Sorry for coming up this late, but I just discovered  this behaviour is important also for some functions I am working on (and also to why distributed compute point locations did not work) and causes this:

https://github.com/dealii/dealii/issues/7532 .

What's happening is that in the past find active cell around point would would also return artificial cells. Now something changed: in the past the functions related to non_matching, and distributed point locations, relied on the fact that if the cell was anywhere on the mesh it would be returned. Then it could be checked by hand

if the returned cell was "of the coarse grid cell" or not.

I shall open an issue about this, as at least in some cases it would be desiderable for the "point not found" exception to be ignored.

[Wolfgang Bangerth]

On 1/8/19 12:43 AM, 'Maxi Miller' via deal.II User Group wrote:
>
> I assume it tries to find the point, which is not available due to being
> within a cell located in another thread. I assume I have to modify the cell
> iterator in that case. Is that correct?

You mean owned by some other MPI process? I don't know whether you
interpretation is correct. The exception is found here:

// The first time around, we check for vertices in the hint_cell. If
// that does not work, we set the cell iterator to an invalid one, and
// look for a global vertex close to the point. If that does not work,
// we are in trouble, and just throw an exception.
//
// If we got here, then we did not find the point. If the
// current_cell.state() here is not IteratorState::valid, it means that
// the user did not provide a hint_cell, and at the beginning of the
// while loop we performed an actual global search on the mesh
// vertices. Not finding the point then means the point is outside the
// domain.
AssertThrow(current_cell.state() == IteratorState::valid,
ExcPointNotFound<spacedim>(p));

Do you think you can come up with a small testcase?
Best
W.