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;
> }
>
>
------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
