It sounds like this is an inherent problem of the present algorithm for determining if a point is in an element??? especially for complex elements such as HEX27. Is there any other more robust solution in computational geometry?
I also tried to use the 1st order element, it passed the test for this simple case, although I am not sure if this really solves the problem. Is it possible to add another contains_point(Point& pt, ORDER order) that allows using lower order shape function? In CFD with mixed finite elements, we sometimes use higher order interpolation for velocity and lower order ones for pressure. But when to determine which element a point moves into, we can still use first order algorithm as an approximation. On Mon, Jul 18, 2016 at 1:26 PM, John Peterson <jwpeter...@gmail.com> wrote: > > > On Mon, Jul 18, 2016 at 12:13 PM, Roy Stogner <royst...@ices.utexas.edu> > wrote: > >> >> On Mon, 18 Jul 2016, Xujun Zhao wrote: >> >> Yes, the points are not quad-points but arbitrary points, which are moving >>> particles in the domain. It mainly happened when the points are out of >>> the >>> element. >>> >> >> Oh, yes, looking for points far outside an element can give you >> inverse map failures. If you imagine extending the master->physical >> element map for, say, a QUAD4 trapezoid with one very narrow side on >> top, you can see how the map stops being invertible a little way above >> that top. We avoid ever hitting many such cases because we precede >> the inverse_map with a bounding box test, but that's a performance >> optimization, not a guarantee of an invertible map. We pass >> secure=false to inverse_map to, in theory, avoid screaming and dying >> when handling a case that might be expected to fail. >> >> Looks like there's a bug in the 3D case, though? We rely on >> TypeTensor::solve() in the Newton iteration in 3D, but we *don't* warn >> solve() (or have a way to warn solve()?) that in the secure=false case >> we don't necessarily expect a solution to exist. >> > > I can fix this. The zero determinant asserts could be changed to > libmesh_error()'s (or explicit exception throws) and then we could catch > the exception in inverse_map(). It should be an "exceptional" case, so it > seems reasonable to me. > > Other options would be to return an error code, or pass an extra flag to > TypeTensor::solve()? In these cases you'd have to check the flag every > time... > > -- > John > ------------------------------------------------------------------------------ What NetFlow Analyzer can do for you? Monitors network bandwidth and traffic patterns at an interface-level. Reveals which users, apps, and protocols are consuming the most bandwidth. Provides multi-vendor support for NetFlow, J-Flow, sFlow and other flows. Make informed decisions using capacity planning reports.http://sdm.link/zohodev2dev _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
