Finite element mesh interface

Dear all,

Do you know a paper, book or any published material that attempts to define an interface for general finite element meshes?

I am in the process of finalizing a proposal interface, and would like to share it with the community for feedback. If you have experience in this area, I appreciate if you can share material and comments.

The interface should be similar to what Tables.jl provides for tabular data.

3 Likes

I’ve built an interface based on tables.jl in GeometryBasics.jl, which I hope to have working with fem… But the integration process got a bit stuck;)

2 Likes

You might wish to check out for instance https://www.scorec.rpi.edu/pumi/. But it depends very much on what your goals are: interface to the mesh can be trivial or unbelievably complex or anything in between, depending on the objectives (what should be supported and where and how, …).

3 Likes

Thank you @sdanisch, I will definitely take a closer look. :100: From quickly scanning the code, it looks like you have a nice implementation of mesh elements and a mesh type with metadata on vertices and faces. I am brainstorming the interface so that we could plug in any mesh element with different properties. For example, we could consider using Polyhedra.jl mesh elements and solve optimization problems directly on meshes, like structural optimization in civil engineering, or mesh refinement with the optimization frameworks we all love :slight_smile:

Yeah, so the idea is, that GeometryBasics offers the interface for such use cases :wink: I guess I haven’t done much yet to document what’s there

2 Likes

Thank you @PetrKryslUCSD, I fully agree it depends on the goals. The initial goal is a very simple interface to enable efficient geostatistical estimation/simulation of properties over general mesh types with GeoStats.jl. Ideally, we could simply perform geostatistical estimation/simulation of properties, pass the resulting mesh with estimated properties to FEM solvers, and plot the results of a physical simulation directly with Makie.jl.

I will read the source carefully, but I couldn’t find for example an interface to query the volume of different elements, the surface_area of the element, the baricenter. Also, it would be nice to have very specialized mesh types for regular grids that don’t take any storage other than the origin, spacing, and dims tuples. That way we could treat general Julia arrays (or images) as regular grid meshes.

1 Like

Feel free to open issues/prs etc! Would be nice to evolve GeometryBasics to fit the case!

2 Likes

Currently I am working on creating mesh interface based on DirectSum.jl and Grassmann.jl instead of using GeometryTypes. This new mesh interface will be based directly on simplicial complex of DirectSum and Grassmann packages instead of relying on GeometryTypes, since these foundations greatly expand the capabilities of the simplex algebra. Later today I will be releasing v0.5 of both those packages, to help facilitate the mesh design I am planning.

1 Like

Thank you @chakravala, could you please elaborate on the advantages of the theory in practical applications? Also, do you have an example of this interface in action that is easy enough to grasp?

Does this mean that the mesh types you are considering are simplicial in the sense that elements cannot assume more general forms like quadrangles, hexahedrons, etc?

1 Like

Yes & no. Indeed, there are some practical advantages to the SubManifold representation of a simplicial mesh, which makes the approach better suited to a simplicial complex instead of the quadrangles and hexahedrons. It is still possible to work with more general forms like quadrangles and hexahedrons, but such forms would have an extra layer of complexity in its sub-chain represenation.

For FEM solvers, we need more general element types with lower approximation error. So it would be nice to interface with the mesh types you have in mind, but also with other types of meshes. This is why we need an abstract interface to talk the same language. :+1:

Currently, the DirectSum package provides a SubManifold{V,N,B} representation of a simplex. This versatile type is extended by Grassmann algebra so that it can be combined into point clouds and chain complexes, (including quadrangles which can be obtained by combining oriented simplices).

The fundamental building blocks of a quadrangle can be represented in terms of a simplicial chain.

I would certainly like to hear more about what people think is a better approach for unifying the simplicial building blocks with the more general elements like quadrangles. The DirectSum package provides the foundation for any Simplex based on a SubManifold representation, while the Grassmann package provides an extension to combine SubManifold into Chain and MultiVector simplicial complex.

Please do share any ideas about interfacing if you have any suggestions.

In CompScienceMeshes, I have interfaces for:

  • Meshes: essentially iterable containers of cells. skeleton and connectivity give lower dimensional entities and how they fit together.
  • Charts: from the topological constituents of mesh, you can build geometric cells by e.g. chart(mesh, cell)… For now only simplices (in any dimension) are implemented. For some dimensions I have implementations for the most common geometric predicates (overlap, isinclosure). volume(simplex) also works.
  • Neighborhoods: Given a chart and a tuple of parameter values, you can get the value under the chart at the parameter. The structure is enriched with tangents, jacobians, and normals as required in FEM/BEM assembly.
  • Quadpoints: for the simplical chart and the spherical chart, quadpoints computes through pull-backs onto the domain quadrules for these domains. I.e. sum(w*f(p) for (p,w) in quadpoints(chrt)) gives the integral.

Quite a lot of common and not so common finite element spaces built on top of this can be found in BEAST.jl. For most combinations ‘trace(X)’, ‘divergence(X)’ etc. just work.

5 Likes

That is great @krcools :100: the interface in the docs of CompScienceMeshes.jl is very aligned with what I had in mind originally. Let’s join efforts to put a lightweight package together with a good set of traits for general meshes. We started the discussion here: https://github.com/JuliaGeometry/GeometryBasics.jl/issues/15

Such lightweight (interface) package could be imported by all packages implementing meshes and so we could all communicate across different frameworks. A first draft is available in that issue and here as a repository, please let me know if the proposed interface serves to build connectivity matrices efficiently for FEM solvers. Feedback is very welcome.

Indeed, my packages also have the skeleton functionality and so on. After thinking some, I think I figured out a new type object to add to my packages:

I intend to add a feature which generalizes the notion of point cloud to higher dimensional simplices. Specifically, CellChain{M,n} <: Manifold{n} object, which will be able to represent an iterable set of sub-Chain elements from a Manifold, kind of like a point cloud generalized to chains of Grassmann multilinear Simplex.

This single new type would provide sufficient new type information to build iterables over sub-chains, which is how a more general cell-chain would be represented.

This package is very nicely done. However it appears to only support simplex meshes, of the first order.
Has any thought been given to extending this arbitrary finite element shapes, of arbitrary order?

To be honest not much has been down beyond non-curved simplices. A co-worker has been successful in using the framework for subdivision surfaces. The geometric predicates - needed in methods where trial and test functions live on non-matching grids - become a nightmare in the general case. The rest of the interface is more easily adopted to this case.

Can this be handled by integrating/applying a k-form ( k=2 for surface, 3 for volume ) over the simplex? Using Discrete Exterior Calculus or something similar? Then all the higher order information is contained in the use of k-forms.

In the current implementation in BEAST.jl, k-forms are evaluated by first computing them on the flat reference simplex and then using an appropriate push-forward formula. So yes, any higher order geometry will be automatically inherited from the definition of the curved simplex as a Chart (in the CompScienceMeshes meaning of that concept).

BTW, also some of the integration routines I use, such as SauterSchwabQuadrature.jl, are defined in terms on the parameter space. So also these will keep working for curved geometries.

On the other hand, increasing the maximal degree of the finite element (as opposed to the degree of the parametric representation of the patches) does require a certain amount of work.