Building blocks for solving a vector valued problem

Question

This question is a follow-up of this previous one. I decided to solve the linear elasticity \begin{cases}- \nabla \sigma(u)=f \\ u=0 \text{ on } \partial \Omega\end{cases}

with P1 Lagrangian finite elements on the square $[-1,1]^2$. Here I want to highlight and ask for a confirmation the crucial steps, that were totally new to me. The convergence plot seems to indicate order $2$ in $L^2$ and $1$ in $H^1$, w.r.t the number of DoFs. I'll attach as less code as possible, but what I want to be sure is if my design is correct.

Discaimer: I know it's not efficient, etc. The purpose of this is to understand how a vector valued problem can be solved.

• I decided to use vector basis functions where each vector function has only a non-zero component. Therefore, in the reference triangle I have, for each node, 2 basis functions. Take as example the first node, indexed by $1$ (I'm using MatLab). There, I have $\Phi_1=[\phi_1(\boldsymbol{x}),0]^T$ and $\Phi_2=[0,\phi_1(\boldsymbol{x})]^T$. Since I have three nodes for $P1$, I have 6 basis functions in the reference.

• I know I have to use a linear index in order to enumerate all these functions. In particular, now I have twice the DoFs. In C++ (or Python), the choice I am aware is that the new indices are given by $$n=2N +d$$ where $d \in\{0,1\}$, $N$ is the global node number.

In MatLab, I changed this with $$n=2(N-1) + d \qquad (\star)$$ with $d \in \{1,2\}$

• The contributions coming from nodes $(i,j)$ are the classical ones: $$\int_\Omega \bar{\varepsilon}(\Phi_i) C \bar{\varepsilon}(\Phi_j)$$ where $C$ is the matrix $$C = \frac{E}{ 1 - \nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0& 0& \frac{1 - \nu}{2}\end{bmatrix}$$

and $ \bar{\varepsilon}=[ \varepsilon_{xx},\varepsilon_{yy},\varepsilon_{xy}]$ where we have the usual relations between strain and displacement $u=[u_1,u_2]$, namely: $$\varepsilon_{xx}(u) = \partial_x u_1 \\ \varepsilon_{yy}(u)= \partial_y u_2 \\ \varepsilon_{xy}(u)= \partial_{x} u_2 + \partial_{y} u_1$$

Assembly routine

• I loop over all triangles, and for each triangle $K$ I get the global indices for the three DoFs. Then, I convert those indices using $(\star)$ and hence my local_to_global has 6 components that have to be distributed.

• Then, inside the quadrature loop, I collect the values and the gradients of each vector basis function $\Phi_i$. To achieve what I wrote at the beginning, I basically have $$\Phi_i = \phi_{i//2}(\boldsymbol{x}) \boldsymbol{e}_{i\%2}$$ In code, this is (fix does the integer division)

    %S are the values of shape functions at quadrature points
    %dSdx,dSdy gradients of shape functions compute via isoperimetric map
    for comp=1:6
        if(mod(comp,2)==0)
            value = [0;S(fix(comp/2))]; %test function is [0,v_comp]
            gradient = [0,0;dSdx(fix(comp/2)),dSdy(fix(comp/2))];
        elseif(mod(comp,2)==1)
            value = [S(fix(comp/2)+1);0]; %test function is [v_comp,0]
            gradient = [dSdx(fix(comp/2)+1),dSdy(fix(comp/2)+1);0,0];
        end
        Phis = [Phis,value];
        Gradients = [Gradients,gradient];
    end
  

So, $\Phi_3=[\phi_2,0]$ and $\Phi_4=[0,\phi_2]$ are the two shape functions associated to the 2nd node on the reference triangle.

• Then, I looped over $i,j$ running both from 1 to 6 and distributed the contributions. In MatLab code, it looks like this:

        for i=1:6
        phi_i = Phis(:,i);
        grad_phi_i = Gradients(:,[2*(i-1)+1,2*(i-1)+2]);
        grad_v_1 = grad_phi_i(1,:);
        grad_v_2 = grad_phi_i(2,:);
        for j=1:6
            grad_phi_j = Gradients(:,[2*(j-1)+1,2*(j-1)+2]);
            grad_u_1 = grad_phi_j(1,:);
            grad_u_2 = grad_phi_j(2,:);
  
            %here sigma, sigma_2 are E/(1-nu^2) and (1-nu)/2 coefficients
            val = sigma * (...
                grad_v_1(1) * (grad_u_1(1) + nu * grad_u_2(2))...
                + ...
                grad_v_2(2) * (nu * grad_u_1(1) + grad_u_2(2)) ...
                + ...
                (grad_v_1(2) + grad_v_2(1)) * ...
                (sigma_2 * (grad_u_1(2) + grad_u_2(1)))...
                ) * JxW;
  
            A(l2g(i),l2g(j)) = A(l2g(i),l2g(j)) + val;
        end
        val_rhs = dot(force(xp,yp),phi_i)*JxW;
        F(l2g(i)) = F(l2g(i)) + val_rhs;
    end
  

After the global matrix has been built, I just zero-out rows and columns corresponding to boundary DoFs and place a $1$ in the diagonal and a $0$ in the forcing term as I have homogeneous Dirichlet.

Test

I chose the well known test-case for which the solution is $u=[(x^2-1)(y^2-1), (x^2-1)(y^2 - 1)]$, i.e. $$\boldsymbol{f}=\frac{E}{1 - \nu^2}[-2y^2-x^2+ \nu x^2 - 2 \nu xy - 2 xy + 3 - \nu,-2x^2 - y^2+ \nu y^2 - 2 \nu xy - 2 xy + 3 - \nu]$$ and here's the EOC in $L^2, H^1$ for the first component.

Answer 1

Your assembly procedure seems correct to me. I made some edits in your post to correct some possible typos. Feel free to raise objections. Looking at the graphs and their legends, I see the rate of convergence as $\mathcal{O}(h)$ and $\mathcal{O}(\sqrt{h})$ for the $L^2$ and $H^1$ norms.

I propose the following extensions to test your understanding:

1. Run your code with different values of $\nu$. What do you experience when $\nu \to 0.5$ ?
2. Inhomogeneous Dirichlet BC That will require a bit more modification in the assembly. Here, you can experiment with various assembly techniques:
   1. Modifying the assembled matrix This is what you have done so far for the homogeneous case.
   2. Enforce the boundary condition with Lagrange multipliers.
   3. Weakly enforce the BC using the penalty or the Nitche method and investigate how well the BC is satisfied, and what is the effect of the penalty/stabilization parameter on the conditioning of your linear system.
   4. Modify the element matrices to incorporate the BC. Note that in this case, you do not need to modify the global matrix.
3. The convergence plots are promising when verifying the code, but you could also create another manufactured solution in order to test if the constant stress state can be represented. This is a good indicator to see if the partition of unity property holds.
4. The convergence plots are normally shown for the vector variable, not for the individual components. You could try computing other error indicators, e.g. the energy error, which does not need any post-processing (you already have the assembled matrix)!
5. Change the local dof numbering within the elements. Also try a numbering scheme in which the dofs are not grouped by nodes but by components (i.e. $[u_{1x}, u_{2x}, u_{3x}, u_{1y}, u_{2y}, u_{3y}]$ for the $P_1$ Lagrange element). How does it affect the sparsity pattern of your global matrix?
6. Change the global node numbering. How does it affect the sparsity pattern of your global matrix? Note that when you use a black-box direct solver (e.g. the backslash in MATLAB), it will automatically rearrange your matrix for a favourable sparsity structure, so you do not need to bother much about sparsity in your basics implementation.

If you want to practice assembly more, or you would like to be informed more about implementation issues, you might find the following tasks interesting.

1. Consider a non-overlapping domain decomposition (construct a very simple setup with only a few elements so that you can do it on paper) of your original domain. What changes here is that in addition to the two mappings that existed in the standard FEM, you will now have an intermediate mapping. So it will look similar to this: element-wise node-dof mapping --> map element dof to subdomain --> map subdomain to the whole domain In actual parallel implementations, there are so-called ghost nodes, for communication among the subdomains.
2. Try the assembly for a DG method.
3. Implement a node-based assembly. This non-conventional assembly strategy has the advantage that there is no race condition during the parallel assembly of the global matrix.
4. You mentioned in your previous post that you are not much interested in the matrix notation used by engineers. Nevertheless, you may profit from it in the future if you understand both the index notation (assembly with nested for loops, as done currently) and the matrix formalism. Dealing with block matrices is definitely useful in many situations. Depending on the problem at hand, sometimes the index notation and sometimes the matrix notation is more convenient (readable) or efficient in terms of memory and CPU. For array-languages, such as MATLAB, you have to sacrifice readability to gain speed. This article might interest you.
5. Not strictly related, but the assembly of a global matrix can be skipped if you compute the contributions on-the-fly. This saves a lot of memory (no need to store a global matrix), but requires the recomputation of certain integrals. This strategy is used for low-order elements, where recomputing terms is cheap so it is not the bottleneck.

Hopefully, you can now have a glimpse of the intricacies of the robust implementation of the FEM. And this is just the surface...