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.