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:
Oh, nvm it’s only a few that are hard to find:
How did you verify, that those exist?
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))[1]; 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][1])
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.
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:
https://github.com/JuliaGeometry/GeometryBasics.jl/blob/master/src/geometry_primitives.jl#L24
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[1]],cols[c[2]],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.
MeshCore can derive these adjacency relationships. Feel free to borrow.
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.
