JuAFEM mesh issue

[bocc]

Hello,

sorry in advance for posting this here, as this is most likely not a JuAFEM issue, but a problem with an external mesh. Also apologies for not being able to create a copy-pastable MWE. Any insight is appreciated regardless.

I wanted to be able to use Fenics meshes, so I created a package to import them. It is very limited at the moment (eg. only tetrahedral meshes), but should be sufficient for what I want to do, which is solving the heat equation on a gear mesh, as a proof of concept / learning example. So I can say

grid = import_grid(path,output=:juafem_grid)

and it returns a Grid with the boundary matrix defined, but no facesets. Adding them to the mesh above can be done by

addfaceset!(grid, "bottom", x -> norm(x[3]) ≈ 0.025)
addfaceset!(grid, "top", x -> norm(x[3]) ≈ 0.413226)

This grid should be usable for the heat equation example almost directly, only by setting dim to 3, and handling the boundary conditions differently:

ch = ConstraintHandler(dh);

dbc_bottom = Dirichlet(:u, getfaceset(grid,"bottom"), (x,t) -> 0)
dbc_top = Dirichlet(:u, getfaceset(grid,"top"), (x,t) -> 0)

add!(ch,dbc_bottom)
add!(ch,dbc_top)

close!(ch)

update!(ch, 0.0)

The issue is that trying to assemble the stiffness matrix triggers a det(J) is not positive positive error. If I try it with JuAFEM’s own

grid = generate_grid(Tetrahedron, (10,10,10))

with the same BCs, it works perfectly. I don’t expect anyone to try and reproduce this, as that would be quite time consuming, but maybe someone has an idea just by reading this.

My guess would be (assuming I did everything right) is that the Fenics representation doesn’t observe the same face orientation rules as JuAFEM does, because every row in the connectivity matrix is in ascending order (while grid_generator creates ‘unordered’ tetrahedrons, probably with common faces having the same orientation in both tets).

Thanks for reading this far

[koehlerson]

I reproduced your issue.

Instead of testing it this way, I’d recommend to setup a test script, where you load your grid (“gears.xml”) and then try to export the dirichlet BC you want to set, see https://github.com/KristofferC/JuAFEM.jl/pull/269

so something like

grid = import_grid("gear.xml",output=:juafem_grid)
dh = DofHandler(grid)
push!(dh, :u, 1) 
close!(dh)

ch = ConstraintHandler(dh)
addfaceset!(grid, "bottom", x -> norm(x[3]) ≈ 0.025)
addfaceset!(grid, "top", x -> norm(x[3]) ≈ 0.413226)
dbc_bottom = Dirichlet(:u, getfaceset(grid,"bottom"), (x,t) -> 0)
dbc_top = Dirichlet(:u, getfaceset(grid,"top"), (x,t) -> 0)
add!(ch,dbc_bottom)
add!(ch,dbc_top)
close!(ch)

vtk = vtk_grid("test_grid", grid)
vtk_point_data(vtk, ch)
vtk_save(vtk) 

in a test folder in your repo and people can then reproduce this easily.

I copy pasted this from 2 sources, so be aware of bugs in the code snippet above. Just try to understand the reference in the issue.

Anyhow, I tested the heat equation example and you have already in the first (and I think also in the second) element a det(J) of 1e-6, which is fairly close to zero. So, seems like a systematic error. Closely after, so one of the first elements, already results in a negative det(J) -1.7208e-6. Therefore, some elements are oscillating around 0 in terms of their det(J). You can analyse this in detail with Debugger.jl

Unfortunately, I handled my unstructured grids by gmsh, therefore I’m not such a great help, but what I do know is, that you find here some references about the global-local indexing of 2d tets
https://github.com/KristofferC/JuAFEM.jl/issues/277

and as far as I remember, the node number ordering coincides with the one from gmsh: Gmsh 4.10.5

[koehlerson]

Actually, nevermind. The test script yields a result, which should be expected I guess?

Should the mesh look like this? The elements on the “top” and “bottom” teeth surface + “top” and “bottom” base circle look a bit messed up

[chakravala]

If your only issue is a negative volume on a stiffness assembly, then why not just use the absolute value?

[kristoffer.carlsson]

You would get normals pointing inwards then so it is best to error since otherwise you will likely get bogus answers when you apply a surface load.

[chakravala]

I agree that orientation matters, that’s why I created a tensor algebra package Grassmann.jl to handle it for me, I’m just trying to offer a possible work around to the OP if it’s indeed only a matter of the stiffness assembly, where orientations should all have the same sign.

[kristoffer.carlsson]

The check is within JuAFEM (the library used) so it isn’t trivial to override.

[chakravala]

In my finite element assembly Adapode.jl, I actually fully switched to using the absolute value of volumes (for things like stiffness assembly) because for an embedded triangle inside a tetrahedral mesh, the volume comes out as a linear combinations of planes (basically isomorphic to a quaternion value). Since i need the volume to be a scalar and not a quaternion or other submanifold, it actually turns out to be correct to take the absolute value of the calculated volume. The volume of a full dimensional simplex is an n-submanifold, which for an embedded triangle in a higher dimension comes out as a quaternion for example. By taking the absolute value of the n-submanifold you get the scalar volume of the element. They should all be having the same orientation anyways, so it makes sense to be working with absolute value for stiffness assembly. That’s why I would recommend switching to absolute value of the volume for stiffness assembly in general. It’s mathematically coherent and correct to do so, although you may need orientation elsewhere.

[Leooop]

@koehlerson

Unfortunately, I handled my unstructured grids by gmsh, therefore I’m not such a great help, but what I do know is, that you find here some references about the global-local indexing of 2d tets

I’m currently thinking about importing Gmsh meshes into the JuAFEM framework. Could you point to a repo containing such code ?

[PetrKryslUCSD]

Gridap has gmsh import.

[Leooop]

Thanks. I’m going to dig into it!

[koehlerson]

https://github.com/koehlerson/Catalyst/blob/master/src/Parser.jl

This will read a .msh and returns a JuAFEM.Grid. Only tested with tets. In general the parsing of the domain itself should work for any kind of discretization if you change https://github.com/koehlerson/Catalyst/blob/master/src/Parser.jl#L147. However, this was my very first project in Julia so it’s really scrub code. If you’re interested to generalize this/extend it, then feel free to write me.

Other than that, you can also use the Gmsh API from https://github.com/JuliaFEM/Gmsh.jl

[Leooop]

Thanks a lot, I’ll definitely write to you it I need any guidance !

[koehlerson]

To come back to OPs problem:

I tested for you the simpler mesh box_with_dent.xml, which looks with a quite similar test script (different BCs) something like this

When I try to apply the heat example to this grid, you get:
ERROR: LoadError: ArgumentError: det(J) is not positive: det(J) = -1.0

Therefore, you are doing something systematically wrong when parsing the .xml, probably related to node ordering

@bocc Could you maybe explain how FEniCS handles node ordering/ connectivity matrix by an example?

[bocc]

Hello Maximilian,

thank you for your insights, they came just at the right moment, as right now I am looking into a linear elasticity problem on a 2d quadrilateral grid that came from Abaqus, and I ran into the same negative Jacobi determinants issue. I will look into gmsh, and yes, that gear and it’s top and bottom surfaces should look like this.

In the meantime, I suspected that the problem is node ordering, so I made this example to reproduce it:

using JuAFEM, Tensors

# ccw
c1 = [(0.0,0.0), (0.0,1.0), (1.0,1.0), (1.0,0.0)]
# cw
c2 = c1[end:-1:1]
# X
c3 = [(0.0,0.0), (1.0,1.0), (0.0,1.0), (1.0,0.0)]


function cv_reinit(cc)
    coords = [Tensor{1,2,Float64}(x) for x in cc]

    dim = 2
    ip = Lagrange{dim, RefCube, 1}()
    qr = QuadratureRule{dim, RefCube}(2)
    cellvalues = CellScalarValues(qr, ip);

    reinit!(cellvalues, coords)

    cellvalues
end

cv_reinit(c1)
cv_reinit(c2)
cv_reinit(c3)

The Jacobi det is positive only in the clockwise case. I think this should be clarified (perhaps on a case by case basis, gmsh is fine, invert order for Abaqus, etc) if we are to use imported meshes with JuAFEM.

[koehlerson]

Hey there,

ah yes good catch. The gear question popped up, because I was wondering initially, if it’s maybe caused by those strongly stretched elements. However, as soon as I tested the simpler .xml I was sure that it’s coming from something else.

I tried to change https://github.com/bocc/Femutils.jl/blob/master/src/import.jl#L78 to p[end:-1:2], however I still get a negative det(J). Any other idea?

[bocc]

I think the element vertices are neither cw nor ccw here, like in the x case above, so just mirroring the order doesn’t help. It seems to me that some software store their elements in ad-hoc order, like ascending node indices, and I assume they take care of ordering at a later step. The transformation being done is described here quite well, it seems reasonable that there is only one correct order. If the nodes are in reverse order, then this involves a ‘reflection’, so a negative Jacobi det is kind of expected.

So it seems to me that it is sufficient to figure out node orders, by for example selecting an internal point in each face, and making sure -internal point to edge first node- and edge cross products have the same sign/direction for all faces. I will try to do this for my current 2d mesh, and then see if it transfers to more complicated elements.

[chakravala]

If your elements are simplicies, it should be sufficient to apply a single transposition of the indices to flip the orientation of a particular element.

[koehlerson]

what do you mean with a single transposition of indices?

I tried to change https://github.com/bocc/Femutils.jl/blob/master/src/import.jl#L78 to p[end:-1:2] , however I still get a negative det(J).

is somewhat what I understand under transposition of the indices, p stores the node indices

[chakravala]

@koehlerson

p[end:-1:2]

That’s reverse, which is not necessarily an odd number of transposition. A reverse of an odd set can be an even number of transpositions, leaving the orientation unchanged. Try restricting to a single transposition.

Cyclic permutation | Transpositions