I am writing a simple differential equation solver using the finite element method. While I have done this in the past, I often find that there is a big chasm between the mathematical form and the implemented code – specifically in the ‘‘assembly’’ of the linear system of equations. I wind up with a piece of code for the Laplace problem, for linear elasticity, etc. Of course, large sections of code are reused, but I am looking now to bridge that ‘‘last mile’’ of connectivity. A large part of the inspiration behind me pursuing this is due to the fenics project (fenicsproject.org).
I can illustrate this with an example. The math is simple enough, but it is not really essential to follow it in detail. In the parlance of the finite element method, say I am looking for a function u that satisfies an integral equation for all ‘‘test’’ functions v.
If we approximate u and v as a linear combination of a finite function basis, we can ask for the coefficients of this basis such as to satisfy the above equation.
So the equation above specifies a sort of ‘‘recipe’’ to assemble the linear system of equations.
Going from math to code
I may be wrong, but I feel like this is a job that can be achieved via metaprogramming. I have defined types for u and v (in the parlance of the finite element method,
TestFunction respectively.) I have defined the action of operators like \nabla in terms of how they manipulate my types.
I also have the basic infrastructure of solving such a problem in place. i.e. given a mesh discretization, I can compute all the things necessary like my basis functions and their derivatives, gaussian quadrature etc. I know the exact form of a Julia expression that I can evaluate on one ‘‘element’’ of my mesh to obtain its contribution to the linear system.
Now I want to write a (hopefully?) simple macro that can take my math expression and turn it into a corresponding Julia expression.
Ideally, I would like the macro to behave something like this for the example above:
would give me
K[I,J] = basis[I][k]*basis[J][k]
(the repeated indices in the
@assemble statement just tell me that the operation is a ‘dot-product’ and is similar to the
@tensor macro from
Now, after all this long-winded chit-chat here are some more specific questions:
- Is a macro the way to go for something like this? If yes…
- Linearity: Is there a way that I can convert
:(u(v + w))into
:(uv + uw)? This will allow me to write more general equations as input without having to simplify them down by hand first. From the mathematical side, this is also essential, so it would also act as a good check that the input is valid. As a simple example, the \nabla operator is linear, so if I encounter a call to
+after a call to
∇, then I can switch their order. How can I do such a thing in my AST?
- my variables
vfrom above have specific types. The fields of this type tell me, for example, what order derivatives I should use in my calculations. What is the most effective way of searching through the AST of my expression above looking for instances of my types?
Apologies if my question is a little too long winded, but I thought some context might help