FEM for vector valued problems: reference request

Question

I'm a PhD working in a computational mechanics lab. I come from a Math department, and I have a good background for what concerns the basics of finite elements, like inf-sup conditions, DG, non-conforming FEs. I also did some basic codes, like Poisson on a square/L-shape/etc. (or other elliptic problems) using the classical building blocks of a code, like local to global mapping, quadrature formulas on the reference, etc.

However, there's something that I really don't know how to do, and this is how to solve vector valued problems. I know how to solve them in deal.II, for instance, and I watched the related video-lecture about it. But I do really need to re-invent the wheel (in 2D) this time. Let me try to explain: when we're using our software(s) to solve scalar problems, I really know what is going on inside the assembly routine and I know I would be able to replicate it (of course not in terms of efficiency), but that's not the case for the linear elasticity equation, for instance. I know that the test functions are vector valued, but I'm lacking the ability to put this into code. In particular, I'm interested in the approach where each vector basis function has only one non-zero component.

So I am looking for some reference (lecture notes, books, whatever) where I can find a simple and didactic explanation about the way the code has to be organised for this equation and the relative code. If there's something with MatLab that would be perfect, as I only need to understand the basic building blocks. However, I can read with no problems also Python and C++ codes.

Answer 1

Short answer: Just replicate the vector of interpolation functions into a block-diagonal matrix, as showed e.g. on page 5 in this lecture note.

Detailed answer: Mathematically oriented texts typically do not bother about detailing the implementation, that is why you haven't met it. On the other hand, engineering FEM do detail the implementation (but lack the theory). The good news is that understanding the implementation after the theory is significantly easier than the other way around.

Linear elasticity is a good example because (in the compressible case), it is typically implemented by Lagrange elements. You can find hundreds of lecture notes and engineering-oriented finite element books (e.g. K. J. Bathe: Finite Element Procedures, Zienkiewicz: The Finite Element Method: Its Basis and Fundamentals), which go into details about forming the element matrices, performing the assembly, doing the numerical integration, etc. To react to your specific questions:

1. What changes in vector-valued problems (e.g. elasticity equation) compared to scalar-valued problems (e.g. heat equation) is that you need to construct the approximation for each unknown function. If you haven't done yet, check the vector-valued problems in deal.ii for a thorough discussion. When you have independent fields to approximate (e.g. in mixed formulations), you often want to apply different interpolation for the different fields. On the other hand, when you have a single vector-valued field (e.g. the displacement vector), you most frequently use the same interpolation for all the components of the vector unknown. Hence, the latter problem is a special case of the former. What you do in both cases is writing the approximation for each field/component. In mathematical-oriented texts, this is done with indices, see e.g. the elasticity tutorial of deal.ii, while engineers form the matrices too (cf. second slide on page 5 of the lecture note referred at the beginning).
2. The other difference in vectorial problems is the inner product. For instance, in elasticity, you need to evaluate $\int_{\Omega}\sigma:\varepsilon \,\mathrm{d}\Omega$, where $\sigma$ and $\varepsilon$ are second order tensors. What you do in this case is called vectorization. For tensors, the term you should look for is the Voigt notation. This is done on the third slide on page 4 of the lecture note. Once you have rewritten all the fields in this matrix notation (in engineering texts, it is referred to as matrix formulation), you can factor out the element matrix from the inner product (first slide on page 8 of the lecture note). The assembly then goes on as for the scalar-valued case.

Footnote: These were just the core ideas, which help you understand the method. However, I must tell you that the devil is in the details. There are so many other things to take care of when implementing the FEM, that your implementation will almost surely be suboptimal in terms of maintainability, performance and scaling. Unless you want a throw-away code for a quick demonstration purpose only, I recommend you not to implement it by yourself (or at least not in 3D). I made this mistake in my PhD, wasting precious time. When I wrote my own FE code in MATLAB, I was partly inspired by mFEM. However, I suggest that you rely on robust, high-quality and actively-maintained libraries and frameworks, such as deal.ii or FEniCS (there are dozens of others, the choice should depend on your use case). Nevertheless, deriving the discrete formulation on paper is a good idea to completely understand your problem at hand.

Answer 2

We can show how it works on the example of linear elasticity. In classical finite elements formulations, on every node, we will have a scalar shape (base) function to which we have associated number coefficients equal to the dimension of the problem. Coefficients at nodes are interpreted as physical nodal displacements.

Let displacements are are approximated by scalar base functions, \begin{equation} u^h_\alpha = \sum_{\alpha=0}^N \phi^\alpha c^\alpha_i, \end{equation} where $\phi^j$ conforming base function on node $\alpha=0..N$, where $N$ is number of nodes and $c^\alpha_i$ are cofficients (i.e. nodal displacements) and $i=0..dim$ are directions in the Cartesian coordinate system, wher $dim$ is dimension of problem. Then gradient of displacements is \begin{equation} u^h_{i,j} = \sum_{\alpha=0}^N \phi^\alpha_{,j} c^\alpha_i, \end{equation} where $()_{,j}$ is derivative, that is \begin{equation} \phi^\alpha_{,j}=\frac{\partial \phi^\alpha}{\partial x_j}. \end{equation}

With above at hand, we can calculate the stiffness matrix from element virtual work, as follows \begin{equation} \begin{split} \delta W^e &= \delta c^\alpha_i F^\alpha_i \\&= \delta c^\alpha_i \left( \int_{V^e} \phi^\alpha_{,j} \sigma^h_{ij} \textrm{d}V \right) \\&= \delta c^\alpha_i \left( \int_{V^e} \phi^\alpha_{,j} D_{ijkl} u^h_{k,l} \textrm{d}V \right) \\&= \delta c^\alpha_i \left( \int_{V^e} \phi^\alpha_{,j} D_{ijkl} \phi^\beta_{,l} \textrm{d}V \right) c^\beta_k \\&= \delta c^\alpha_i \overline{M}^{\alpha\beta}_{ik} c^\beta_k, \end{split} \end{equation} where we exploit minor symmetry of material stiffens tensor $D_{ijkl}$, $V^e$ is volume of element and $\delta c^\alpha_i$ virtual nodal displacements at node $\alpha$ and direction $i$. $F^\alpha_i$ is internal nodal force vector. Einstein summation convention applies for repeating indices. For problem in $dim$ dimensions, local element matrix $M_{IJ}$, will have size $I,J=0..dim*N$. Thus we understand $\overline{M}^{\alpha\beta}_{ij}$ as subblock $dim$ by $dim$ of element matrix, as follows \begin{equation} M_{IJ} = \overline{M}^{\alpha\beta}_{ij}, \quad\textrm{where}\; I=dim*\alpha+i,\,J=dim*\beta+j. \end{equation}

Note: Indexes start from "$0$".

Note: We can define material stiffens tensor $\mathbf{D}$ for the Hooke model as follows \begin{equation} D_{ijkl} = G \left[ \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} \right] + A (K - \frac{2}{3} G) \left[ \delta_{ij} \delta_{kl} \right], \end{equation} where $K$ and $G$ are bulk and shear modulus, respectively. For 3D problem or planse strain porblem $A=1$, and for plane strain 2D problem whereas for plane stress it takes the following form: \begin{equation} A=\frac{2 G}{K+\frac{4}{3} G}. \end{equation}

Note: Look at Linear Elasticity Video link where my colleague explains how to implement general dimension implemented case for linear/nonlinear elasticity using forms integrator implemented as above. Disclaimer: I am the developer of MoFEM.

Note: Example of form integrator implemented as above is under following link. Implementation using FTensor library, and exploiting properties of minor symmetries $\mathbf{D}$ tensor.