Usually I write the weak form by hand for my FEM code, but it's a little annoying and mechanic sometimes.
So, I wonder, is there any way to generate the symbolic weak form in Mathematica? For instance, if I have this strong form:
$$\nabla\cdot(\nabla u) = f$$
I'd like to have an automation of multiplying by a shape function $\phi$, and integrating it obtaining
$$\nabla \phi \ \nabla u = \phi f$$
So the procedure is as follows.
- take the left hand side and left-multiply with $\phi$: $\phi \; (\nabla\cdot(\nabla u))$
- integrate by parts (usually): $\int \phi \; (\nabla\cdot(\nabla u)) = \int \nabla\phi \cdot \nabla u - \int_\partial \phi \cdot \nabla u$
- by definition $\phi$ is zero on the boundary $\int_\partial \phi \cdot \nabla u = 0$
- hence only the left hand side (integrated) remains $\int \nabla\phi \cdot \nabla u$
The very same procedure could be applied to the right hand side.
This would be great to achieve symbolically, as I am planning to implement a FEM software, not using
NDSolve or other functions.
The top would be even generating the symbolic sum with Jacobians and Gauss integration points, but this is maybe asking too much! In a very classical FEM formulation, as for instance here in deal.II for a simple Laplace problem, $u$ and $\phi$ are just interpolated with a shape function, for instance $\phi$ itself, as $u(x) = \sum_j U_j \phi_j(x)$. The integral is then reduced to a sum of integrals over a triangulated domain, but as I said, I'd be glad to just extract the weak form in general (not just for linear problems, but also for instance for $\nabla \cdot ( A(x) \nabla u) = -f$ for a general $A(x)$.