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.

  • $\begingroup$ See "The finite element method: Linear static and dynamic finite element analysis" by TJR Hughes. It has a very good explanation of finite elements for linear elasticity, including implementation. $\endgroup$
    – NNN
    Commented Dec 27, 2021 at 3:40

2 Answers 2


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.

  • $\begingroup$ Zoltàn, thanks a lot for your answer and the advice, it's really precious. Indeed, I just want a throw-away code in 2D, to be sure I really know what's going on. The approach I'd like to follow is the one where each basis function of the vector element has one non-zero component, like the step-8 you linked (I don't like the way engineers form matrices, to be honest). Do you have some reference also for that? (the books you mentioned follow the engineering approach) $\endgroup$
    – bob_bill
    Commented Dec 26, 2021 at 21:41
  • $\begingroup$ I don't remember any reference (book, lecture note etc.) for the vectorial basis function approach. All I saw considered scalar basis functions and vectorial degrees of freedom. However, building the element matrices is analogous to what I described in my answer. What changes is the approximation of the unknown (displacement field), but the rest (forming the derivatives of the unknown field, and then the terms of the inner product) is the same. $\endgroup$ Commented Dec 26, 2021 at 22:14
  • 1
    $\begingroup$ You can always check the source code of deal.ii to see what is going on. However, note that deal.ii is a general FE framework, so for your purpose (throwaway code, just for understanding) it will be too verbose. If your purpose is understanding, it might be better to form the element matrix for your specific problem as I described above, and implement it yourself. That's how you will learn the most. $\endgroup$ Commented Dec 26, 2021 at 22:20
  • $\begingroup$ Sure, deal.II is way too verbose for my simple task :-) My problem with the approach you suggest is that it is not "general", meaning that if I want to solve the mixed laplace equation, it doesn't work. That's why I want to stick with the other one $\endgroup$
    – bob_bill
    Commented Dec 26, 2021 at 22:28
  • $\begingroup$ Could you elaborate why taking scalar basis/shape functions and vectorial degrees of freedom do not generalize to mixed formulations? $\endgroup$ Commented Dec 26, 2021 at 22:44

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.

  • 1
    $\begingroup$ Thanks a lot for your answer and resources, they've been really useful ! As a beginner, I wanted to stick with the approach presented in deal.II, where they use vector basis functions that have only one non-zero component. $\endgroup$
    – bob_bill
    Commented Dec 29, 2021 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.