## Re: [Libmesh-users] Wrong result of bool Elem::contains_point (const Point & p )

 Re: [Libmesh-users] Wrong result of bool Elem::contains_point (const Point & p ) From: Roy Stogner - 2005-12-01 17:43:24 ```On Thu, 1 Dec 2005, Michael Povolotskyi wrote: > The point does not belong to the element. However, the method returns 1. Would you see what happens with the inverse_map of that point? Try tossing a "std::cout << mapped_point;" in the Elem::contains_point() routine and see where the inverse map thinks the point would lie on the master element. Better still, change the "false" argument at the end of inverse_map() to "true", and let's see if the Newton iteration reports divergence. My guess as to what's happening is that your physical point is so far from your element that the inverse map onto the master element isn't unique, so the Newton's method to find a master point is failing and quitting on an iterate that happens to fall inside the element. If that's the problem, I'm not sure what the fix would be. We could add "sanity checks" before calling inverse_map() to make sure the point is inside the element's bounding box, and that would catch this particular case, but it's easy to come up with cases that wouldn't be fixed. --- Roy ```

 [Libmesh-users] Wrong result of bool Elem::contains_point (const Point & p ) From: Michael Povolotskyi - 2005-12-01 16:39:19 ```Dear Libmesh developers, I found a wrong result produced by method bool Elem::contains_point (const Point & p ). The case is as follows. The element is of Lagrange type. It has 8 vertexes with such coordinates: x y z 11.8263 3.79213 0 -0.536113 1.015 0 -3.15014 -11.2129 0 14.5533 -13.1308 0 11.8263 3.79213 2.41439 -0.536113 1.015 2.41439 -3.15014 -11.2129 2.41439 14.5533 -13.1308 2.41439 The point I want to check has such coordinates: -28.266 23.1273 0; The point does not belong to the element. However, the method returns 1. Thank you, Michael. ```
 Re: [Libmesh-users] Wrong result of bool Elem::contains_point (const Point & p ) From: Roy Stogner - 2005-12-01 17:43:24 ```On Thu, 1 Dec 2005, Michael Povolotskyi wrote: > The point does not belong to the element. However, the method returns 1. Would you see what happens with the inverse_map of that point? Try tossing a "std::cout << mapped_point;" in the Elem::contains_point() routine and see where the inverse map thinks the point would lie on the master element. Better still, change the "false" argument at the end of inverse_map() to "true", and let's see if the Newton iteration reports divergence. My guess as to what's happening is that your physical point is so far from your element that the inverse map onto the master element isn't unique, so the Newton's method to find a master point is failing and quitting on an iterate that happens to fall inside the element. If that's the problem, I'm not sure what the fix would be. We could add "sanity checks" before calling inverse_map() to make sure the point is inside the element's bounding box, and that would catch this particular case, but it's easy to come up with cases that wouldn't be fixed. --- Roy ```
 Re: [Libmesh-users] Wrong result of bool Elem::contains_point (const Point & p ) From: Michael Povolotskyi - 2005-12-01 18:32:33 ```Hello Roy, unfortunatly now I don't have enough time to recompile the library. But when I called relative_point = FEInterface::inverse_map(dim, fe_type, elem, point); I've got something like this: WARNING: Newton scheme has not converged in 11 iterations: physical_point=(x,y,z)=( -28.266, 23.1273, 0) physical_guess=(x,y,z)=(-2.54221, 0.51462, 0) dp=(x,y,z)=(-5.63852, 5.29642, 0) p=(x,y,z)=(-0.679189, 0.889012, -1) error=7.73596 I coded suggested by you the "sanity check" before calling contains_point method and this solved my problem. I'd like to say that for Lagrange elements one can code a non-iterative algorithm in order to find if a point is inside a polyhedral or not. Also one can substitute the existing Newton method that is very sensitive to the initial guess by a more stable one that is implemented in PETSc. Michael. Roy Stogner wrote: > On Thu, 1 Dec 2005, Michael Povolotskyi wrote: > >> The point does not belong to the element. However, the method returns 1. > > > Would you see what happens with the inverse_map of that point? Try > tossing a "std::cout << mapped_point;" in the Elem::contains_point() > routine and see where the inverse map thinks the point would lie on > the master element. Better still, change the "false" argument at the > end of inverse_map() to "true", and let's see if the Newton iteration > reports divergence. > > My guess as to what's happening is that your physical point is so far > from your element that the inverse map onto the master element isn't > unique, so the Newton's method to find a master point is failing and > quitting on an iterate that happens to fall inside the element. If > that's the problem, I'm not sure what the fix would be. We could add > "sanity checks" before calling inverse_map() to make sure the point is > inside the element's bounding box, and that would catch this > particular case, but it's easy to come up with cases that wouldn't be > fixed. > --- > Roy > -- ------------------------------------------------------------ Michael Povolotskyi, Ph.D. University of Rome "Tor Vergata" Department of Electronic Engineering Viale Politecnico, 1 - 00133 Rome - Italy Phone + 39 06 72597781 Fax + 39 06 2020519 http://www.optolab.uniroma2.it/pages/moshe/moshe.html ------------------------------------------------------------- ```
 Re: [Libmesh-users] Wrong result of bool Elem::contains_point (const Point & p ) From: Roy Stogner - 2005-12-01 20:17:05 ```On Thu, 1 Dec 2005, Michael Povolotskyi wrote: > I'd like to say that for Lagrange elements one can code a non-iterative > algorithm in order to find if a point is inside a polyhedral or not. Does this apply even if the region is not polyhedral or polygonal? It's easy for me to see how to quickly test if a point is within a collection of flat sides, but with second-order elements it's possible for libMesh users to use curved sides, and with hexes and prisms it's almost impossible for users to avoid curved sides. --- Roy ```