Treat bases as additional trailing axes

[adtzlr]

Hi @kinnala !

Have you ever considered to treat the data of the bases (e.g. the gradient data) simply as additional trailing axes in addition to the cells and quadrature-points axes? Recently, I tried to outline the assembly in minimal scripts and this seemed like an obvious choice for me (at least for an isoparametric hexahedron).

The snippet is undocumented (yet), but I'm sure you get the idea in SolidBody() immediately.

3D Hyperelastic Finite Element Analysis in 150 Lines of Python Code

(dhdX are the basis function gradients and F is the deformation gradient with δF (or grad(v)) and ΔF (or grad(u) in scikit-fem notation).

Thank you for your constant improvements to scikit-fem!

[kinnala]

I think I might've tested something similar and eventually settled to the current approach in scikit-fem. However, numpy has improved a lot and most likely the conclusions I have made from the performance are not valid anymore.

As it might interest you as a FEM developer, I have lately done some experiments by dynamically compiling FEM assembly routines using gcc/clang from Python. In particular, I have written some multithreaded assembly routine templates in C and the relevant parts of the routine are replaced with basis function expressions and weak form expressions by some simple Python code. The resulting C code gets written into /tmp/foo.c, compiled using gcc/clang and then imported into Python dynamically. This is all done within Python by a simple call such as

form = BilinearForm("ux * vx + uy * vy + uz * vz", ElemTetP1())  # the string is written directly into the C source code template
form.compile()  # this will call gcc/clang and afterwards form.assemble() works as usual

The template C source code is programmed carefully to not depend on any C level packages/libraries so the only true dependency is a standard C compiler such as gcc/clang. It uses pthreads.h for multithreading which is also available on Apple laptops.

There are at least three reasons for these experiments. For instance,

1. implementing complicated forms such as those in hyperelasticity will be slow at Python level because of lots of back and forth between Python and numpy low level code
2. high order basis functions that require lots of integration points will be slow using the current approach in scikit-fem (a different approach would be required for h-FEM and p-FEM)
3. it is challenging to do proper multithreading in Python

Even in some simple cases (e.g., 3D tetrahedral MINI element) I have gotten 100x improvement in assembly time. I currently cannot foresee this work becoming part of scikit-fem due to backward incompabilities. It would most likely be another package which is very close to scikit-fem in its API and even more focused to assembly only.

[kinnala]

Thanks to some "recent" developments in sympy, it is possible to create a symbolic language for finite element weak forms pretty easily. This was not possible when I started scikit-fem, but now I can implement the following in very few lines:

u, v = define(u=0, v=0)  # these are sympy expressions, 0 means scalar
form = dot(grad(u), grad(v))  # dot, grad, etc. are defined by me
form.to_c()  # this will output the string "ux * vx + uy * vy + uz * vz"

[adtzlr]

Hi @kinnala,

sorry for not replying but Influenza caught me a few weeks ago and I totally forgot to answer when I was healthy again. Thank you very much for your detailed insights, very inspiring!

I could've imagined that you already tried something like that.

Unfortunately my coding knowledge is very limited beside Python/Fortran/Matlab, but the C-compilation sounds promising. Casadi does something similar if I remember correctly.

Also thanks for the hint on Sympy. Will have a look at it.

In FElupe, the assembly is more or less backend agnostic. It is implemented in two stages, a) by a constitutive material formulation with the gradients/hessians of the hyperelastic potential (strain energy function). And b) the integration+assembly of the gradient/hessian vectors/matrices. But by doing so, the flexibility of defining forms in a more general way like scikit-fem does is lost. This is not so important if one sticks to (hyperelastic) solid bodies like me but a more important concern for a more general library.