negative Jacobian tests in 2D

[WilkAndy]

Refs #9740

I need help with this. @jwpeterson ?

This PR contains a single 2D input file that illustrates:
(1) MOOSE does not fail if an element is inverted, if the kernels use_displaced_mesh=false
(2) MOOSE does fail (correctly) if an element is deformed such that Jacobian=0, if the kernels use_displaced_mesh=true

But, i can't get MOOSE to fail if:
(3) an element is deformed such that Jacobian < 0, if the kernels use_displaced_mesh=true

There is already a 3D test (in the same directory) that illustrate (1), (2) and (3) work in 3D.

The question is: why is Libmesh not producing an error in 2D when Jacobian < 0. From line 736 of libmesh/src/fe/fe_map.C i believe that it should.

[WilkAndy]

Is LIBMESH_DIM=2 or LIBMESH_DIM=3 for this 2D input file?

[moosebuild]

Job Documentation on cf2d540 wanted to post the following:

View the site here

This comment will be updated on new commits.

[bwspenc]

Thanks for setting this test up, @WilkAndy! I'd love to get to the bottom of this.

[jwpeterson]

Is LIBMESH_DIM=2 or LIBMESH_DIM=3 for this 2D input file?

LIBMESH_DIM is a compile-time constant which is always equal to 3. I'm not even sure if libmesh currently compiles if you set LIBMESH_DIM==2, and I'm almost positive it wouldn't run. This constant primarily controls the size of the Point objects.

Your question about about mixed dimensions got me thinking though. When computing the Jacobian determinant for a 2D element (potentially) living in 3D space, we are talking about the determinant of a non-square matrix, which of course does not exist. So what we do, on line 797 of fe_map.C, is compute the determinant of the square matrix T' T which is formed by multiplying the rectangular Jacobian by its transpose. And this determinant can never be negative, see e.g. the comments from the code below.

// Compute the Jacobian.  This assumes a 2D face in
// 3D space.
//
// The transformation matrix T from local to global
// coordinates is
//
//         | dx/dxi  dx/deta |
//     T = | dy/dxi  dy/deta |
//         | dz/dxi  dz/deta |
// note det(T' T) = det(T')det(T) = det(T)det(T)
// so det(T) = std::sqrt(det(T' T))

This determinant can of course still be zero, which is consistent with what @WilkAndy reported above. I suspect this is why we never see negative Jacobians for 2D problems no matter how badly we invert elements, but I don't really know what the fix is, short of somehow detecting that the last row of T is 0 and falling back on the 2D-element-living-in-2D code path somehow. Maybe @roystgnr has some ideas... whatever we do cannot be expensive, since compute_single_point_map() gets called billions of times...

[roystgnr]

I'm not even sure if libmesh currently compiles if you set LIBMESH_DIM==2, and I'm almost positive it wouldn't run.

BuildBot actually used to test --enable-2D-only and even --enable-1D-only regularly; they might not have succumbed to bitrot since then yet...

I don't really know what the fix is, short of somehow detecting that the last row of T is 0 and falling back on the 2D-element-living-in-2D code path somehow.

Even that wouldn't be a valid fix, I'm afraid, because that would be a false positive in a supported use case! We deliberately allow users to wrap 2D manifolds around in 3D space. If one of @pbauman's stretchy membrane simulations has a parachute where one element ends up flat-but-exactly-inverted (which seems improbable a priori but might happen frequently in symmetric initial conditions) then we want that calculation to chug on, not error out.

Maybe @roystgnr has some ideas... whatever we do cannot be expensive, since compute_single_point_map() gets called billions of times...

If this is for sanity checking then we don't need to put the check in every single mapping computation; you could do a single sanity-check-only sweep over elements before assembly even begins, and then not bother to re-test at every point in every residual evaluation. And if you know a priori that you're not in an "element legitimately flipped" use case then you can just run the 2D jac calculation in user space and test for non-positive.

[jwpeterson]

If this is for sanity checking

Ideally we'd like to be able to detect "negative Jacobian" (whatever this means for non-square Jacobian matrices) elements in the middle of e.g. solid mechanics computations with a displaced mesh, and throw an exception that would then allow us to cut the time step and retry the most recent solve.

We could get part of the way there by the pre- or post-step check of the elements that you suggest, backing up an entire time step if an inverted element was detected. But even this would require a reliable algorithm for detecting negative Jacobian, non-planar manifold elements... which I am not sure how to even define at the moment.

[jwpeterson]

@WilkAndy I have added a quick and dirty "2D_in_3D" branch in my libmesh fork that detects when an Elem is xy-planar (even though it "lives" in 3D) and instead uses the standard 2D formula for computing the element Jacobian determinants. When I uncomment your new test and run ./run_tests --re=jacobian_2D, the misc/jacobian.jacobian_2D_negative test prints:

Time Step  2, time = 1.5
                dt = 0.7
  Elem Information
   id()=0, unique_id()=4, processor_id()=0
   type()=QUAD4
   dim()=2
   n_nodes()=4
    0  Node id()=0, processor_id()=0, Point=(x,y,z)=(       0,        0,        0)
    DoFs=(0/0/0) (1/0/0) (1/1/1) 
    1  Node id()=1, processor_id()=0, Point=(x,y,z)=(    -0.5,        0,        0)
    DoFs=(0/0/1) (1/0/2) (1/1/3) 
    2  Node id()=2, processor_id()=0, Point=(x,y,z)=(    -0.5,        1,        0)
    DoFs=(0/0/2) (1/0/4) (1/1/5) 
    3  Node id()=3, processor_id()=0, Point=(x,y,z)=(       0,        1,        0)
    DoFs=(0/0/3) (1/0/6) (1/1/7) 
   n_sides()=4
    neighbor(0)=NULL
    neighbor(1)=NULL
    neighbor(2)=NULL
    neighbor(3)=NULL
   hmin()=0.5, hmax()=1.11803
   volume()=0.5
   active()=1, ancestor()=0, subactive()=0, has_children()=0
   parent()=NULL
   level()=0, p_level()=0
   refinement_flag()=DO_NOTHING
   p_refinement_flag()=DO_NOTHING
   infinite()=0
   DoFs=
ERROR: negative Jacobian -0.125 at point (x,y,z)=(-0.105662, 0.211325,        0) in element 0

which I believe is close to (the value is off by a factor of 2) the expected result?

Running the rest of the test suite with this version of libmesh also results in quite a few other failures:

1221 passed, 66 skipped, 0 pending, 227 FAILED

that would have to be addressed on a case by case basis. So I'm not sure where we go from here, but this is at least a start.

[jwpeterson]

@gardnerru was also interested in this thread.

[jwpeterson]

Oops, there was an extra sqrt(det) in my first attempt... Fixing that makes the MOOSE tests look a lot better :-)

MOOSE: 1451 passed, 58 skipped, 0 pending, 6 FAILED
MODULES: 2154 passed, 178 skipped, 0 pending, 8 FAILED

I also updated the link above so that no one grabs the wrong commit by mistake...

[WilkAndy]

Thanks for looking into this everyone.

Yes, the above ERROR: negative Jacobian -0.125 at point (x,y,z)=(-0.105662, 0.211325, 0) in element 0 is what i expected to see.

Where can we see the test results (6 moose failures, 8 modules failures)? If i have to update libmesh i will be compiling all day (literally) on my little MacBookAir.

I'll leave it up to you guys to decide what to do here. I'm not just being lazy or passing-the-buck - this is a pretty fundamental issue and i don't know what's best from a MOOSE perspective. Do we want MOOSE to fail in 2D simulations when elements get inverted (i think so, but will leave it up to you)? My initial reason for this Issue was that i noticed MOOSE suddenly wasn't failing (in 3D) with use_displaced_mesh=false and that was VERY GOOD for me :-) So i certainly don't want to go back to the old days of MOOSE failing for inverted use_displaced_mesh=false elements.

[jwpeterson]

Where can we see the test results (6 moose failures, 8 modules failures)? If i have to update libmesh i will be compiling all day (literally) on my little MacBookAir.

I just ran them locally, didn't push a PR or anything for testing.

Other than your new test (which only failed due to a slightly different expected error message), the remaining MOOSE tests that failed for me were all mortar related:

misc/jacobian.jacobian_2D_negative.................................................. FAILED (NO EXPECTED ERR)
mortar/2d.non-conforming.................................................................... FAILED (EXODIFF)
mortar/2d.conforming......................................................................... FAILED (ERRMSG)
mortar/2d.conforming_two_var................................................................. FAILED (ERRMSG)
mortar/2d.conforming-2nd-order............................................................... FAILED (ERRMSG)
mortar/2d.equalgradient...................................................................... FAILED (ERRMSG)

I haven't looked into it in detail, but I think those tests may be using a mesh with an inverted element somehow, and we just never caught it because of the problem with the 2D Jacobian calculation.

[WilkAndy]

And were the modules failures in tensor mechanics?

[jwpeterson]

A couple of them were, yeah:

combined/test:gap_heat_transfer_htonly.RZ................................................... FAILED (EXODIFF)
combined/test:inelastic_strain/elas_plas.elastic_plastic_sm................................. FAILED (EXODIFF)
tensor_mechanics/test:CylindricalRankTwoAux.test............................................ FAILED (EXODIFF)
combined/test:gravity.rz_quad8.............................................................. FAILED (EXODIFF)
phase_field/test:mesh.mortarperiodic_hex8.................................................... FAILED (ERRMSG)
phase_field/test:mesh.mortarperiodic_hex27................................................... FAILED (ERRMSG)
tensor_mechanics/test:generalized_plane_strain.generalized_plane_strain_squares.............. FAILED (ERRMSG)
combined/test:contact_verification/patch_tests/single_pnt_2d.mu_0_2_kin_sm.................. FAILED (CSVDIFF)

The plane_strain_squares is a crash due to a negative Jacobian:

  Elem Information
   id()=0, unique_id()=22, processor_id()=0
   type()=QUAD4
   dim()=2
   n_nodes()=4
    0  Node id()=0, processor_id()=0, Point=(x,y,z)=(     0.5,        1,        0)
    DoFs=
    1  Node id()=1, processor_id()=0, Point=(x,y,z)=(    0.75,        1,        0)
    DoFs=
    2  Node id()=2, processor_id()=0, Point=(x,y,z)=(    0.75,      0.5,        0)
    DoFs=
    3  Node id()=3, processor_id()=0, Point=(x,y,z)=(     0.5,      0.5,        0)
    DoFs=
   n_sides()=4
    neighbor(0)=NULL
    neighbor(1)=1
    neighbor(2)=2
    neighbor(3)=NULL
   hmin()=0.25, hmax()=0.559017
   volume()=0.125
   active()=1, ancestor()=0, subactive()=0, has_children()=0
   parent()=NULL
   level()=0, p_level()=0
   refinement_flag()=DO_NOTHING
   p_refinement_flag()=DO_NOTHING
   infinite()=0
   DoFs=
ERROR: negative Jacobian -0.03125 at point index 0 in element 0

The CylindricalRankTwoAuxone is a relatively small exodiff.

[WilkAndy]

Wow, i have no idea what that plane_strain_squares is trying to do! The original author is @hchen139 - perhaps they can shed light on why it fails?

[pbauman]

Even that wouldn't be a valid fix, I'm afraid, because that would be a false positive in a supported use case! We deliberately allow users to wrap 2D manifolds around in 3D space. If one of @pbauman's stretchy membrane simulations has a parachute where one element ends up flat-but-exactly-inverted (which seems improbable a priori but might happen frequently in symmetric initial conditions) then we want that calculation to chug on, not error out.

Just to chime in on this point, it definitely can happen. When we first started playing with those meshes, we actually had to do some cleanup because some elements didn't have a consistent outward normal (some pointed inside the parachute, some outside). That required some reorienting of some elements.

[jwpeterson]

My current thinking on this issue is: for a 2D element living in 3D it is not possible to reliably detect an inverted/bad element (via Jacobian determinant sign) without some kind of external information, i.e. a "reference normal", being provided by the user. We can't use the element's own surface normal, because a twisted element will have a correspondingly inverted normal and its area would be self-consistent with this. Therefore, we might as well go with the current implementation which works in the general case, but which we now know is too permissive in that it allows all kinds of tangled and badly-shaped elements.

We would then additionally provide a sanity check like FE::signed_area(const Point & reference_normal); where the user provides the reference normal. Application codes like MOOSE could then loop over the elements, call this function as a post-processing step, and decide what to do (retry the previous timestep, etc.) if inverted elements are detected. This sanity check would work for 2D meshes living in any plane, and could also work for things like the surfaces of spheres where the reference outward normal direction is a known function of (x,y,z).

Not sure if this will fix the case that @WilkAndy was originally trying to address since you won't be able to throw an exception immediately (i.e. in the middle of computeResidual()) when a negative Jacobian is detected, but I think it would be better than nothing...

[roystgnr]

We can't use the element's own surface normal, because a twisted element will have a correspondingly inverted normal and its area would be self-consistent with this.

Well, a "twisted" quad would have a normal at some quadrature points which was opposite the normal at others, but I would expect tensor mechanics simulations on quad meshes to risk wholly-inverted elements, not just twisted elements, and on TRI3 meshes you would only ever get wholly-inverted elements. So a sanity check that could only detect twisted quad/TRI6 elements would be little more than false confidence.

Therefore, we might as well go with the current implementation which works in the general case, but which we now know is too permissive in that it allows all kinds of tangled and badly-shaped elements.

We would then additionally provide a sanity check like FE::signed_area(const Point & reference_normal); where the user provides the reference normal.

So, a third option just now occurred to me: should MOOSE (or libMesh, for that matter) provide applications a way to save the current normals for all top-level elements, and then reuse those saved normals (looking them up via top_parent() in the case of AMR) for sanity checking later? When the user reads in their mesh, whether it's something generated a priori or a restart mesh from a previous (sanity-checked!) run, they ought to know that that mesh has no inverted or twisted elements. And if it has no inverted or twisted elements initially, then they can rely on those initial normal vectors to determine whether each element or any of its descendants has been inverted or twisted subsequently.

This would only work if you are solving meshes with 2D manifolds where the element movement is constrained to remain within the manifold (because if not then e.g. one of @pbauman's parachutes can "fold" around a stiffener and that's not non-physical) and where the manifold doesn't curve too much (because if an element can "slide" around a 90-degree bend then old_normal * new_normal < 0 would throw a false positive), so it might not make good default behavior, but it ought to be a complete sanity check for the majority of users.

[jwpeterson]

Well, a "twisted" quad would have a normal at some quadrature points which was opposite the normal at others

You're right, bad terminology by me. I was thinking primarily of the wholly inverted case while writing this.

So, a third option just now occurred to me: should MOOSE (or libMesh, for that matter) provide applications a way to save the current normals for all top-level elements, and then reuse those saved normals

This is an interesting idea. It would allow initial meshes which were "bad" for whatever reason (this sometimes happens when people build meshes by hand or try to implement a reader for their own mesh format) but I suspect that to be a small percentage of use cases. Storing the extra information in a user-friendely, parallel, elegant, way consistent with the rest of the library sounds like it might be a bit of a chore though.

[permcody]

Storing the extra information in a user-friendely, parallel, elegant, way consistent with the rest of the library sounds like it might be a bit of a chore though.

Maybe, maybe not: libMesh/libmesh#1439

[bwspenc]

Currently, at least for mechanics simulations, all of the 2D problems that we solve in MOOSE have 2D elements in a 2D domain. These models also are currently always defined in the xy plane, although I imagine they could be defined in other planes. This is the case that's really important for us to have an inversion check for currently. Detecting inversion of a 2D element in 3D space is a far more exotic problem, and probably far less likely to occur because elements can distort out of plane when they are pushed in ways that tend to make them invert. We also don't even have shell elements in MOOSE yet, so it's not even relevant to us yet (although we do have plans to add such a capability in the near future).

I just looked at some of our 2D meshes, and I think we usually (maybe always) have 'num_dim=2' set in the exodus mesh. Couldn't we just have some way to either detect that the mesh is living in 2D (by looking at num_dim in the case of exodus meshes) or have an API to tell libMesh that the mesh is in 2D, and then modify the inversion check to take that into account? I envision some kind of enum that stores the dimensionality of the mesh, as well as information about the plane (or axis for the 1D case) it lives in. This could be used for inversion checks of 1D elements in 1D space as well. It could have values like 1DX, 1DY, 1DZ, 2DXY, 2DXZ, 2DYZ, 3D. With that information available, I imagine that we would be able to do efficient inversion checks within fe.reinit() for all these cases.

[permcody]

@jwpeterson - What's the final resolution on this, or is there one? Do we need to talk about this some more. This is nearly two months old now.

[jwpeterson]

I think the best solution is to not change the current behavior of FE::compute_single_point_map() since that would negatively impact the performance of every MOOSE application and probably cause quite a few diffs along the way.

Instead, we could add a sanity check based on either a specified reference normal, as described in my previous comment, or something more heavy-handed (like described in Roy's followup) that would alert users when an element's Jacobian goes bad. This sanity check could eventually be used to cut back timesteps or in some other way allow user's to recover from previous solutions once elements become inverted.

[WilkAndy]

It's late on Friday afternoon (=beer) and i just noticed this negative-Jacobian thing again. Thanks for all your contemplating, everyone. I hope you find it fun to fix, when you eventually get round to it :-)

[permcody]

Closing due to inactivity but we need to fix this soon!

[WilkAndy]

:-(