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. One case it failed is as follows:
pt[186] = (-11.545647,22.442960,12.336682)
Elem id = 449, and volume is 383.684010
Element Nodes are :
node 0 = (-18.750000,7.911975,7.911975)
node 1 = (-18.750000,3.750062,7.806584)
node 2 = (-18.750000,3.834322,12.613942)
node 3 = (-18.750000,8.099410,12.760054)
node 4 = (-30.163472,12.728141,12.728141)
node 5 = (-31.084410,6.216985,12.942030)
node 6 = (-29.400453,6.012310,19.778966)
node 7 = (-28.660566,12.380462,19.504553)
node 8 = (-18.750000,5.831019,7.859280)
node 9 = (-18.750000,3.792192,10.210263)
node 10 = (-18.750000,5.966866,12.686998)
node 11 = (-18.750000,8.005692,10.336015)
node 12 = (-24.456736,10.320058,10.320058)
node 13 = (-24.917205,4.983523,10.374307)
node 14 = (-24.075226,4.923316,16.196454)
node 15 = (-23.705283,10.239936,16.132304)
node 16 = (-30.623941,9.472563,12.835085)
node 17 = (-30.242432,6.114647,16.360498)
node 18 = (-29.030509,9.196386,19.641760)
node 19 = (-29.412019,12.554301,16.116347)
node 20 = (-18.750000,5.898942,10.273139)
node 21 = (-24.686971,7.651791,10.347183)
node 22 = (-24.496216,4.953420,13.285381)
node 23 = (-23.890255,7.581626,16.164379)
node 24 = (-24.081009,10.279997,13.226181)
node 25 = (-29.827225,9.334474,16.238423)
node 26 = (-24.288613,7.616708,13.255781)
On Mon, Jul 18, 2016 at 11:48 AM, John Peterson <jwpeter...@gmail.com>
wrote:
>
>
> 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
