# What FEM solver should be used for matrix-valued FE spaces?

I am pretty new to using FE solvers. I am trying to solve a system of (up to) 9 complex equations. We write these as a matrix equation (here), (with the implied sum over $$j$$, for each component identified by $$\mu$$, $$i$$), and known parameters (or at least constant wrt position): $$K_i$$, $$\beta_i$$, and $$\alpha (T)$$; $$K_1 \partial^2_j A_{\mu i} + K_{23} \partial_i (\partial_j A_{\mu j}) = 2 \beta_1 Tr(AA^T) A_{\mu i} + 2 \beta_2 Tr(AA^\dagger)A_{\mu i} + 2 \beta_3[AA^TA]_{\mu i} + 2 \beta_4[AA^\dagger A]_{\mu i} + 2 \beta_5[AA^TA]_{\mu i} + \alpha(T) Tr(AA^\dagger) = (\text{rhs})_{\mu i},$$ assuming, $$A = \begin{pmatrix} A_{uu} & A_{uv} & A_{uw} \\ A_{vu} & A_{vv} & A_{vw} \\ A_{wu} & A_{wv} & A_{ww} \end{pmatrix}.$$

Thanks to uranix, I have the weak form of the equations, $$-K_1 \int_{\Omega} (\partial_j \psi_{\mu i}) (\partial_j A_{\mu i}) d\Omega -K_{23} \int_{\Omega} (\partial_i \psi_{\mu i}) (\partial_j A_{\mu j}) d\Omega\\ +K_1 \int_{\Gamma} \psi_{\mu i} n_j (\partial_j A_{\mu i}) d\Gamma +K_{23} \int_{\Gamma} \psi_{\mu i} n_i (\partial_j A_{\mu j}) d\Gamma = \int_{\Omega} \psi_{\mu i} ({\rm rhs})_{\mu i} d\Omega.$$

Question: What FE solver would be the best to solve this set of equations?

I have been learning FreeFEM, and it seems like it can do this by using "vectorial FE spaces" (like they briefly mention here), and define functions to convert between arrays and matrices. I have been unsuccessful in declaring the trial and test functions as members of such a "vectorial FE space." Does FreeFEM have anything else like it that I could use? I have also done a little with MOOSE, but I'm not sure how I could make it work. I looked into deal.II, and it looks like I could totally make it work, even if I have to build my own class in C++, but it doesn't build right with CMake (complains that it can't find some serialization library in Boost). The last one I've looked at is NGSolve, and one person talked about a system of 5 equations, and here the documentation mention a "matrix-valued function." I just don't want to download and build yet another FEM solver if it might not work.

• For what concerns your problem with CMake, you may ask for help in the deal.II mailing list (groups.google.com/g/dealii)
– VoB
Jun 12 at 18:13

This is a fairly standard elliptic problem (the operator is of the form $$-\Delta - \text{grad}\,\text{div}$$). It is vector-valued (in fact matrix-valued) in the same way as the elasticity equations are vector-valued.
• You state that $$A$$ is a matrix quantity, but is there any requirement that it be a symmetric matrix? This substantially complicates the choice of finite element spaces.
• deal.II did look very promising! I'll have to see if I can find help with the CMake problems. And unfortunately, $A$ is not always symmetric. The systems that we are solving have some kinds of symmetry (like just in the form/structure of $A$), but $A_{uw}$ and $A_{wu}$ sometimes will be quite different. Jun 13 at 19:32
• The fact that you are looking for non-symmetric solutions actually makes it easier because enforcing symmetries is hard. As for coordinate systems: Yes, that only leads to a slightly different weak formulation (every integral will be over $r,\phi$ instead of $x,y$ and will have an extra $2\pi r$ as a weight in the integral. Jun 15 at 11:44