Many initial boundary value problems require solving for unknown vector fields, such as solving for displacements in a mechanical problem. Discretization of weak forms of such problems leads to higher-order linear systems which need to be reshaped to be solved by regular linear solvers. There are also more indices involved than a scalar problem, which can be confusing. In this post, I’ll try to elucidate the procedure by deriving for a basic higher-order system and giving an example.

The weak formulation of a linear vectorial problem reads

Find $\Bu\in V$ such that

for all $\Bv\in V$.

Discretizations of vectorial problems requires the expansion of vectorial quantities as linear combinations of the basis vectors $\Be_i$:

where $\cbr{u_i}_{i=1}^{\ndim}$ are the components corresponding to the basis vectors and $\ndim=\dim V$. Here, we chose Cartesian basis vectors for simplicity.

We can therefore express its discretization as

Substituting discretized functions in the weak formulation, we obtain

and

We define the following arrays

Hence we can express the linear system

as

For arbitrary $\Bv_h$, this yields the following system of equations

for $i=1,\dots,\ndim$ and $I=1,\dots,\nnode$.

We reshape this higher-order system as shown in the previous post Reshaping Higher Order Linear Systems:

by defining a map $i_d$ that maps original indices to the reshaped indices

where we used 1-based indexing of the arrays. We set

and write

The inverse index mapping can be obtained as shown in the previous post.

## Example: Linear Elasticity

Our initial-boundary value problem is

Find $\Bu\in V$ such that

for all $\Bv\in V$ where $V=H^1(\Omega)$.

We apply integration by parts on the left-hand side

and apply the divergence theorem to the first resulting term:

Substituting the linear stress $\Bsigma=\IC:\Bvareps=\IC:\nabla\Bu$, we obtain the following variational forms:

We have the following discretizations of the unknown function and test function

With the given discretizations, the matrix corresponding to \eqref{eq:linelastdiscretebilinear} can be calculated as

and finally obtain

The vector corresponding to \eqref{eq:linelastdiscretelinear} is calculated as

which yields