An interface for meshes

[juliohm]

Hi @SimonDanisch , this is a follow up from yesterday's Discourse post where we shared that it would be nice to have an interface for general meshes, from simple Julia arrays (treated as regular meshes with zero storage for the coordinates) to general mixed-element meshes for various different purposes.

As an initial target, we could try to prove the concept that such interface would enable (1) geostatistical estimation/simulation of properties over meshes with GeoStats.jl, (2) physical simulation with available FEM solvers, and (3) plotting with Makie.jl

I understand that (3) is already done, and that perhaps (2) as well? Regarding (1), I would like to have your feedback on the interface copied/pasted below:

module Meshes

using StaticArrays: MVector

#--------------
# MESH TRAITS
#--------------
ndims(::Type{M}) where M = @error "not implemented"

ctype(::Type{M}) where M = @error "not implemented"

cbuff(m::Type{M}) where M = MVector{ndims(m),ctype(m)}(undef)

isstructured(::Type{M}) where M = Val{false}()

isregular(::Type{M}) where M = Val{false}()

COMPILE_TIME_TRAITS = [:ndims, :ctype, :cbuff, :isstructured, :isregular]

# default versions for mesh instances
for TRAIT in COMPILE_TIME_TRAITS
  @eval $TRAIT(::M) where M = $TRAIT(M)
end

elements(mesh::M) where M = @error "not implemented"

nelms(mesh::M) where M = length(elements(mesh))

#----------------------
# MESH ELEMENT TRAITS
#----------------------
coords!(x::AbstractVector, mesh::M, elm::E) where {M,E} =
  @error "not implemented"

function coords!(X::AbstractMatrix, mesh::M,
                 elms::AbstractVector) where M
  for j in 1:length(elms)
    coords!(view(X,:,j), mesh, elms[j])
  end
end

function coords(mesh::M, elm::E) where {M,E}
  x = cbuff(M)
  coords!(x, mesh, elm)
  x
end

function coords(mesh::M, elms::AbstractVector) where M
  X = Matrix{ctype(M)}(undef, ndims(M), length(elms))
  coords!(X, mesh, elms)
  X
end

coords(mesh::M) where M = coords(mesh, elements(mesh))

vertices(mesh::M, elm::E) where {M,E} = @error "not implemented"

nverts(mesh::M, elm::E) where {M,E} = length(vertices(mesh, elm))

#------------------------
# ELEMENT VERTEX TRAITS
#------------------------
coords!(x::AbstractVector, mesh::M, elm::E, vert::V) where {M,E,V} =
  @error "not implemented"

function coords!(X::AbstractMatrix, mesh::M, elm::E,
                 verts::AbstractVector) where {M,E}
  for j in 1:length(verts)
    coords!(view(X,:,j), mesh, elm, verts[j])
  end
end

function coords(mesh::M, elm::E, vert::V) where {M,E,V}
  x = cbuff(M)
  coords!(x, mesh, elm, vert)
  x
end

function coords(mesh::M, elm::E, verts::AbstractVector) where {M,E}
  X = Matrix{ctype(M)}(undef, ndims(M), nelms(mesh))
  coords!(X, mesh, elm, verts)
  X
end

end # module

It is a quite simple interface that covers all algorithms currently implemented in GeoStats.jl as far as I can tell. It consists of basic information about the dimension of the mesh and type of coordinates, an access function to an iterator of elements (e.g. elements could simply be the range 1:n), and functions to access the coordinates of the elements themselves (some pre-defined baricenter or centroid for different element types), and the coordinates of the vertices of the elements.

Does the current implementation in GeometryBasics.jl (and associated interface) cover this use case already? If not, can we work on it together to create a separate package (I suggest Meshes.jl) with a set of traits without implementation that covers the use case?

After that we could start extending this interface to include volume, surfacearea and other useful properties of elements that are constantly used in algorithms.

[SimonDanisch]

Btw, the whole idea of the metadata system in GeometryBasics was, to easily add caching like that...

[PetrKryslUCSD]

Yes, but this is not caching of local information. For certain incidence relations, the information to compute (and cache) is global.

[SimonDanisch]

Yes, but I also need global metadata, probably don't have any examples right now ;)

[PetrKryslUCSD]

Neat, I'll check it out.

[SimonDanisch]

I'll need to create tests & examples for it - likely doesn't work out of the box ;) Sorry if I was misleading. What I meant was, that the current metadata system has been written to allow exactly this, but it hasn't been added to all types yet.

You think that will mostly cover the needed funcitonality?

[PetrKryslUCSD]

I see. Might be...

[PetrKryslUCSD]

https://github.com/PetrKryslUCSD/MeshCore.jl has been redesigned from the ground up. Meshes are now based on incidence relations and shape collections, and the geometry is an attribute of the shape collection of vertices. Understanding mesh topology as an incidence relation made the conceptual view of the mesh self-consistent and easier to understand, in my opinion at least. The entire spectrum of the (unambiguous) topological relations is now covered.

[juliohm]

Getting back to this issue, finally! I am excited to start reviewing this immense body of work on meshes. I will try to prepare a document as I do this review so that we are all aligned with design choices adopted in our ecosystem, missing features, and possible improvements.

[PetrKryslUCSD]

Sounds good. Just in case this is useful: here's the submitted paper on the https://github.com/PetrKryslUCSD/MeshCore.jl library. 102557_0_art_file_394454_qggdgd_sc.pdf

[juliohm]

@SimonDanisch is it possible to create a 3D regular grid mesh with GeometryBasics.jl? Could you please share a code snippet? I've noticed that the codebase evolved a lot since the last time I looked into it.

[juliohm]

I've started a document this weekend to collect/organize the information shared in this thread: https://docs.google.com/document/d/1eFE4MSwARDJdQ_eX1UYLlN5gdayigjSa2PAOGfJRvsQ/edit?usp=sharing Could you please review it to see if it makes sense?

I am finding it challenging to map some of the concepts in the document to the GeometryBasics.jl implementation. If someone can point out how I can construct 3D meshes (without data) with the current API that would help a lot.

[juliohm]

@SimonDanisch did you have a chance to take a look into this? It would be really nice to be able to construct some simple example meshes to try implement/evolve the interface.

[SimonDanisch]

You can create a Mesh from any AbstractVector{<: PolyTope}. The most common one is FaceView(positions::Vector{AbstractPoint}, faces::Vector{AbstractFace}), which the Mesh constructor has a convenience constructor for: Mesh(positions, faces). I just notice, that I don't have a convenience constructor for FaceView since it's covered by connect and Mesh:

# this is equivalent to the non existing FaceView constructor mentioned above
fv = connect(rand(Point3f0, 10), TriangleFace{Int}[(1, 2, 3), (3, 4, 1)])::FaceView{Triangle}

The FaceView will use the faces to index into the positions to construct the polytopes. E.g.:

@assert fv[1] isa Triangle
m = Mesh(fv)
m[1] == fv[1]

For flat geometries e.g. Lines/Triangles/Quads you use these face type: https://github.com/JuliaGeometry/GeometryBasics.jl/blob/master/src/basic_types.jl#L36 But you can use any arbitrary face type, e.g. a SimplexFace for defining Tetrahedrons, or create a new Face type inhering from AbstractFace for e.g. a Face type with varying elements. You can destructure the mesh with these functions:

faces(m) / faces(fv) == the face vector
coordinates(m) / coordinates(fv) == the positions

[SimonDanisch]

Btw, there are quite a few utitilities in place to keep the memory layout of existing meshes. E.g. consider having an Int array for faces:

positions = rand(Point3f0, 4)
faces = [1, 2, 3, 3, 4, 1]
m1 = Mesh(positions, connect(faces, TriangleFace{Int}))
# or 0 based indices can also stay that way in memory:
faces = [1, 2, 3, 3, 4, 1] .- 1
m2 = Mesh(positions, connect(faces, TriangleFace{OffsetInteger{-1,Int}}))
@assert m2[end] == m1[end]

[sjkelly]

@juliohm I did some work with mesh evolutions using Tetgen and GeometryBasics here: https://github.com/juliageometry/DistMesh.jl . The current limitation is that it is 3D only, and doesn't preserve properties. I can send you my paper and the original thesis if you are interested.

[juliohm]

Thank you all for the updates. I will try to look into it today.

@SimonDanisch is there a helper function to build a N-dimensional regular grid as well? I will play with your code snippets to see if I am able to grasp the functionality.

@sjkelly I'd be happy to learn more. How would you like to share your thesis and paper? Maybe over Zulip/Slack?

[SimonDanisch]

is there a helper function to build a N-dimensional regular grid as well? I will play with your code snippets to see if I am able to grasp the functionality.

I dont know what exactly that would create, so I guess no :D But maybe? This creates a 50x50 quadmesh on a rectangular grid:

julia> t = Tesselation(FRect2D(0, 0, 2, 2), (50, 50))
julia> GeometryBasics.mesh(t, facetype=QuadFace{Int})

[juliohm]

julia> t = Tesselation(FRect2D(0, 0, 2, 2), (50, 50))
julia> GeometryBasics.mesh(t, facetype=QuadFace{Int})

Thank you @SimonDanisch , how to visualize this mesh with Makie.jl? That will help grasping the concepts.

[juliohm]

I tried Makie.mesh and Makie.plot on the object, but the renderer returned an error. Stable version of the packages.

[SimonDanisch]

Lol...

This calls for some tests & bug fixes^^

[SimonDanisch]

Wireframe already works:

m = GeometryBasics.mesh(t, pointtype=Point2f0, facetype=QuadFace{Int})
s = wireframe(m)

[SimonDanisch]

Alright, fixes in: https://github.com/JuliaGeometry/GeometryBasics.jl/pull/62 https://github.com/JuliaPlots/AbstractPlotting.jl/pull/450

[juliohm]

Awesome @SimonDanisch . The wireframe is great actually for visualizing the meshes. I am trying to shape the interface in Meshes.jl and implement it for the meshes in GeometryBasics.jl as planned.

[juliohm]

I think I've reached a point where I can start replacing the spatial domain types in GeoStatsBase.jl by the mesh types defined in GeometryBasics.jl. Currently, I am only using 4 function names from the Meshes.jl proposal: ndims, coordtype, elements, and ecoords!.

Below are proof-of-concept implementations of the interface for (1) general mesh types provided by GeometryBasics.jl, (2) regular grids (a special case with auxiliary specs [origin, spacing, dims]), and (3) point sets (degenerated meshes).

using GeometryBasics
import Meshes

# ----------------------------
# GENERAL MESHES
# ----------------------------
Meshes.ndims(::Type{GeometryBasics.Mesh{N,T,E,V}}) where {N,T,E,V} = N
Meshes.coordtype(::Type{GeometryBasics.Mesh{N,T,E,V}}) where {N,T,E,V} = T
Meshes.elements(m::GeometryBasics.Mesh) = m # Mesh is iterable already
function Meshes.ecoords!(x::AbstractVector, m::GeometryBasics.Mesh, e)
  points = coordinates(e)
  x .= sum(points) ./ length(points)
end

# ------------------------------------------------
# REGULAR GRID (SPECIAL CASE)
# ------------------------------------------------
struct RegularGrid{T,N}
  dims::Dims{N}
  origin::NTuple{N,T}
  spacing::NTuple{N,T}

  # state fields
  mesh::GeometryBasics.Mesh

  function RegularGrid{T,N}(dims, origin, spacing) where {N,T}
    rect = Rect{N,T}(origin, dims.*spacing)
    tess = Tesselation(rect, dims.+1)
    mesh = GeometryBasics.mesh(tess, facetype=NgonFace{2^N,Int})
    new(dims, origin, spacing, mesh)
  end
end
RegularGrid(dims::Dims{N}, origin::NTuple{N,T}, spacing::NTuple{N,T}) where {N,T} =
  RegularGrid{T,N}(dims, origin, spacing)
RegularGrid{T}(dims::Dims{N}) where {N,T} =
  RegularGrid{T,N}(dims, ntuple(i->zero(T), N), ntuple(i->one(T), N))
RegularGrid(dims::Dims{N}) where {N} = RegularGrid{Float32}(dims)

Meshes.ndims(::Type{RegularGrid{T,N}}) where {T,N} = N
Meshes.coordtype(::Type{RegularGrid{T,N}}) where {T,N} = T
Meshes.elements(g::RegularGrid) = Meshes.elements(g.mesh)
Meshes.ecoords!(x::AbstractVector, g::RegularGrid, e) = Meshes.ecoords!(x, g.mesh, e)

# ---------------------------------------------------
# POINT SET (DEGENERATED MESH)
# ---------------------------------------------------
struct PointSet{T,N}
  coords::Matrix{T}
end
PointSet(coords) = PointSet{eltype(coords),size(coords,1)}(coords)

Meshes.ndims(::Type{PointSet{T,N}}) where {T,N} = N
Meshes.coordtype(::Type{PointSet{T,N}}) where {T,N} = T
Meshes.elements(p::PointSet) = eachcol(p.coords)
Meshes.ecoords!(x::AbstractVector, p::PointSet, e) = (x .= e)

I have a few questions @SimonDanisch regarding next steps:

1. Do you think these 4 function names are ok? Would you like to suggest a better name for the coordinates of the element centroid? I came up with ecoords but I totally understand if you prefer a more verbose name or even center.
2. I've noticed that the construction of regular grids in GeometryBasics.jl is currently limited to 2D. Even though I used a NgonFace{2^N,Int} in the constructor of the RegularGrid{N,T}, the mesh does not compile. Could you please confirm this is the case?
3. I had to add 1 to the dims of the regular grid to construct a grid with dims x dims elements. Would it make sense to adopt a cell-centered interface in the Tesselation call or it is like this because it needs to handle more general cases where the number of elements in each dimension isn't known a priori?
4. There is value in storing the regular grid specs (origin, spacing,...) in the mesh object for later use. Where do you think this RegularGrid mesh type should live? Would it be appropriate to add something like it in GeometryBasics.jl or should I maintain special cases as a separate effort in GeoStatsBase.jl?
5. Could you please share with us which function names from the proposed interface are required by Makie.jl to plot meshes in general? Do you think you could isolate a reduced set of function names from the proposal and implement the plot recipes in terms of these functions only?

The next step after I experiment with the new domain types in GeoStatsBase.jl is to start brainstorming the interface for metadata (spatiotemporal fields over the meshes). This next step will likely influence (and be influenced by) the GeoTables.jl effort that @visr is mentoring.

[juliohm]

I've renamed ecoords to elmcoords and vcoords to vertcoords to be more explicit. Appreciate if you can review the proposal in Meshes.jl and the other questions (2 to 5) above.