Fix: solving issues of find active cell around point(cache) with points at the edge #5512

GivAlz · 2017-11-21T14:03:17Z

Dear everyone,
working on #5411 I've found a problem with the Cached version of compute point locations: the new version is failing, there's probably some issue when looking for points which are on the edge of the Triangulation.

This test currently fails: it shall pass once the bug is found and corrected ( @luca-heltai I guess that's also something I shall need to do).

Best,
Giovanni

GivAlz · 2017-11-23T21:01:53Z

Update: I found the problem, it was with tolerance: points at the edge where sometimes trashed.
(As soon as I update I noticed a bug in my "patch", so here's the comment to the new version):

What I did was adapting what happens in the old version of find active cell around point: every time a cell is checked, if the point is not inside we check the distance of the transformed point to the unit cell. If it is less than the tolerance, we keep track of this "good" cell and lower the tolerance... if, in the end, no "inside" cell was found but we found a good cell, this is returned.

I really think this shall make (finally) #5411 work :)

bangerth

This is a change in the behavior of this function. I think it should be documented.

bangerth · 2017-11-24T03:14:57Z

source/grid/grid_tools.cc

+        // It is possible the vertex is close
+        // to an edge, thus we add a tolerance
+        // setting it initially to 1e-10
+        // to keep also the "best" point


"best" point -> "best" cell

bangerth · 2017-11-24T03:15:01Z

source/grid/grid_tools.cc

+                    best_distance = dist;
+                    cell_and_position_approx.first = *cell;
+                    cell_and_position_approx.second = p_unit;
+                    approx_cell = true;


I think you lack a break statement here. Or is your intention to keep searching through the other cells? If so, I think a comment in front of this block would be nice of the sort "The point is not inside this cell, but let us see how far outside it is and whether we want to use this cell as a backup if we can't find a cell within which it is."

I'd say that it would also be easier to read if in the block above (the one where p_unit is inside the cell) you'd set approx_cell = false, just to make clear that you are discarding the approximate cell.

My intention is to keep searching: maybe the point is close to the edge but inside another cell. Thus I want to return it only if I'm sure it's inside a cell or it's the best approximation among all possible cases, to avoid "false positives".

I believe this is also the thinking behind
find_active_cell_around_point (const Mapping<dim,spacedim> &mapping, const MeshType<dim,spacedim> &mesh, const Point<spacedim> &p, const std::vector<bool> &marked_vertices)
but, for how that function is made, I believe currently best_distance is not really used (because, as soon as a "close enough point" is found we exit the cycle by setting found to true... so it doesn't make sense to set best_distance = dist; but this is for another PR).

Thank you for your comments: updating soon

GivAlz · 2017-11-24T09:10:00Z

@bangerth for the change in behaviour which should be documented: actually the documentation of all find_active_cell_around_point is of the like "a version of the previous function with...". The original one is the following function:

find_active_cell_around_point (const Mapping<dim,spacedim> &mapping, const MeshType<dim,spacedim> &mesh, const Point<spacedim> &p, const std::vector<bool> &marked_vertices = std::vector<bool>());
which has a note explaining how, for points on the edge, we can only approximate things. So I believe the documentation was already there, only new versions didn't take into account "problems at the edge"
... I can always put a remark on the other versions of the function, if you think it's better!

Also: the way I shall return the cell it's slightly different from this original find_active_cell_around_point version where, as soon as the point is close enough, it is returned...
I personally prefer to check all cells until we finish them or we found the point inside, so I'd rather change the old find_active_cell_around_point function (maybe in the future), so that the outputs are as similar as possible, but this is clearly only my opinion.

bangerth · 2017-11-25T05:43:48Z

Thanks for digging into the documentation -- if it's already mentioned elsewhere, then that's good enough. No need for further changes!

/run-tests
OK to merge once the tester is happy.

GivAlz · 2017-11-25T08:42:34Z

Ok, grid/find_active_cell_around_point_03.debug is failing which is bad (but the release version is passing, at least on my machine); I'm looking into this problem right now

GivAlz · 2017-11-25T09:25:28Z

Ok, problem fixed: the debug failed because I used something like
auto my_pair_1 = find_active_cell_around_point
with the old version of the function which has two versions, differing only in output, and the "auto" worked differently in the debug and release version (also on my machine).
I also removed some unnecessary headers from the test and changed the assert to a gentler if(cells are different) deallog...

I think I shall update this in the documentation and fix a few things in the old find active cell around point (just not right now: my schedule it's a bit too busy for the next week).

… and test

bangerth · 2017-11-25T15:42:22Z

On 11/25/2017 02:25 AM, Giovanni Alzetta wrote: Ok, problem fixed: the debug failed because I used something like |auto my_pair_1 = find_active_cell_around_point| with the old version of the function which has two versions, differing only in output, and the "auto" worked differently in the debug and release version (also on my machine).

I don't think I understand this. Can you elaborate?

GivAlz · 2017-11-25T16:22:04Z

@bangerth it's actually easy, so sorry if I managed not to be clear immediately!
There are two functions (find_active_cell_around_point) which have the same input parameters, but different outputs (in one case a pair, in the other just the second element of the pair).

In the first version of the test, the one which failed, instead of writing
std::pair<...> output = find_active_cell (... )
I wrote
auto output = find_active_cell (...)
Thus what happened is that the compiler decided which function to use and, apparently, when compiling for the debug version it chose the one I didn't want to use while for the release version it chose the version I wanted to use.
I believe this possible problem should be pointed out in the documentation...and I'd also change a couple of things in the "old" find active cell function (when I shall have time to do a new PR about this).
Thanks,
Giovanni

bangerth · 2017-11-26T23:05:26Z

On 11/25/2017 09:22 AM, Giovanni Alzetta wrote:

@bangerth https://github.com/bangerth it's actually easy. There are two functions (find_active_cell_around_point) which have the same input parameters, but different outputs (in one case a pair, in the other just the second element).

No, that can't be it. C++ does not allow this. Functions need to differ in input argument types or a compiler can not determine which function to call. (C++ does not consider the return type as part of overload resolution.)

In the first version of the test, instead of writing
|std::pair<...> output = find_active_cell (... )|
I wrote
|auto output = find_active_cell (...)|
Thus what happened is that the compiler decided which function to use and, apparently, when compiling for the debug version it chose the one I didn't want to use.

That also can't be it -- we do not let the compiler choose a different function call in debug than in release mode. There must be another reason for your issue.

GivAlz · 2017-11-27T10:49:18Z

@bangerth then I really don't know where the problem was and, I guess, it was just fixed by luck: if you take line 61-62 of the grid/find_active_cell_around_point_03 from
std::pair< typename Triangulation<dim,spacedim>::active_cell_iterator, Point< dim > > res_1 = GridTools::find_active_cell_around_point (mapping, dof_handler, p1);

to auto = GridTools::find_active_cell_around_point(mapping, dof_handler, p1); the test fails and you get the following error:

2054: In file included from /home/my_nick/opt/src/dealii_fork/include/deal.II/grid/tria_accessor.templates.h:26:0,
2054:                  from /home/my_nick/opt/src/dealii_fork/include/deal.II/grid/tria_accessor.h:3409,
2054:                  from /home/my_nick/opt/src/dealii_fork/include/deal.II/grid/tria.h:3898,
2054:                  from /home/my_nick/opt/src/dealii_fork/include/deal.II/fe/mapping.h:22,
2054:                  from /home/my_nick/opt/src/dealii_fork/include/deal.II/fe/fe.h:25,
2054:                  from /home/my_nick/opt/src/dealii_fork/include/deal.II/hp/fe_collection.h:20,
2054:                  from /home/my_nick/opt/src/dealii_fork/include/deal.II/dofs/dof_handler.h:32,
2054:                  from /home/my_nick/opt/src/dealii_fork/tests/grid/find_active_cell_around_point_03.cc:23:
2054: /home/my_nick/opt/src/dealii_fork/include/deal.II/grid/tria_iterator.h: In instantiation of ‘dealii::TriaRawIterator<Accessor>::TriaRawIterator(const dealii::TriaActiveIterator<OtherAccessor>&) [with OtherAccessor = dealii::CellAccessor<1, 1>; Accessor = dealii::DoFCellAccessor<dealii::DoFHandler<1, 1>, false>]’:
2054: /home/giovy/opt/src/dealii_fork/tests/grid/find_active_cell_around_point_03.cc:82:20:   required from ‘void test() [with int dim = 1; int spacedim = 1]’
2054: /home/my_nick/opt/src/dealii_fork/tests/grid/find_active_cell_around_point_03.cc:95:14:   required from here
2054: /home/my_nick/opt/src/dealii_fork/include/deal.II/grid/tria_iterator.h:904:23: error: no matching function for call to ‘dealii::DoFCellAccessor<dealii::DoFHandler<1, 1>, false>::DoFCellAccessor(const dealii::CellAccessor<1, 1>&)’
2054:    accessor (i.accessor)

My guess was confusion between the following functions (which, actually, do not have the exact same input arguments, but in the second one the last variable has a default value); but honestly I don't know C++ enough and so, if you say this simply can't be...it means I really don't know why the change I made makes everything work...

  template <int dim, template <int, int> class MeshType, int spacedim>
#ifndef _MSC_VER
  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
#else
  std::pair<typename dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
#endif
  find_active_cell_around_point (const Mapping<dim,spacedim>  &mapping,
                                 const MeshType<dim,spacedim> &mesh,
                                 const Point<spacedim>        &p,
                                 const std::vector<bool>      &marked_vertices = std::vector<bool>());

and

  template <int dim, int spacedim>
  std::pair<typename hp::DoFHandler<dim, spacedim>::active_cell_iterator, Point<dim> >
  find_active_cell_around_point (const hp::MappingCollection<dim,spacedim> &mapping,
                                 const hp::DoFHandler<dim,spacedim>        &mesh,
                                 const Point<spacedim>                     &p);

masterleinad · 2017-11-27T16:38:32Z

2054: /home/my_nick/opt/src/dealii_fork/include/deal.II/grid/tria_iterator.h:904:23: error: no matching function for call to ‘dealii::DoFCellAccessor<dealii::DoFHandler<1, 1>, false>::DoFCellAccessor(const dealii::CellAccessor<1, 1>&)’

The problem is simply that you can't convert a CellAccesor to a (richer) DoFCellAccessor while the converse works. The return type of the non-cached version is of DoFCellAccessor type while your version only provides a CellAccessor.
Just swapping the arguments in your comparisons would also work with using auto for all the variables in question.

masterleinad · 2017-11-27T16:41:53Z

Since, @bangerth didn't request any further changes, this looks good to me and all the test pass, I merge.

GivAlz · 2017-11-27T16:50:00Z

Thank you @masterleinad , I really coudln't figure out why it was now working!!

islent · 2018-09-05T13:26:11Z

The problem is simply that you can't convert a CellAccesor to a (richer) DoFCellAccessor while the converse works. The return type of the non-cached version is of DoFCellAccessor type while your version only provides a CellAccessor.
Just swapping the arguments in your comparisons would also work with using auto for all the variables in question.

@masterleinad Sorry for bothering you, but I seem to be facing the same problem as above.

I'm using find_active_cell_around_point as following

GridTools::Cache<dim, spacedim> cache(dof.get_triangulation(), mapping);

const std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator, Point<spacedim>>
    cell_point = GridTools::find_active_cell_around_point(cache, point);`

I've tried tons of different methods, only to get the errors again and again

/usr/local/include/deal.II/grid/tria_iterator.h:938:15: error: no matching function for call to 
‘dealii::DoFCellAccessor<dealii::DoFHandler<3, 3>, false>::DoFCellAccessor(const 
dealii::CellAccessor<3, 3>&)’ : accessor(a)

I couldn't figure out where the program convert the CellAccesor to DoFCellAccessor when I simply call a cached version of find_active_cell_around_point

Could you please to explain

Just swapping the arguments in your comparisons would also work with using auto for all the variables in question.

in more details?

And is there any examples about how to use a cached version of find_active_cell_around_point or how to use the kdtree version of it (One of your developer said it did have a kdtree version)?

Thanks a lot :)

jppelteret · 2018-09-05T13:29:29Z

I couldn't figure out where the program convert the CellAccesor to DoFCellAccessor

There's an entry in the FAQ on this: Can I convert Triangulation cell iterators to DoFHandler cell iterators?

islent · 2018-09-05T13:49:43Z

@jppelteret Thanks a lot, but I'm still confused of when the conversion happens and how to avoid it.

Do you mean that I should do the conversion manually or something?

drwells · 2018-09-05T14:00:46Z

Would it be possible to move this discussion either to the slack channel or the mailing list? We should try to solve this, but this GitHub issue is not the best location to do so.

jppelteret · 2018-09-05T14:01:01Z

Do you mean that I should do the conversion manually or something?

Yes, I guess so. If you have a look at the signature of the specific function you are calling, returns a std::pair< typename Triangulation< dim, spacedim >::active_cell_iterator, Point< dim > >. This Tria::active_cell_iterator cannot be automatically cast "up" to a DofHandler::active_cell_iterator, so you must convert it yourself.

It should be possible to implement another more general version of this function, say

template<int dim, template< int, int > class MeshType, int spacedim>
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > 
  GridTools::find_active_cell_around_point  ( 
    const Cache< dim, spacedim > & cache,
    const MeshType< dim, spacedim > & mesh,
    const Point< spacedim > &  p,
    const typename MeshType< dim, spacedim >::active_cell_iterator & cell_hint = typename MeshType<dim, spacedim>::active_cell_iterator(), 
    const std::vector< bool > & marked_vertices = std::vector<bool>() )

that does the conversion automatically. If you would like to implement this feature then we'd be happy to include it in the next release.

islent · 2018-09-05T14:35:17Z

@drwells Sorry for that, I'm a newer both to dealii and github, and didn't realized that I might have bothered a lot of people. I've been trapped by this error for almost a month and when came to this related issue I just couldn't help putting forward my questions.

And thank you @jppelteret for your patience and kindness, I would try this out and attempt to implement the automatic conversion later on. It would be a great honor to have any contributions to the dealii project!

drwells · 2018-09-05T14:47:32Z

Don't worry about it. Its more important that your questions get answered (but it does help if things stay sorted).

davydden added the Triangulation label Nov 21, 2017

GivAlz force-pushed the Failing_find_active branch from 3e7d4f2 to cdc85a3 Compare November 23, 2017 20:34

GivAlz changed the title ~~Added test for corner cases find active cell around point~~ Fixing: adding tolerange and test for corner cases of find active cell around point Nov 23, 2017

GivAlz force-pushed the Failing_find_active branch 2 times, most recently from 5eac006 to 164591f Compare November 23, 2017 20:58

GivAlz changed the title ~~Fixing: adding tolerange and test for corner cases of find active cell around point~~ Fixing: adding tolerance and test for corner cases of find active cell around point Nov 23, 2017

GivAlz changed the title ~~Fixing: adding tolerance and test for corner cases of find active cell around point~~ Fix: solving issues of find active cell around point(cache) with points at the edge Nov 23, 2017

bangerth reviewed Nov 24, 2017

View reviewed changes

bangerth added the Reviewed and ready to merge label Nov 25, 2017

GivAlz force-pushed the Failing_find_active branch from 8d8e67b to b313818 Compare November 25, 2017 09:20

Added tolerance and test for GridTools::find_active_cell_around_point…

b313818

… and test

bangerth mentioned this pull request Nov 25, 2017

Make fe_field_function compute::point locations use the gridtools version #5411

Merged

masterleinad merged commit 18cda4d into dealii:master Nov 27, 2017

islent mentioned this pull request Sep 14, 2018

Cached find_active_cell_around_point() incompatible with cell iterator of Triangulation in point_gradient() #7188

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix: solving issues of find active cell around point(cache) with points at the edge #5512

Fix: solving issues of find active cell around point(cache) with points at the edge #5512

GivAlz commented Nov 21, 2017

GivAlz commented Nov 23, 2017 •

edited

bangerth left a comment

bangerth Nov 24, 2017

bangerth Nov 24, 2017

GivAlz Nov 24, 2017 •

edited

GivAlz commented Nov 24, 2017 •

edited

bangerth commented Nov 25, 2017

GivAlz commented Nov 25, 2017

GivAlz commented Nov 25, 2017

bangerth commented Nov 25, 2017 via email

GivAlz commented Nov 25, 2017 •

edited

bangerth commented Nov 26, 2017

GivAlz commented Nov 27, 2017 •

edited

masterleinad commented Nov 27, 2017

masterleinad commented Nov 27, 2017

GivAlz commented Nov 27, 2017

islent commented Sep 5, 2018

jppelteret commented Sep 5, 2018

islent commented Sep 5, 2018

drwells commented Sep 5, 2018

jppelteret commented Sep 5, 2018

islent commented Sep 5, 2018

drwells commented Sep 5, 2018

Fix: solving issues of find active cell around point(cache) with points at the edge #5512

Fix: solving issues of find active cell around point(cache) with points at the edge #5512

Conversation

GivAlz commented Nov 21, 2017

GivAlz commented Nov 23, 2017 • edited

bangerth left a comment

Choose a reason for hiding this comment

bangerth Nov 24, 2017

Choose a reason for hiding this comment

bangerth Nov 24, 2017

Choose a reason for hiding this comment

GivAlz Nov 24, 2017 • edited

Choose a reason for hiding this comment

GivAlz commented Nov 24, 2017 • edited

bangerth commented Nov 25, 2017

GivAlz commented Nov 25, 2017

GivAlz commented Nov 25, 2017

bangerth commented Nov 25, 2017 via email

GivAlz commented Nov 25, 2017 • edited

bangerth commented Nov 26, 2017

GivAlz commented Nov 27, 2017 • edited

masterleinad commented Nov 27, 2017

masterleinad commented Nov 27, 2017

GivAlz commented Nov 27, 2017

islent commented Sep 5, 2018

jppelteret commented Sep 5, 2018

islent commented Sep 5, 2018

drwells commented Sep 5, 2018

jppelteret commented Sep 5, 2018

islent commented Sep 5, 2018

drwells commented Sep 5, 2018

GivAlz commented Nov 23, 2017 •

edited

GivAlz Nov 24, 2017 •

edited

GivAlz commented Nov 24, 2017 •

edited

GivAlz commented Nov 25, 2017 •

edited

GivAlz commented Nov 27, 2017 •

edited