Mesh representation (for FEM)

[SimonDanisch]

I wanted to refactor the mesh type long ago in GeometryTypes, and now that I'm working on propper FEM integration into Makie it becomes even more relevant to rethink the mesh representation.

My current thinking is, that at the core one has a set of nodes (vertices) and a vector of cells (faces) that connect those nodes. This can be represented as view(nodes, faces), since the faces do actually represent a view into those nodes. Like this, going from a e.g. tedrahedron to a visualizable surface mesh is just a matter of changing the face type.

So why should this part of the FEM infrastructure live in GeometryTypes and what will be the advantage?
Well, everything graphic related works pretty nicely with GeometryTypes.
I have a decompose infrastructure in place, which supports decomposing arbitrary geometries into its components.
E.g., this makes it very easy to visualize all kind of wireframes, meshes etc, since I can just do:

scatter(x::GeometryPrimitive) = scatter(decompose(Point, x)) # decomposes it into its positions
# extract positions and 2 segment connections between those points
wireframe(x::GeometryPrimitive) = lines(view(decompose(Point, x), decompose(Face{2, Int}, x)))
# creates a mesh with normals for shading
mesh(x::GeometryPrimitive) = mesh(Mesh(decompose(Point, x), decompose(Face{3, Int}, x), decompose(Normal, x))

Node, that decompose(Face{3, Int}, polygonmesh) already works to triangulates meshes with arbitrary faces.
So the Makie recipes could just live in GeometryTypes and work across all meshing/fem packages.

This discussion becomes more relevant, since I just made a big step forward in making Tetgen available as a package, so I expect to use advanced meshing a lot more in the future:
JuliaGeometry/TetGen.jl#7

The problem with the use of view as above is, that node[face] == NTuple{N, Node}, since one face/cell connects multiple nodes.
But in principle, that's most of what we need to represent a mesh. If we want more attributes per node, we can put them into the node type.
I wonder if we should pirate the meaning of view, or just have something like:

struct Mesh{T, N} <: AbstractVector{NTuple{N, T}}
nodes::Vector{T}
cells::Vector{Cell{N}}
end
mesh(nodes, cells) = Mesh(nodes, cells)

We might want to define the node type like this, to contain arbitrary attributes:

struct Node{T <: NamedTuple}
data::T
end
Node((color = RGB(0,0,0), position = Point(0, 1, 2), displacement = 0.0))

The above would work nicely with a trait system!
The beauty is, we won't need to pin down the Node/Vertex type right away, and should just be able to keep using already present structures and slowly converge to one representation.

What does everyon think about this? Ready to move some fem meshing types into GeometryTypes?

Best,
Simon

CC: @KristofferC, @ChrisRackauckas , @louisponet, @Kevin-Mattheus-Moerman, @mohamed82008

[mohamed82008]

I think there needs to be a base mesh data structure without any toppings. This can be defined inside GeometryTypes or preferably in the FEM packages, more on this in the next paragraph. This base data structure will be holding the information that would be stored in a plain mesh file used in FEM applications. For unstructured meshes this would be nodes as a bunch of coordinates and cells as tuples (or the likes) of node indices. For structured meshes, this is a 2D or 3D array of coordinates. The toppings can then be node or cell attributes, which can be numbers, tuples, tensors, or time- and space- dependent callables returning numbers, tuples or tensors.

I think there needs to be a minimal "interface" that any mesh type can define and be automatically visualized. So I recommend defining Abstract versions of types and concrete ones. The visualization functionality can be defined for the Abstract versions assuming the concrete types will respect the defined interfaces. That way if FEM package developers want to define their own Hexahedron type for any reason, they can just subtype GeometryTypes.AbstractHexahedron and visualization will work seamlessly with the new Hexahedron.

The interface for AbstractUnstructuredMesh can be:

1. cells(mesh) which returns some <:AbstractVector{<:GeometryTypes.AbstractCell}, and
2. nodes(mesh) which returns some <:AbstractVector{<:GeometryTypes.AbstractNode}.

The interface for GeometryTypes.AbstractCell can be:

1. cell[local_idx] to get the global node index from the local node index, and
2. The cell must follow VTK's convention for nodes' order.

The GeometryTypes.AbstractNode interface can be node[dim] to get the coordinate along a certain dimension.

Making the interface minimal is important and even more importantly to not expose any visualization internals to the FEM developer to make it more developer-friendly. The interface can be formalized by traits but this may be overkill because it will burden the developer with having to define the traits on their types. I don't see a merit to traits here but correct me if I am wrong.

Attributes can/should probably use traits, since they can return numbers, tuples, tensors, be time- and/or space- dependent, etc. So dispatching correctly on all these "super-types"/"traits" is important when visualizing. GeometryTypes is probably the wrong place for this though. This part can use more ideas.

One point worth highlighting here is that of higher order elements. An <: AbstractQuadraticTetrahedron for instance is probably expected to curve quadratically when visualized. This means that a finer triangulation would have to be made. That's the kind of stuff that will be hidden away from a FEM developer who merely subtypes AbstractQuadraticTetrahedron to define their NewAwesomeQuadraticTetrahedron.

[mohamed82008]

CC-ing @haampie

[ChrisRackauckas]

CC-ing @ysimillides

[PetrKryslUCSD]

Some thoughts on the proposals above:

Some (imprecise) definitions

• Coordinates: an ordered collection of real numbers.
• Node: the basic attribute is Coordinates.
• Connectivity: ordered collection of node numbers (vector, tuple).
• Finite element set: ordered or unordered collection of finite elements of the same type,each element defines Connectivity. The type of the element may be used to dispatch methods specific to the element shape, order, and so on. The constraint is that each element must define connectivity of the same length (the same number of nodes).
• Mesh: ordered or unordered collection of finite element sets, and ordered collection of nodes. The ordered collection of nodes needs to guarantee O(1) access time to an arbitrary node by its number. Therefore the ordered collection should probably be a vector. The constraint is that each node must define Coordinates with the same number of space dimensions.

Discussion

In order to elucidate the potential role that the standardized mesh representation could play, it might be useful to consider how it might be used in a finite element program. What do writers of FE programs expect of the mesh entities?

An incomplete list of typical actions is:

• Find out the manifold dimension of the finite elements set.
• Associate attributes with the finite element set (for instance to associate thickness with a surface).
• Evaluate the basis function values and the gradients of the basis functions at some parametric point. This enables the calculation of the Jacobian matrix, the Jacobian, and the basic interpolation of quantities located at the nodes. This depends on the type of the finite elements that, and therefore these calculations would dispatch on the type.
• Extract the boundary of some finite element set or a subset of the finite element set. This calculation is based on the connectivity, but it needs some understanding of which type of finite element is created by extracting the boundary. In other words there should be some notion of a hierarchy of finite element types based on dimensionality. For instance, the boundary of the tetrahedron with four nodes consists of three node triangles, the boundary of three node triangles consists of two-node line segment elements, and the boundary of a line segment element consists of two "point" finite elements.
• Extract a subset of the nodes based on some criteria.
• Extract the subset of a finite element sets based on some criteria.
• Compute the normal to a hypersurface-like finite elements set.
• Evaluate integral of some function over the mesh.
• Represent results on the mesh, both at the nodes (globally continuous), and elementwise discontinuous.
• Transfer results from one mesh to another (not necessarily compatible).
• Refine mesh by element subdivision.
• Associate some attribute with selected finite elements set.

If the proposed mesh representation only defines the data structure for the mesh, the reason to adopt it in finite element programs is in my opinion not strong enough.

Defining an interface to some basic operations on the mesh (such as those above) might be a sufficiently compelling reason to build the FE program around such mesh representation. Especially if there were reasonable default implementations of this interface that would be easy to layer on top and customize.

[mohamed82008]

Evaluate the basis function values and the gradients of the basis functions at some parametric point. This enables the calculation of the Jacobian matrix, the Jacobian, and the basic interpolation of quantities located at the nodes. This depends on the type of the finite elements that, and therefore these calculations would dispatch on the type.
Extract the boundary of some finite element set or a subset of the finite element set. This calculation is based on the connectivity, but it needs some understanding of which type of finite element is created by extracting the boundary. In other words there should be some notion of a hierarchy of finite element types based on dimensionality. For instance, the boundary of the tetrahedron with four nodes consists of three node triangles, the boundary of three node triangles consists of two-node line segment elements, and the boundary of a line segment element consists of two "point" finite elements.
Extract a subset of the nodes based on some criteria.
Extract the subset of a finite element sets based on some criteria.
Compute the normal to a hypersurface-like finite elements set.
Evaluate integral of some function over the mesh.

Because of the different styles that these can be done in inside the FEM package and the technical FEM knowledge and careful testing they need, I recommend keeping this "standardized mesh" representation to a minimal scope, i.e. for visualization. That is to have a either a set of functions to be defined for your mesh/cell to be visualized by Makie, or to have a basic mesh type visualizable by Makie that can be then made part of a bigger fancier FEMMesh struct in the FEM package with all sorts of data gadgets needed to get the FEM job done.

This de-scoping is as much for the benefit of a visualization developer as it is for a FEM developer so each can do their specialized job without knowing much about the other end of the world, only what they need to know. That's my opinion anyways.

[KristofferC]

I don't think many people will use the mesh struct directly in their FEM code base, but instead use it as a target for converting to when visualizing.

[ysimillides]

I agree with @mohamed82008. While similarities between FEM packages may exist, everyone will have a slightly different focus, and thus structure. From the point of view of Makie/plotting, wouldn't it be better to have a plotting/mesh/trisurf function which takes as an input the node coordinates/connectivity (and values for when plotting solution), and then individual packages can overload this function with their specific FEM structures?

Edit : This seems to be what @KristofferC is also suggesting?

[mohamed82008]

An interesting suggestion which comes to mind from Kristoffer's comment, may or may not be what he meant, is to essentially do the opposite of my first suggestion, that is define the mesh type in GeometryTypes and define the functions needed in the FEM package, e.g. iterators over the cells and nodes of cells, which is really the main mesh-related feature used in JuAFEM for instance.

The following are some of the things that might be asked for by a FEM developer:

1. Give me an iterator over the elements -> e.g. for assembly
2. What is this element's type? -> To know how many geometric shape functions to use
3. For a given element, give me an iterator over the nodes' coordinates -> To find the Jacobian of the geometric shape functions in the element

There may be more things to query, but that's an interesting approach to add to the list.

[mohamed82008]

I can see this last suggestion kind of working for pure Julia FEM packages since one can entertain the fact that the mesh is owned by GeometryTypes and just ask for what one needs, but perhaps not so much for FEniCS where the mesh is owned by Python/C++ and converting all of it to a GeometryTypes compatible form may be too heavy on the memory, since this mesh would still have to be triangulated, i.e. even more memory usage. It would be nice if triangulation can be made as lazy views, that way memory is not replicated, not sure if that's feasible though in the GL world.

[PetrKryslUCSD]

@mohamed82008 , @KristofferC , @ysimillides :
As I outline in the second-to-last paragraph, if all this mesh interface did was to define a structure for the mesh, there is little compelling reason to adopt this structure by programmers as a basis of their libraries/programs.

As @mohamed82008 just now commented, and as my discussion in the last paragraph above suggested, it might be possible to define an abstract interface to such a mesh structure: the operations are really common to all finite element programs, as everybody does it in some way (but their own distinct way). Therefore having an interface might be an acceptable compromise: comply with the interface, but your implementation is up to you. For some basic operations it might be possible to supplement the abstract interface with some sample implementations (dispatched on abstract types).

I don't find such a mesh-type-definition package very appealing for visualization: I agree with @ysimillides that facilities for conversion of volume-, surface-, wire-, and point-like finite elements to displayable entities might be more appropriate for visualization than an actual mesh data structure definition. So then what should be supported is for instance generation of cut and iso (hyper-)surfaces, subdivision of finite element shapes for the purpose of refining the geometry (for instance to support NURBS and higher-order polynomial finite elements), and so on.

[mohamed82008]

I think the main questions to answer here are:

Should the mesh type be owned by GeometryTypes or the FEM package?

1. If owned by GeometryTypes, what is structure of this (these) mesh type(s) and what is the interface that a FEM package may use from such mesh types, i.e. iterators and adding attributes and such?
2. If owned by the FEM package, what is the minimal interface/trait-group required to visualize this mesh in Makie?

I am not sure which path @PetrKryslUCSD or @KristofferC are suggesting.

[mohamed82008]

The main compelling reason I can see for adopting even a basic mesh type and building stuff on top of it, is if it provides nice utility functions, e.g. iterators, and it can be visualized. Also it would be nice if it can be GPU-stored or distributed on multiple machines, that way one can ask for the part of the mesh on this machine, or apply a series of GPU-map function on subsets of non-interacting elements, etc. I think this last option assumes we went with the first path in my previous comment, since then Simon can define the mesh structure himself, including all the GPU mapping and coloring stuff.

[KristofferC]

The most efficient storage of the mesh might depend a lot on your application. For example, if one does hierarchial refinements of quadrilaterals (and then take care of the hanging nodes) one probably want some multilevel format, where each refinement level of the elements are stored contiguously. While for other types of adaptivity, other storage formats can likely be better.

My point is that it feels quite unlikely that there can even be a common "base" format that would be used directly in the fem code. But a common format to act as an entry point to nice visualization (similar as to how you now would export to VTK) is likely more approachable.

[mohamed82008]

With the exception of mesh refinement which I have no experience in, VTK itself supports multiple base mesh types:

1. Unstructured mesh (can be homogeneous or heterogeneous),
2. Structured mesh (quads and hexas),
3. Rectilinear (non-uniform pixels or voxels), and
4. Uniform rectilinear (uniform pixels or voxels)

Being able to exploit the memory efficiency and convenience of structured or rectilinear meshes, including the cheap boundary and neighborhood information, is important when iterating over the elements or applying boundary conditions. I may be wrong, but I think any adaptively refined mesh will have to be one layer of abstraction above these base types, so there is nothing to lose by starting somewhere and moving on from there.

Yes the main focus here is visualization, I assume this is why we are even talking about this. While the above list does not cover all use cases for mesh types, I think the heterogeneous unstructured mesh type is a pretty convergent type, because of its flexible definition. Essentially nothing stops you from throwing in there heterogeneous, overlapping elements of different dimensions, e.g. some linear lines, quadratic triangles, vertices and cubic hexahedrons.

So all the unsupported multi-level specialized mesh types can still be visualized so long as they are converted to a HetergenousUnstructuredMesh struct and passed to Makie. So I guess the first step is to support visualizing some useful mesh types, then things can naturally move on from there. Whether or not these types are attractive enough to be used inside the FEM package for more than just visualization depends on the expectations the FEM package/developer has from the mesh type. If these expectations can be provided for in GeometryTypes, great, then it can be adopted, and I think it can. If not, then at least we have a working FEM mesh visualization tool.

[SimonDanisch]

@ysimillides

and then individual packages can overload this function with their specific FEM structures?

I created this issue specifically to discuss how to make this overloading easier ;) Because, if we have a clever interface in GeometryTypes, one possibly doesn't need to overload anything in the FEM packages.
And if we do it right, we'd just need to move a very small part of the infrastructure to GeometryTypes.

There are things that are purely geometric, that don't need to be reinvented in all FEM packages.
E.g. triangulating a cell, getting all geometric nodes of a mesh etc.
This functionality would fit very well into GeometryTypes.

Another reason I opened this issue is, that now with Tetgen, which I consider as a purely geometric package not specific to FEM, I also need a more refined mesh type in GeometryTypes with nodes + cells and holes etc ;)
So I will have a more advanced mesh type in GeometryTypes, which could also serve as a basic mesh type in any FEM package.

In the end, I probably just need those interface functions defined:

decompose(Point, mesh)::Vector{Point} # the nodes
decompoe(Face{2, Int}, mesh)::Vector{Face{2, Int} # line segments for wireframe view
decompose(Face{3, Int}, mesh)::Vector{Face{3, Int} # triangulated view for surfaces

Since tetrahedrons etc are well defined geometric structures, we could have those decompositions into Face{2}/Face{3} implemented in GeometryTypes.
All those Cell types here already remind me of the face type in GeometryTypes:
https://github.com/KristofferC/JuAFEM.jl/blob/master/src/Grid/grid.jl#L24
So we might as well merge those concepts, put a simple, generic mesh type on top of that, and ultimately make plotting + tetgen integration easier.

[mohamed82008]

I think some basic PR that supports Makie would be a good starting point. Then it becomes obvious how to extend it to different types of meshes and elements. I want to contribute but I need a backbone program first.

[Kevin-Mattheus-Moerman]

This is an interesting discussion. I'm not a Julia user yet (currently use MATLAB and GIBBON) so I do not fully understand all that is suggested here.

I just wanted to add that there are many different types of element descriptions. For instance there are 4, 10, and 15 node tetrahedral elements, and 8, 20, 27 node hexahedral elements. The basic 4-node tetrahedron has 3-node triangular faces while the 10-node tetrahedron has 6 node triangular faces. The connectivity defining the element varies from one implementation to another. Therefore also the way one composes faces from the element description.

Perhaps you can demand that users convert their element descriptions to face descriptions themselves. This way you avoid issues with differences in element descriptions. For a tet4 element one would create a list of 3-node triangles, for a tet10 description a list of 6-node triangles. In GIBBON I have conversion functions e.g. element2patch which converts a list of elements and element data (e.g. stresses) to a list of patches (faces) and equivalent face data. I then visualize this with a generic mesh viewer (patch function) which uses the triplet: faces,vertices,colordata, colordata can be defined for faces or for nodes (if scalar its colormap driven, if nx3 then it is assumed to be RGB data).

[PetrKryslUCSD]

@Kevin-Mattheus-Moerman : The "streaming" idea is attractive. It needs some generalization though. It would be good to stream for instance volume elements with the purpose of extracting an iso surface, or to stream surface elements with the purpose of extracting the boundary.