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.

[visr]

Discourse thread for reference: https://discourse.julialang.org/t/finite-element-mesh-interface/33325

Just want to mention mesh types are becoming more frequent in geospatial context as well, so I'm quite interested in having an interface that could help there as well.

See for instance MDAL that @evetion make a start at wrapping in MDAL.jl.

Mesh data are (currently at least) not part of GeoInterfaceRFC.jl, but it is relevant since it is an interface that we hope GeometryBasics will subscribe to, to allow interoperability with other geospatial packages. Since GeometryBasics is close to many of the ideas in the interface, the implementation is very simple:
JuliaGeometry/GeometryBasics.jl@master...visr:rfc.

[juliohm]

That is a quite important use case as well @visr , thanks for sharing. We can start by creating a minimum interface for discretising geometrical objects into mesh elements as above, and then start adding more traits for projections, geometry, etc to include the use case with geo maps. Hopefully we can find a nice set of traits that gives us the power to express both 3D models of the subsurface as well as 2D projections of geographical objects in spheroids. That would be amazing.

What I will try to do next is take the mesh type defined here in GeometryBasics.jl and make it adhere to the interface. This will consist of simply making sure that access functions for the elements and coordinates of elements are defined for each type of element.

Will keep you guys posted of any updates.

[SimonDanisch]

Some of these are already implemented with a different name... Btw, I'd much prefer writing out the names, e.g. a novice user likely (or actually me, I can only guess) won't know what ctype or cbuff does ;)

[juliohm]

Thank you @SimonDanisch perhaps a doc string solves the issue? Please let me know if it is not clear and we can rename to something more explicit.

I've put a draft here organized by different types of traits, inspired by the existing efforts: https://github.com/juliohm/Meshes.jl

I will try to make the types in GeometryTypes.jl and CompScienceMeshes.jl adhere to the interface as a proof of concept. Please let me know if the repo should be migrated to JuliaGeometry, it will certainly be a better home than my personal account.

I've separated the traits into "mesh traits" containing general information about the mesh object/type, "element traits" containing information about mesh elements, and "vertex traits" containing a simple access function for the coordinates. The source is organized as follows:

src/
--core.jl ---> core API that all mesh types can and should implement
--basic.jl ---> basic API derived from core API (i.e. fallback implementations)
--geotraits.jl ---> geometric properties of the mesh and its elements (e.g. centroid ,volume)

We can continue extending the interface in separate files to fit all use cases that were raised. @visr @evetion do you have suggestions about projections and geographical coordinates? How do you see them interacting with the current verbs defined in the draft?

[SimonDanisch]

Thank you @SimonDanisch perhaps a doc string solves the issue? Please let me know if it is not clear and we can rename to something more explicit.

I'm roughly sticking to the YASGuide, and I also find it super useful for reading code by someone else to have more explicit names without starting to look up docstrings (e.g. when reading code on github, you can't even look them up... And even if you have a running julia session, it's not always clear where the function comes from and therefore unclear how to look up the doc string ;) )

[visr]

@visr @evetion do you have suggestions about projections and geographical coordinates? How do you see them interacting with the current verbs defined in the draft?

Not at this point, and I'm not sure these should be part of the mesh interface. As long as there is a solid, well though out mesh interface that supports a wide variety of mesh types, I'm sure it will be useful.

[PetrKryslUCSD]

For the mesh interface design, wouldn't it be helpful to start with a verbal description of what is to be expected of such an interface. What types of operations should be supported, roughly at what level, and so on. I think writing code at this point is a bit premature.

[juliohm]

@PetrKryslUCSD I have some ideas of what is expected from the interface coming from a statistical perspective, but certainly will miss other views from FEM solver writers and visualisation package writers. I want to get this draft going so that we at least have a proof of concept that this interface includes two existing implementations designed with these other goals in mind.

[PetrKryslUCSD]

Perhaps this may help? https://petrkryslucsd.github.io/FinEtools.jl/latest/guide/guide.html#Finite-element-1

[juliohm]

Awesome :) one more reference to dive in.

…

On Wed, Jan 15, 2020, 22:08 Petr Krysl ***@***.***> wrote: Perhaps this may help? https://petrkryslucsd.github.io/FinEtools.jl/latest/guide/guide.html#Finite-element-1 — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <https://github.com/JuliaGeometry/GeometryBasics.jl/issues/15?email_source=notifications&email_token=AAZQW3PXWSCHYEDXEOSJCOLQ56XP3A5CNFSM4KG4PJS2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJCMMAY#issuecomment-574932483>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZQW3JNLCW4Q6MVCZXOGSDQ56XP3ANCNFSM4KG4PJSQ> .

[PetrKryslUCSD]

This is interesting: https://codedocs.xyz/CEED/FMS/. I found it especially so since it uses several concepts included in my FinEtools. :)

[PetrKryslUCSD]

I wonder if this conversation could be revived? Again, I am asking: what is the purpose of the interface?

To answer my own question I will try to describe a scenario.

Let us say there are application libraries A, B, C (as you say, perhaps statistical modeling, finite element modeling, visualization). Then there is a package the user writes to define the mesh, D. The package D defines the structure for the mesh, functions to access this data, and the data for some particular meshes. The mesh data then should be available for processing in the libraries A, B, C. This can be done provided these libraries (A, B, C) are all aware of the interface: in other words, the interface must be defined in its own package, I, and the libraries A, B, C be using I.

So in this way the packages A, B, C are insulated from the actual implementation and storage of the mesh data in package D. However, they all need to be aware of I. And so do all the packages of the type of D.

Is that in any way consistent with what this conversation is about in the minds of those who already participated?

[juliohm]

That is pretty much it @PetrKryslUCSD. The ultimate goal is to be able to seamlessly transfer mesh types across different frameworks and have them work out of the box, no matter the internal storage layout, the different complexities involved, etc.

Sorry I didn't get back yet regarding this interface initiative. I've migrated the repository to the organization already, but got busy with a submission deadline for unrelated research. The issue is on my radar, and it will certainly be addressed to some extent before my next paper submission, which depends on it.

It would be really nice if we could prove the concept with the GeoStats.jl + FEM solver + Makie.jl combo. Will keep you posted of any progress or if I have questions regarding the design of the traits. My plan, as I shared above, is to first refine the interface to accommodate the existing mesh types in the ecosystem, starting with the mesh types defined in GeometryBasics.jl and CompScienceMeshes.jl.

[PetrKryslUCSD]

OK. So I take it your plan is to work on I? It seems to me that in order for this to make sense, package I should not reflect the realities of current mesh implementations (GeometryBasics.jl and CompScienceMeshes.jl). To begin with, they are fairly limited in what they can support (mostly only simplexes?).

I believe that a good place to start would be to list the functionality that this interface should support. Evaluation of interior and boundary integrals, location of points, evaluation of various fields (and differential operators on fields), topological operations (merging, splitting, boundary taking), and so on.

[PetrKryslUCSD]

I guess the alternative to an interface package with some functionality (which should work without any need for access to implementations other than through the interface) would be a bare-bones access to plane-vanilla mesh data (coordinates, connectivity), leaving everything up to the packages using the interface.

[juliohm]

Yes, the plan is to work on the set of traits in "I". You can find a very rough initial draft from that discussion in the beginning of the month here: https://github.com/JuliaGeometry/Meshes.jl

I certainly don't have the expertise to cover the FEM requirements well, and your help is much appreciated. To give you some context, my experience with FEM is limited to writing a simple 3D heat transfer solver with standard Galerkin methods for a specific PDE in C++.

What I would like to do before we move on with more advanced use cases for this interface, is to get the core API right (the coordinates and connectivity that you mentioned). My understanding is that the GMSH and VTK specifications, cover most mesh types of interest. After we have this core API set, we can add iterator traits for faces, edges, etc, and possibly default implementations based on the core traits.

[PetrKryslUCSD]

I think traits is the right tool. In addition to the mesh itself, I think also the "element" types should be equipped with trait dispatch. It other words, imagine you have an FE with four nodes. Obviously it must be handled differently when it is a quad as opposed to a tet.

[PetrKryslUCSD]

By the way, I don't think it is enough to consider homogeneous meshes. Mixing triangles with quadrilaterals, or wedges with hexahedra is quite common.

Also, there are meshes for the interior and meshes for the boundary. Obviously they are related by compatibility requirements. Should these meshes be stored separately or together?

Quite often multi-material domains require some treatment to distinguish between the materials. Separate meshes, or perhaps labels.

One other point: Incompatible meshes with tie constraints are also a possibility. Again, this may require storing meshes composed of different element types. So either separate meshes, or meshes that are heterogeneous.

[juliohm]

Fully agree. We have to provide a set of traits for elements of the mesh, and that is in the plans. The traits should be such that the different branches of code are determined at compile time without runtime penalties.

What taxonomy of elements you currently have in mind? The draft could start by listing the different element trait types like:

abstract type Element end
struct Simplicial <: Element end
struct Quadlateral <: Element end
...

This taxonomy will interplay with the connectivity traits and vertex ordering assumptions.

[PetrKryslUCSD]

FinEtools separates FEs by manifold dimension. Then by individual characteristics (how to compute basis functions, how to determine boundary, how to determine point inside/outside, etc.).

[SimonDanisch]

Btw, the mesh type in here already supports arbitrary connections & vertex types, and it also supports arbitrary edge/face/vertex metadata!

[PetrKryslUCSD]

Here is a sample implementation of the mesh interface and a package that wraps FinEtools meshes.
Only a little bit of the functionality is actually implemented.
https://github.com/PetrKryslUCSD/TestMeshInterface.jl.git
https://github.com/PetrKryslUCSD/MeshInterface.jl.git

[juliohm]

Thank you @PetrKryslUCSD , I am collecting all the resources shared so far to tackle this interface seriously in the near future. When time allows, I will try to list the major use cases behind the interface.

[Paulms]

Hi, I have been learning recently the virtual element method VEM (works in polytopal heterogeneous meshes) and replicated some of my class examples by modifying the JuaFEM.jl code. The base mesh type that I found useful for VEM is the following:
https://github.com/Paulms/jFEMTools.jl/blob/master/src/mesh.jl
I hope you find it useful as an additional perspective for the interface.

[evetion]

@visr @evetion do you have suggestions about projections and geographical coordinates? How do you see them interacting with the current verbs defined in the draft?

For me, adding a CRS attribute would be enough to have this working, you could refer to the CoordinateReferenceSystemFormat type here: https://github.com/JuliaGeo/GeoFormatTypes.jl/blob/master/src/GeoFormatTypes.jl

Mesh data are (currently at least) not part of GeoInterfaceRFC.jl, but it is relevant since it is an interface that we hope GeometryBasics will subscribe to, to allow interoperability with other geospatial packages.

To be fair, Simple Features does specify some mesh support (PolyhedralSurface, such as TINs), see my issue here: JuliaGeo/GeoInterfaceRFC.jl#3. I would like to support the Meshes implemented here.

Lastly, there are proposals (disclosure, partly made by my organization @Deltares), to support meshes (output from our unstructured numerical models) in structured formats such as NetCDFs, called UGRID. The API is actually quite relevant and written out completely: https://ugrid-conventions.github.io/ugrid-conventions/.

[PetrKryslUCSD]

I have designed a lightweight mesh library. I believe it could accommodate the desirable features listed above. Even though it currently implements only linear shapes (simplex and cuboid types), adding quadratic and higher-order shapes is trivial. It may even be possible to accommodate shapes with implicit geometry (such as for topologically regular grids). https://github.com/PetrKryslUCSD/MeshCore.jl

[juliohm]

Thank you @PetrKryslUCSD , it seems like you have an implementation of meshes for FEM, correct? I see some structs defined there, and some actual implementations of how things behave.

I am almost done with that journal paper, and will start reviewing the prior art soon. Thank you for sharing yet another resource. 💯

[SimonDanisch]

Cool! There is lots of overlap with GeometryBasics, how do we want to handle that?

[juliohm]

I maintain that viewpoint that we can work together on a lightweight Meshes.jl interface (without implementations), and let people implement the interface for their own mesh types.

P.S.: I bought this book on Voronoi tessellation that arrived yesterday, and I hope that I will find the time to start working on the interface by the end of the month. Do you have suggestions of books on meshes and meshing algorithms in general?

[PetrKryslUCSD]

There are good papers I'd recommend. I will make up a list and post it here.

There are also poor mesh topology designs: The one underlying FENICS (as implemented in Dolfin) is for instance wrong.

[mohamed82008]

Yes so one can have arguments specifying which mappings get pre-computed and stored at mesh construction time. Dispatch should then be able to do the 0 -> 3 query for example in O(1) time if possible. Otherwise, an O(N) fallback can be used. If the query is done once or a few times, the user may not need the additional caching. But if an algorithm will call 0 -> 3 many (~N) times, then the user can opt in to the additional caching to avoid a quadratic complexity.

[mohamed82008]

I think this automatic caching and iterator interface plus visualization are enough reasons for me to switch to the new mesh interface in TopOpt.jl. It will simplify a fair bit of code that I have over there.

[PetrKryslUCSD]

Yes, I see the value. At this point I am redesigning the lib using incidence relations consistently. It should be checked in later today. That will fill in all those incidence relations in the "guide" table that make sense (not the d -> d, d > 0; those need definitions through a third entity. That is another place where the Logg mesh interface of FENICS does not make sense.)

[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:
JuliaGeometry/GeometryBasics.jl#62
JuliaPlots/AbstractPlotting.jl#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.