Types of adjacencies

[PetrKryslUCSD]

Hi,

Congratulations on an impressive piece of work!

I started studying the package, but I don't see so far if it is possible to derive adjacency relationships such as faces to nodes, nodes to faces, faces to faces, etc. And what the definitions should be: for instance face to face. One to many?

Thanks.

P

[j-fu]

Hi, thanks! The package has some pre-history from our work in C++. And it would be significantly less comprehensive without the contributions of @chmerdon who does research in finite elements (See GradientRobustMultiphysics.jl).

That said, the package so far is being developed in order to get our "frontend packages" GradientRobustMultiphysics.jl and VoronoiFVM.jl to speed (which means that we can do simulaition projects and publish papers; without that it would be hard to justify the efforts we put into Julia), and I am well aware that it is quite eclectic and could be better documented...

So answering your question might help to improve this situation.

Some core decisions behind this package:

• Store points, adjacencies etc. as matrices instead of vectors of points or tuple like structures (this allowed to carry over our code from C++ rather quickly)
• Lazy creation of adjacencies: if you create a grid g, you just call g[FaceNodes] and the adjacency is created on the fly if it was not there yet.
• Storing grid data in a dict with types as keys. The idea was to be able to dispatch getindex etc on types. Not 100% sure now if this is really clever but it works ok for sure. One point is here that with storing element type information in an array of types we ran into dynamic dispatch of the corresponding getindex.

The first item sets up a dividing line to the style of other Julia packages. One of the ideas in the back of my head is to provide interfaces which would allow to bridge this divide. reinterpret might allow this quite conveniently.

And I am also aware of FinEtools.jl which seems to have grown quite impressively as well since I first looked.

[PetrKryslUCSD]

The reason I started looking into ExtendableGrids was because I was attempting to understand GradientRobustMultiphysics (@chmerdon). :-) And I couldn't quite understand how the assembly worked, which meant that I had to understand extendable grids.

The packed vector a la sparse matrix is a good idea. I was contemplating to introduce that in https://github.com/PetrKryslUCSD/MeshCore.jl for incidence relations that were "irregular". In the end I decided to use a vector of vectors, because of convenience. But the packed vector would be a much better choice because of the reduced memory fragmentation.

An interesting aspect is: when the collection of the entities is heterogenous (aka different element types), one has to store the type (which you do), or dispatch with a switch statement. In either case it seems to me that during assembly, a context switch is needed. For performance reasons I tried to avoid this both in https://github.com/PetrKryslUCSD/FinEtools.jl and in https://github.com/PetrKryslUCSD/Elfel.jl. Then, when looping over finite elements, the quadrature rule is fixed, and a number of things can be precomputed once before the loop starts. I'm trying to understand what happens when the collection of finite elements is heterogenous: how is the context switching handled, and what the performance impact might be. I have had some ideas on whether and how to extend MeshCore to heterogenous collections. The advantage of MeshCore is that it can compute arbitrary incidence relations needed in more complicated FEM (for instance when degrees of freedom live on edges). Its disadvantage is currently that it can only work with homogeneous collections. This is the consequence of not wanting to switch contexts when working with the collection.

The lazy creation of the adjacencies (or, as I call them, incidence relations) is a great idea. I couldn't quite see how the creation of arbitrary adjacencies (for instance face to edge) was done. I must have missed something.

[chmerdon]

Thanks for your interest in GradientRobustMultiPhysics :) I try to give some insights into some aspects of the assembly:

1. How the assembly works in principle can be seen in one of the low level API examples, for example in the Pluto notebook for Navier-Stokes: https://chmerdon.github.io/GradientRobustMultiPhysics.jl/stable/nbhtml/LowLevelNavierStokes/
   Here you can see how the assembly loops work and what low level objects are generated prior the loop (like the QuadratureRule or FEEvaluators that evaluate an operator of an finite element in the respective quadrature points,
   and already takes care of different FE continuities (H1, Hdiv, Hcurl) and orientation etc.). Finite element spaces are encoded by FESpaces and also have lazy evaluators for CellDofs, BFaceDofs, FaceDofs etc., so these arrays are generated when first requested. This low level example shows that the assembly loop can be, in principle, allocation-free. Unfortunately this is not the case, yet, for the high level API.

2. The high level API allows to generate more or less arbitrary LinearForms, BilinearForms, NonlinearForms (called assembly patterns) and arrange them inside a PDEDescription (together with boundary data and global constraints). Then, these are assembled automatically when the solve is called. This assembly also takes into account different ElementGeometries (with separate pre-allocated quadrature rules etc.). As you described this requires some additional checks and dispatches in the assembly loop. If you want to dig deeper you may want to start in https://github.com/chmerdon/GradientRobustMultiPhysics.jl/blob/master/src/assemblypatterns/bilinearform.jl to see the loop for a general bilinearform. The assembly loop starts in line 226, but the geometry search is hidden in the update_assembly! call in line 232, and unfortunately the code is not very readable from there ^^;; Also, I don't really use this mixed-geometry feature and not many finite elements are implemented that support non-simplices. If I find time I want to rewrite the assembly from scratch as there are also several small allocations happening that slow down the code for large problems.

3. The creation of FaceEdges (like CellFaces, BFaceFaces, FaceNodes etc.) is implemented in https://github.com/j-fu/ExtendableGrids.jl/blob/master/src/derived.jl. One can also transpose them to e.g. get NodeCells or EdgeFaces (never needed that one, though, I think) with the atranspose function of ExtendableGrids. Faces to faces does not exist (only BFaceFaces to get the face number for the bfaces).

[PetrKryslUCSD]

Cool, thanks for the guide, Christian. I will follow those breadcrumbs...