# JuAFEM mesh issue

**URL:** <https://discourse.julialang.org/t/juafem-mesh-issue/41467>\
**Category:** Numerics\
**Tags:** juafem\
**Created:** [June 15, 2020, 8:08pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467 "2020-06-15T20:08:54Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![bocc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bocc/32/7546_2.png) [@bocc](https://discourse.julialang.org/u/bocc)\
**Post date:** [June 15, 2020, 8:08pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/1 "2020-06-15T20:08:54Z")

</div>

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](https://github.com/bocc/Femutils.jl) 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](https://fenicsproject.org/pub/data/meshes/gear.xml), as a proof of concept / learning example. So I can say

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

```julia
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](https://kristofferc.github.io/JuAFEM.jl/dev/examples/heat_equation/) example almost directly, only by setting `dim` to `3`, and handling the boundary conditions differently:

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

```julia
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 😉

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [July 17, 2020, 11:19am UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/2 "2020-07-17T11:19:15Z")

</div>

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](https://github.com/KristofferC/JuAFEM.jl/pull/269)

so something like

```julia
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](https://github.com/JuliaDebug/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](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](https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering)

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [July 17, 2020, 11:30am UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/3 "2020-07-17T11:30:49Z")

</div>

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

 ![test_1](https://global.discourse-cdn.com/julialang/original/3X/d/4/d4f3b495c8d2db116581f023e35160795bca512e.png) ![test_2](https://global.discourse-cdn.com/julialang/original/3X/4/f/4f8abca664efda23feed083314371ba09868426c.png)

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

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [July 17, 2020, 2:27pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/4 "2020-07-17T14:27:38Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [July 17, 2020, 2:46pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/5 "2020-07-17T14:46:18Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [July 17, 2020, 3:04pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/6 "2020-07-17T15:04:00Z")

</div>

I agree that orientation matters, that’s why I created a tensor algebra package [Grassmann.jl](https://github.com/chakravala/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.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [July 17, 2020, 3:08pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/7 "2020-07-17T15:08:51Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [July 17, 2020, 3:17pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/8 "2020-07-17T15:17:42Z")

</div>

In my finite element assembly [Adapode.jl](https://github.com/chakravala/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.

---

<div class="post-metadata">

**Author:** ![Leooop](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leooop/32/15154_2.png) [@Leooop](https://discourse.julialang.org/u/Leooop)\
**Post date:** [July 17, 2020, 3:24pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/9 "2020-07-17T15:24:21Z")

</div>

> [@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 ? 🙂

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [July 17, 2020, 3:31pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/10 "2020-07-17T15:31:09Z")

</div>

Gridap has gmsh import.

---

<div class="post-metadata">

**Author:** ![Leooop](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leooop/32/15154_2.png) [@Leooop](https://discourse.julialang.org/u/Leooop)\
**Post date:** [July 17, 2020, 3:35pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/11 "2020-07-17T15:35:21Z")

</div>

Thanks. I’m going to dig into it!

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [July 17, 2020, 4:07pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/12 "2020-07-17T16:07:49Z")

</div>

[https://github.com/koehlerson/Catalyst/blob/master/src/Parser.jl](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](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](https://github.com/JuliaFEM/Gmsh.jl)

---

<div class="post-metadata">

**Author:** ![Leooop](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leooop/32/15154_2.png) [@Leooop](https://discourse.julialang.org/u/Leooop)\
**Post date:** [July 17, 2020, 4:27pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/13 "2020-07-17T16:27:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [July 17, 2020, 4:28pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/14 "2020-07-17T16:28:39Z")

</div>

To come back to OPs problem:

I tested for you the simpler mesh [box\_with\_dent.xml](https://fenicsproject.org/pub/data/meshes/box_with_dent.xml), which looks with a quite similar test script (different BCs) something like this

 ![other_test_bc](https://global.discourse-cdn.com/julialang/original/3X/a/f/af1644dce7a699463ac18864187430b9d9eb403a.png)

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?

---

<div class="post-metadata">

**Author:** ![bocc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bocc/32/7546_2.png) [@bocc](https://discourse.julialang.org/u/bocc)\
**Post date:** [July 17, 2020, 5:40pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/15 "2020-07-17T17:40:04Z")

</div>

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [July 17, 2020, 5:59pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/16 "2020-07-17T17:59:34Z")

</div>

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](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?

---

<div class="post-metadata">

**Author:** ![bocc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bocc/32/7546_2.png) [@bocc](https://discourse.julialang.org/u/bocc)\
**Post date:** [July 18, 2020, 9:09am UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/17 "2020-07-18T09:09:20Z")

</div>

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](https://www.iue.tuwien.ac.at/phd/orio/node48.html) 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.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [July 18, 2020, 10:02am UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/18 "2020-07-18T10:02:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![koehlerson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/koehlerson/32/13108_2.png) [@koehlerson](https://discourse.julialang.org/u/koehlerson)\
**Post date:** [July 18, 2020, 12:49pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/19 "2020-07-18T12:49:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [July 18, 2020, 1:07pm UTC](https://discourse.julialang.org/t/juafem-mesh-issue/41467/20 "2020-07-18T13:07:36Z")

</div>

> [@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](https://en.m.wikipedia.org/wiki/Cyclic_permutation#Transpositions)**
>
> A cycle with only two elements is called a transposition. For example, the permutation 
>   
>     
>       
> π
> =
>         
>           
> (
>             
>               
>                 
> 1
>                 
>                 
> 2
>                 
>                 
> 3
>                 
>                 
> 4
>                 
>               
>               
>                 
> 1
>                 
> ...