# Practical way to build element stiffness matrix for 3D FEM simulations

The element stiffness matrix in 3D FEM problems is build as follows:

$$K = \int\limits_{[-1,1]^3} B^T C B\, |J| \mathrm{d}r\, \mathrm{d}s\, \mathrm{d}t$$

The integral can be solved using e.g. Gauss quadrature method.

But for me it seems to computational expensive if you perform the above equation for, e.g., 500.000 hexahedral elements. Thus I was wondering if there is an additional approximation or something else done to reduce the computational effort.

Question: Is there some kind of approximation method used to solve the two matrix-matrix products as well as the determinant or are they computed explicitly for each element?

• I don't understand you question. Have you implemented stiffness matrices before? For example, in 1D. – nicoguaro Apr 18 at 23:08
• Yes, I have implemented stiffness matrix befor as described in my post above. But for me it seems to computational expensive if you perform the above equation for e.g. 500.000 hexahedral elements. Thus I was wondering if there is an additional approximation or something else done to reduce the computational effort. – vydesaster Apr 19 at 7:19
• @vydesaster -- take a look at the tutorials of any of the large finite element libraries. They all do essentially this. It's a small fraction of what it costs to solve a linear system with the large matrix $K$. – Wolfgang Bangerth Apr 20 at 20:45
• Thanks for the information :) – vydesaster Apr 21 at 7:38

## 2 Answers

No that I'm aware of. In general the Jacobian is a function and it takes a different value for each Gauss point.

On the other hand, if your elements have the same shape, you can compute the integral once and then just assemble. This implies that you have a regular mesh and limits the geometries that you can model. One can, naively, think that for this case you can compute the integrals analytically but I don't think that it gives any advantage, as discussed in a previous answer.

Since you mention hexahedral cells in particular, there are optimizations you can do to reduce the computational cost by exploiting the tensor product structure of the element. See for example Antolin et al., "Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization" and references therein. This isn't so much an approximation as much as a clever way to group floating-point operations to reduce overall cost.