On Mon, Jul 18, 2016 at 9:42 AM, Xujun Zhao <xzha...@gmail.com> wrote:

> Guys,
>
> I met this error when I try to determine if a point is inside an element
> using the function constains_point(). I tracked the code using the
> dbg version and found that this happened when computing the
> FE::inverse_map(), in which it calls TypeTensor::inverse() and leads to
> Jacobian = 0. Here are the error messages:
>
> Assertion `my_det != static_cast<T>(0.)' failed.
>
> my_det = 0
>
> static_cast<T>(0.) = 0
>
>
> Stack frames: 12
>
> 0: 0   libmesh_dbg.0.dylib                 0x000000010b72db6c
> libMesh::print_trace(std::ostream&) + 472
>
> 1: 1   libmesh_dbg.0.dylib                 0x000000010b729122
> libMesh::MacroFunctions::report_error(char const*, int, char const*, char
> const*) + 387
>
> 2: 2   libmesh_dbg.0.dylib                 0x000000010bde82f2
> libMesh::TypeTensor<double>::inverse() const + 666
>
> 3: 3   libmesh_dbg.0.dylib                 0x000000010b93c35a
> libMesh::FE<3u, (libMesh::FEFamily)0>::inverse_map(libMesh::Elem const*,
> libMesh::Point const&, double, bool) + 1250
>
> 4: 4   libmesh_dbg.0.dylib                 0x000000010b8a101d
> libMesh::FEInterface::inverse_map(unsigned int, libMesh::FEType const&,
> libMesh::Elem const*, libMesh::Point const&, double, bool) + 3725
>
> 5: 5   libmesh_dbg.0.dylib                 0x000000010bae6f18
> libMesh::Elem::point_test(libMesh::Point const&, double, double) const +
> 1744
>
> 6: 6   libmesh_dbg.0.dylib                 0x000000010bae680a
> libMesh::Elem::contains_point(libMesh::Point const&, double) const + 250
>
> 7: 7   example-dbg                         0x000000010b2e7fd2
> libMesh::PointMesh<3u>::build_elem_neighbor_list() + 1158
>
> 8: 8   example-dbg                         0x000000010b2e8aae
> libMesh::PointMesh<3u>::reinit() + 534
>
> 9: 9   example-dbg                         0x000000010b275e07
> test_cell(libMesh::Parallel::Communicator const&) + 4146
>
>
> I am curious why this happens. It seems to me that J can be zero when the
> volume of an element is zero. But I let elem->volume() to compute the
> volume value of each element, none of those are zeros!
>
> Also, I am using the second order Hex27 element, can this be caused by the
> unconvergence of newton iteration? Thanks a lot.
>


Yeah... the Newton iteration inside FE::inverse_map() appears to have
failed, I suppose this could happen if the Newton method evaluates the
Jacobian at a point outside the element, but I'm a little surprised the
determinant was *exactly* 0.0.

BTW, a non-zero volume() doesn't guarantee a pointwise non-singular
Jacobian (especially for a complicated element like the Hex27).  It just
depends on where one is evaluating it, i.e. quadrature points vs. arbitrary
points within the element.

It would be helpful if you could post the node positions of the element in
question (in the libmesh node ordering for Hex27).  That way we could take
a closer look.

-- 
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

Reply via email to