# Makie bug with tetrahedral meshes

In Makie.jl there is a bug with displaying tetrahedral meshes. Using the example from TetGen.jl

Note that at the corners there are missing edges and faces in the plot.

However, I have verified that those edges/faces/simplices are part of the mesh data.

Therefore, it’s a bug with Makie for displaying the tetrahedral mesh, @sdanisch ?

Looks correct to me:

I might have fixed this issue “accidentally” in the newest release of GeometryBasics, so make sure you’re on the newest version of that Oh, nvm it’s only a few that are hard to find: How did you verify, that those exist?

1 Like

Verified the existence with method I described in Grassmann.jl A\b 3x faster than Julia's StaticArrays.jl

First, you need to convert the mesh into `Grassmann` format:

``````using Grassmann
function mesh_init(m::GeometryBasics.Mesh)
c,f = coordinates(m),faces(m)
s = size(eltype(f)); V = SubManifold(ℝ^s)
p = ChainBundle([Chain{V,1}(1.0,k...) for k ∈ c])
t = ChainBundle([Chain{p,1}(SVector(k)) for k ∈ f])
return (p,t)
end
``````

Then define the following `checkpoint` for determining the tetrahedron containing a point:

``````function checkpoint(P::Chain{V,1},t) where V
p = Manifold(t)
for i ∈ 1:length(t)
P ∈ Chain{V,1}(p[value(t[i])]) && (return i)
end
return 0
end
``````

Then you can determine whether there exists a tetrahedron with point `P` or not (returns 0 if not).

``````julia> p,t = mesh_init(result) # result from TetGen
(Λ¹⟨++++⟩×211, Λ¹Λ¹⟨++++⟩×211×566)

julia> checkpoint(Chain{SubManifold(ℝ^4),1}(1,2,2,2),t)
122
``````

Note that I am using homogeneous coordinates with (1,x,y,z) coordinates for the points, so you will have to use 4 dimensional points with an additional 1 coordinate.

You will find that indeed the mesh data contains simplices having points in the problematic regions.

PS, it would be great if `Makie` could natively support `Grassmann` so I don’t need `GeometryBasics` as a compatibility layer. Instead of `GeometryBasics.Point` I prefer using `Grassmann.Chain{V,1}`, because I need that type in order to do my calculations. As you can see, I am using homogeneous coordinates with `Chain{V,1}` types in `Grassmann` to represent mesh data. Currently, I need to convert this mesh data back into `GeometryBasics.Point` for plotting (using cache), which works fine, but I would much prefer being able to work with `Grassmann.Chain{V,1}` directly without conversions.

``````GLMakie.mesh(t, color=(:blue, 0.1), transparency=true)
wireframe!(ans[end])
``````

As you can see, it’s possible to use `mesh` function on `Grassmann` data as well, but it requires the conversion process with caching I mentioned above. I’d prefer to have `Makie` work with `Grassmann` point data directly instead of converting to `GeometryTypes.Point`.

Interesting… Can you try the same with:

``````tetrahedral_mesh = ...
triangular_mesh = triangle_mesh(tetrahedral_mesh)
... figure out if triangle is present ...
``````

Makie needs to triangulate the mesh before displaying. Since there aren’t many places that could just ignore a few triangles after that step, I’m pretty sure the triangulation is where triangles get skipped.

1 Like

Will try that out, but need to modify and extend my programming a bit for embedded manifolds.

Taking a closer look though, I am noticing a few more problems with the triangulation displayed.

There is definitely an issue with the `triangle_mesh` algorithm, I’ll investigate it a bit further later.

Thank you Yeah, this shouldn’t really come as a surprise… This is backed by some old generic algorithms, and you’re probably one of the first to really use it^^
This is where the magic should happen:

1 Like

There likely needs to be a specialization for `SimplexFace` of that method, since it needs to create a closed mesh from a simplex, while NgonFace tries to create an n-vert polygon.

Alright, I made a new commit to `Grassmann` which allows me to check if a point is inside of an embedded simplex of a higher dimensional manifold (e.g. check if point is inside affine triangle).

Using this method, I found that indeed the triangles are missing from the `triangulate_mesh` result.

Not only are triangles missing, some of them are duplicates! There are triangles which appear twice.

Most likely, I believe that your algorithm has an issue with the orientation of simplices. I have made an algorithm for triangulation of a simplicial complex in `Grassmann` but it’s only for local geometries. It is one of my tasks to generalize this triangulation to global mesh geometries also. Therefore, I will also spend some more time later to investigate the issues with your algorithm to see what’s wrong.

Got an algorithm which directly computes all the edges from a simplex mesh (of any dimension > 1)

``````column(t,i) = getindex.(value(t),i)
columns(t) = column.(Ref(value(t)),Grassmann.list(1,ndims(t)))
function edges(t,cols=columns(t))
np,N = length(points(t)),ndims(Manifold(t))
A,M = spzeros(np,np),points(t)(Grassmann.list(N-1,N)...)
for c ∈ Grassmann.combo(N,2)
A += sparse(cols[c],cols[c],1,np,np)
end
f = findall(x->x>0,LinearAlgebra.triu(A+transpose(A)))
[Chain{M,1}(SVector{2,Int}(f[n].I)) for n ∈ 1:length(f)]
end
``````

Resulting in the correct line segments for the mesh

However, your method involves first computing all the triangles first, and then the line segments. So this algorithm doesn’t directly translate over into your implementation. When I come up with a method similar to yours based on first getting the triangles, I’ll make a pull request.

2 Likes

MeshCore can derive these adjacency relationships. Feel free to borrow.

1 Like

Anyone know how to obtain the boundary triangles from the `GeometryTypes.Mesh` or `TetGen`?

How have I missed the rather wonderful progress you have made with Finite Element methods in Julia?
Last time I looked at FE in Julia things were at quite an early stage.  1 Like