Let $A$ be a $3×3$ matrix, a (possibly) complex function of $x, z \in \mathbb{R}$ representing the Order Parameter in the Ginzburg Landau equations, $$ A = \begin{pmatrix} A_{uu} & A_{uv} & A_{uw} \\ A_{vu} & A_{vv} & A_{vw} \\ A_{wu} & A_{wv} & A_{ww} \end{pmatrix}. $$ $K_i$,$\beta_i$, and $\alpha (T)$ are known parameters. I need to solve the following equation, (with the implied sum over $j$, for each component identified by $\mu$, $i$), $$ 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} $$

I have a code in C++ that implements the FDM with relaxation, but we have found that our mixed derivative approximations have large enough error that the code doesn't always converge to a solution. We are looking at the option of using others' FEM solvers (like FreeFEM or MOOSE), but I'm having a hard time getting the weak form of our set of equations. I am following the description here, and obtained $$ 0 = K_1 \oint_{\Gamma} {\psi \vec{\nabla}A_{\mu i} \cdot \hat{n}} - K_1 \int_{\Omega} {\vec{\nabla} \psi \cdot \vec{\nabla}A_{\mu i}} - \int_{\Omega} {\psi (\text{rhs})_{\mu i}} + K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})}, $$ but I don't know how to get the last term. I saw this post and thought that maybe I'd have to explicitly write out all 9 equations and use 9 test functions?

Edit: I played around with the equations a little more and figured I could write them generally as, $$ (\text{rhs})_{\mu i} = K_1 \vec{\nabla}^2 A_{\mu i} + K_{23} \partial_i \Big[ \vec{\nabla} \cdot \vec{A}_{\mu;r} \Big], $$ where $\vec{A}_{\mu;r}$ is the vector formed from the $\mu$-th row of $A$ (is there a better way to represent that?). Thus, $$ 0 = K_1 \oint_{\Gamma} {\psi \vec{\nabla}A_{\mu i} \cdot \hat{n}} - K_1 \int_{\Omega} {\vec{\nabla} \psi \cdot \vec{\nabla}A_{\mu i}} - \int_{\Omega} {\psi (\text{rhs})_{\mu i}} + K_{23}\int_{\Omega} {\psi \partial_i \Big[ \vec{\nabla} \cdot \vec{A}_{\mu;r} \Big]}, $$ which might be easier (or more obvious how) to integrate by parts?

So, would I still have to use 9 test functions? We expect the solution function for each element to be different, so I might still need to have all 9 separately...

Edit 2: Taking the last term, $K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})}$, explicitly writing out the sum over $j$, and using the product rule, we can say, $$ K_{23}\int_{\Omega} {\psi \partial_i (\partial_j A_{\mu j})} \to K_{23} \sum_{j=u,v,w}{ \Bigg[ \int_{\Omega}{\partial_i (\psi \partial_j A_{\mu j})} - \int_{\Omega}{(\partial_i \psi)(\partial_j A_{\mu j})} \Bigg] } $$

The first term here ($\int_{\Omega}{\partial_i (\psi \partial_j A_{\mu j})}$) looks like the divergence term that we previously rewrote as a boundary integral...except that there are no vectors. Is this now in the right form so that I can use it in a FEM solver?

  • $\begingroup$ You would be better off asking this on the scientific computing stackexchange. $\endgroup$ Jun 6, 2022 at 22:41
  • 1
    $\begingroup$ You need test functions also to be matrices, i.e. $\psi_{\mu i}$ and need summation by $i, \mu$ $\endgroup$
    – uranix
    Jun 7, 2022 at 12:29

1 Answer 1


To start with I strongly suggest not to mix vector and index notation. Let's start with indexed version of the integration by parts formula: $$ \int_{\Omega} (\partial_i u_{\alpha,\beta,\dots}) v_{\lambda,\mu,\dots} d\Omega = -\int_{\Omega} u_{\alpha,\beta,\dots} (\partial_i v_{\lambda,\mu,\dots}) d\Omega + \int_{\Gamma} n_i u_{\alpha,\beta,\dots} v_{\lambda,\mu,\dots} d\Gamma. \tag{*} $$ Indices $\alpha, \beta, \dots$ and $\lambda, \mu, \dots$ are arbitrary and may include $i$. Note that in the last term $\int_\Omega \partial_i $ became $\int_\Gamma n_i$.

The strong form of the equation is: $$ K_1 \partial_j^2 A_{\mu i} + K_{23} \partial_i (\partial_j A_{\mu j}) = ({\rm rhs})_{\mu i}. $$ Weak form is obtained by multiplying with $\psi_{\mu i}$, summation by $\mu, i$ and integration over $\Omega$: $$ K_1 \int_{\Omega} \psi_{\mu i} \partial_j^2 A_{\mu i} d\Omega + K_{23} \int_{\Omega} \psi_{\mu i} \partial_i (\partial_j A_{\mu j}) d\Omega = \int_{\Omega} \psi_{\mu i} ({\rm rhs})_{\mu i} d\Omega. $$ This form is not good enough since it contains second order derivatives. Let's eliminate them using (*) $$ -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. $$

Note that $\int_{\Omega}$ terms are symmetric by $\psi \leftrightarrow A$, which means that their discretization also would be a symmetric matrix.

You did not specify your boundary conditions, they are required to process further the $\int_\Gamma$ terms. In case of Dirichlet boundary conditions they simply vanish (due to $\psi\big|_\Gamma = 0$). In Neumann case they would be $$ K_1 n_j (\partial_j A_{\mu i}) + K_{23} n_i (\partial_j A_{\mu j}) = g_{\mu i} $$ and the integrals migrate to the right hand side as $\int_\Gamma \psi_{\mu i} g_{\mu i} d\Gamma$ term.

  • $\begingroup$ This is very helpful! Just to clarify a few things: the $K_1$ integral over $\Gamma$ should have $d\Gamma$? The $\int_{\Omega}$ terms aren't quite symmetric since one $A$ has an $i$ index while the other has $j$? For mixed boundary conditions, the integrals can be trivially broken up to apply different kinds on different domain boundaries? And the $n_i$ represents what, in this case with matrices? $\endgroup$
    – Izek H
    Jun 7, 2022 at 17:03
  • 1
    $\begingroup$ That was a copy-paste typo, coorected it, thank you. Each of the terms $\int_\Omega$ is symmetric with respect to interchanging of $A$ and $\psi$. It might not be so obvious for the second, but note that $i$ and $j$ are dumb indices. Mixed boundary conditions are similar to Neumman, you eliminate derivatives using boundary conditions. But in this case there remain terms like $\int_\Gamma \psi A d\Gamma$ in the left hand side. The vector $n_i$ is the unit outward normal to $\Omega$, it depends only on the domain and not on the trial or test function. $\endgroup$
    – uranix
    Jun 7, 2022 at 17:42
  • $\begingroup$ Uranix, this is somewhat unrelated, but I am also looking for a good FE solver that can solve these matrix-valued equations. Any suggestions? $\endgroup$
    – Izek H
    Jun 8, 2022 at 20:47
  • 1
    $\begingroup$ @IzekH Sorry, I've only used FEniCS. The documentation states that it has TensorFunctionSpace, so it might already support tensor-valued trial and test functions. But beware, its API changes often so many examples are outdated and do not work anymore. The documentation covers only some simple cases, for others one needs to dig into the autogenerated API docs and framework code itself. $\endgroup$
    – uranix
    Jun 8, 2022 at 20:58

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .