Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project!

## Re: [Libmesh-devel] Elem::contains_point()

 Re: [Libmesh-devel] Elem::contains_point() From: Roy Stogner - 2012-06-19 17:27:15 ```On Tue, 19 Jun 2012, John Peterson wrote: > I was thinking of sending this to the list, but I wanted to be sure I > wasn't being stupid first... Copying to libmesh-devel. If we start playing by "don't send to the list unless we're certain it's not stupid" rules, nobody's going to hear from me until 2013. Sending stuff to me for a sanity check is also not likely to be productive for a few months. > So, for 2D elements living in 3D, it looks like Elem::contains_point() > maybe returns true if the point is in the projection of the element > onto the x-y plane? > In the following code snippet, the mesh contains one triangle, which > is just the "diagonal" face of the reference tet. > > If we check the point (0,0,1/2) which is directly below the face, but > not contained in it, contains_point() returns true. > > The deal seems to be that FEAbstract::on_reference_element() only > looks at xi and eta for the TRI3, which are indeed on the reference > elementn You're correct about both the problem and the cause. > How do we catch the out-of-plane cases properly? It looks like > FE::inverse_map() solves the normal equations for the "2D point living > in 3D" case, so that it returns zero in the third component of its > mapped point? That is what inverse_map does, and similarly for the 1D-in-2D/3D case. Not sure what the best solution is. Add a sanity check by calculating xyz for the results of the inverse_map and making sure the distance from the original point is under tol*h_max? Seems unnecessarily expensive in 90% of case, but I suppose we'd only need to do the check in 1D and 2D cases. --- Roy > #include "libmesh.h" > #include "mesh.h" > #include "elem.h" > > int main (int argc, char** argv) > { > { > LibMeshInit init(argc, argv); > > // Location of libmesh install > std::string libmesh_dir = std::getenv("LIBMESH_DIR"); > > // Create the Mesh > Mesh mesh(2); > > // Read in the reference element > mesh.read(libmesh_dir+std::string("/reference_elements/2D/one_tri.xda")); > > // A pointer to the element > Elem* the_elem = mesh.elem(0); > > // Move around the nodes: > the_elem->point(0) = Point(1., 0., 0.); > the_elem->point(1) = Point(0., 1., 0.); > the_elem->point(2) = Point(0., 0., 1.); > > // A point which is *not* in the triangle, but contains_point() thinks is > Point p1(0., 0., 0.5); > > // A point which is in the triangle > Point p2(.5, 0., .5); > > if ( the_elem->contains_point(p1) ) > std::cout << "Oops, contains point p1" << std::endl; > > if ( the_elem->contains_point(p2) ) > std::cout << "OK, Contains point p2" << std::endl; > > mesh.write("contains_point_test.e"); > > } // LibMeshInit > > > return 0; > } > > ```

 Re: [Libmesh-devel] Elem::contains_point() From: Roy Stogner - 2012-06-19 17:27:15 ```On Tue, 19 Jun 2012, John Peterson wrote: > I was thinking of sending this to the list, but I wanted to be sure I > wasn't being stupid first... Copying to libmesh-devel. If we start playing by "don't send to the list unless we're certain it's not stupid" rules, nobody's going to hear from me until 2013. Sending stuff to me for a sanity check is also not likely to be productive for a few months. > So, for 2D elements living in 3D, it looks like Elem::contains_point() > maybe returns true if the point is in the projection of the element > onto the x-y plane? > In the following code snippet, the mesh contains one triangle, which > is just the "diagonal" face of the reference tet. > > If we check the point (0,0,1/2) which is directly below the face, but > not contained in it, contains_point() returns true. > > The deal seems to be that FEAbstract::on_reference_element() only > looks at xi and eta for the TRI3, which are indeed on the reference > elementn You're correct about both the problem and the cause. > How do we catch the out-of-plane cases properly? It looks like > FE::inverse_map() solves the normal equations for the "2D point living > in 3D" case, so that it returns zero in the third component of its > mapped point? That is what inverse_map does, and similarly for the 1D-in-2D/3D case. Not sure what the best solution is. Add a sanity check by calculating xyz for the results of the inverse_map and making sure the distance from the original point is under tol*h_max? Seems unnecessarily expensive in 90% of case, but I suppose we'd only need to do the check in 1D and 2D cases. --- Roy > #include "libmesh.h" > #include "mesh.h" > #include "elem.h" > > int main (int argc, char** argv) > { > { > LibMeshInit init(argc, argv); > > // Location of libmesh install > std::string libmesh_dir = std::getenv("LIBMESH_DIR"); > > // Create the Mesh > Mesh mesh(2); > > // Read in the reference element > mesh.read(libmesh_dir+std::string("/reference_elements/2D/one_tri.xda")); > > // A pointer to the element > Elem* the_elem = mesh.elem(0); > > // Move around the nodes: > the_elem->point(0) = Point(1., 0., 0.); > the_elem->point(1) = Point(0., 1., 0.); > the_elem->point(2) = Point(0., 0., 1.); > > // A point which is *not* in the triangle, but contains_point() thinks is > Point p1(0., 0., 0.5); > > // A point which is in the triangle > Point p2(.5, 0., .5); > > if ( the_elem->contains_point(p1) ) > std::cout << "Oops, contains point p1" << std::endl; > > if ( the_elem->contains_point(p2) ) > std::cout << "OK, Contains point p2" << std::endl; > > mesh.write("contains_point_test.e"); > > } // LibMeshInit > > > return 0; > } > > ```
 Re: [Libmesh-devel] Elem::contains_point() From: John Peterson - 2012-06-19 19:15:13 Attachments: 0001-Fixing-Elem-point_test-called-by-Elem-contains_point.patch ```On Tue, Jun 19, 2012 at 11:27 AM, Roy Stogner wrote: > > Not sure what the best solution is.  Add a sanity check by calculating > xyz for the results of the inverse_map and making sure the distance > from the original point is under tol*h_max?  Seems unnecessarily > expensive in 90% of case, but I suppose we'd only need to do the check > in 1D and 2D cases. The attached patch implements this fix, it works for my simple test case and a more involved one I have... It could be made a bit faster by adding an FEInterface::map() function in the same vein as FEInterface::inverse_map() which forwards to FE::map(). This seems like it would also be a useful feature in general... An even better optimization would be to add geometric-only specializations of Elem::contains_point(), at least for the linear element types. -- John ```
 Re: [Libmesh-devel] Elem::contains_point() From: John Peterson - 2012-06-19 19:58:04 ```On Tue, Jun 19, 2012 at 1:14 PM, John Peterson wrote: > On Tue, Jun 19, 2012 at 11:27 AM, Roy Stogner wrote: >> >> Not sure what the best solution is.  Add a sanity check by calculating >> xyz for the results of the inverse_map and making sure the distance >> from the original point is under tol*h_max?  Seems unnecessarily >> expensive in 90% of case, but I suppose we'd only need to do the check >> in 1D and 2D cases. > > The attached patch implements this fix, it works for my simple test > case and a more involved one I have... > > It could be made a bit faster by adding an FEInterface::map() function > in the same vein as FEInterface::inverse_map() which forwards to > FE::map().  This seems like it would also be a useful feature > in general... OK, FEInterface::map() is added and Elem::point_test() uses it. Considering that Elem::point_test() already contains a call to FEInterface::inverse_map(), adding one more call to FE::map() makes it only trivially more expensive... -- John ```
 Re: [Libmesh-devel] Elem::contains_point() From: Roy Stogner - 2012-06-19 20:29:07 ``` On Tue, 19 Jun 2012, John Peterson wrote: > On Tue, Jun 19, 2012 at 1:14 PM, John Peterson wrote: >> On Tue, Jun 19, 2012 at 11:27 AM, Roy Stogner wrote: >>> >>> Not sure what the best solution is.  Add a sanity check by calculating >>> xyz for the results of the inverse_map and making sure the distance >>> from the original point is under tol*h_max?  Seems unnecessarily >>> expensive in 90% of case, but I suppose we'd only need to do the check >>> in 1D and 2D cases. >> >> The attached patch implements this fix, it works for my simple test >> case and a more involved one I have... >> >> It could be made a bit faster by adding an FEInterface::map() function >> in the same vein as FEInterface::inverse_map() which forwards to >> FE::map().  This seems like it would also be a useful feature >> in general... > > OK, FEInterface::map() is added and Elem::point_test() uses it. This looks good; thanks! > Considering that Elem::point_test() already contains a call to > FEInterface::inverse_map(), adding one more call to FE::map() makes it > only trivially more expensive... You're probably right. --- Roy```