System of equations for vector valued functions problems I am having trouble about determining the system to solve when using finite element method for vector fields.
Suppose we have the following problem:
$$-\Delta\mathbf{u}=\mathbf{f}\;\text{in}\;\Omega,$$ 
with some Dirichlet boundary conditions. The weak form would be:
$$\int_{\Omega}\nabla\mathbf{u}:\nabla\mathbf{v}=\int_{\Omega}\mathbf{fv},$$
where $\nabla\mathbf{u}$ and $\nabla\mathbf{v}$ are matrices.
If we discretize the functions like this:
$$\mathbf{u}^h=\sum_{j=1}^{N}\mathbf{u}_j\mathbf{\phi}_j,$$ where $N$ are the nodes of the discretization. Here is my first doubt, should be $\mathbf{\phi}_j$ scalar or vector valued function? In the first case I arrive to the following system of equation:
$$\begin{pmatrix}A & 0\\
0 & A\end{pmatrix} \begin{pmatrix}U^1\\U^2
\end{pmatrix}= \begin{pmatrix}f_1\phi_i\\
f_2\phi_i\end{pmatrix},$$where $A = \int\nabla\phi_i\nabla\phi_j$ and $f_1\phi_i,\;i,j=1...N$, is it correct?.
However, if $\mathbf{\phi}_j$ is a vector valued function then the discretization $u^h=\sum_{j=1}^{N}\mathbf{u}_j\cdot\mathbf{\phi}_j$ is not a vector but a scalar so I do not know what to do. Moreover, I do not know if it makes sense that $\mathbf{\phi}_j$ be a vector valued function. 
Thanks!
 A: Recall that the $\phi_i$ come from a specific space of functions which is related to the dimensionality of $\Omega$ and other things. You should be careful to spell out what space these functions come from as it is a necessary part of the finite element formulation.
If your problem is one dimensional then $\phi$ should indeed be scalar valued. If you are interested in solving the problem in higher dimensions the way one normally proceeeds is by multiplying by a seperate test function for each dimension. One can then go through the details and recombine the multiple equations into one later, but you will still have multiple test functions cooresponding to each trial space and each of these new higher dimensional $'\phi_i'$ will be in a product space of the lower dimensional spaces.
Also note that your linear system is missing the variables that you intend to solve for! It is usually notated as $Au=F$, where $u$ is your vector of variables.
A: Finally, I found how to discretize the problem for vector valued functions. The following annotations were inspired by the great book Finite Element Methods for Flow Problems writen by Jean Donea and Antonio Huerta.
Let's define the canonical basis:$$\mathbf{e_1}=(1,0)^T,\;\mathbf{e_2}=(0,1)^T.$$
Then, the solution can be discretized as:$$\mathbf{u^h}=\sum_{i=1}^{2}u_i^h(\mathbf{x})\mathbf{e_i}=\sum_{i=1}^{2}\sum_{j=1}^{N}u_{ij}\phi_j(\mathbf{x})\mathbf{e_i},$$
where $N$ is the number of nodes in the discretization and $u_{ij}$ is the approximation of the solution for the node $j$ in dimension $i$ (not spatial dimension but the dimension of the vector valued function output).
So, the system to solve for the problem $$\int_{\Omega}\nabla\mathbf{u}:\nabla\mathbf{v}=\int_{\Omega}\mathbf{fv}\Rightarrow a(\mathbf{u},\mathbf{v})=l(\mathbf{v})$$is:
$$\sum_{i=1}^{2}\sum_{j=1}^{N}u_{ij}\;a(\phi_j\mathbf{e_i},\phi_s\mathbf{e_c}) = l(\phi_s\mathbf{e_c}),$$
and there are $N\times 2$ equations, that is, and equation for each $e$ and $s$ where $s=1,...,N$ and $c=1,2$.
I would like to say that @MSEoris was right about saying that the problem itself has no physical sense however I asked the question having in mind a Stokes problem and in order to summarize it I only wrote the laplacian part.
