From this mathwork page "c for system", the usual second order PDE is written in tensor form:
$$ -\nabla\cdot(\mathbf{c} \otimes \nabla \mathbf{u})+\mathbf{a}\mathbf{u}=\mathbf{f} $$
and $\mathbf{c}$ is claimed to be a $N$-by-$N$-by-$2$-by-$2$ tensor. Each side of the equation should correspond to a $N$-by-$1$ vector in this case, and now I am confused... I want to express this equation in matrix form, since $\mathbf{c}$ is a 4th rank tensor while $\nabla{u}$ is 2nd order tensor($2$-by-$N$ matrix?), how to define the tensor operation $\otimes$ then? Seems wiki doesn't give the answer.
Cand anybody elaborate the rule here? (Also how to cast the gradient operator in this case?) Any physical insight is honestly appreciated!
Using abstract index notation:
There are two types of indices here: I will use $i,j,\ldots$ for indices relative to the domain (which appears to be $\mathbb{R}^2$) and $A,B,\ldots$ for indices relative to the co-domain (basically from the number of equations) which is $\mathbb{R}^N$. We can write $\mathbf{c}$ as the rank 4 tensor in index notation
$$ c^{A;ij}_{B} $$
(I included a semicolon just for clarity to separate the upper case and lower case indices.) and adopting the Einstein summation convention your equation is
$$ - \nabla_i ( c^{A;ij}_{B} \nabla_j u^B) + a^A_B u^B = f^A $$
Without using indices:
Your tensor $\mathbf{c}$ can be thought of a $(N\times 2)\times (N\times 2)$ matrix, that is a square matrix with $N\times 2$ rows and columns. You can flatten $\nabla \mathbf{u}$ to an $N\times 2$ column vector. Then $\mathbf{c}$ acts on $\nabla \mathbf{u}$ by matrix multiplication.
The key, of course, is to flatten things the same way from $\mathbf{c}$ and $\mathbf{u}$.