Interactive Lecture
Element Types, Shape Functions, and Isoparametric Formulation
by Dr Arnaldo Delli Carri
The power of FEM lies in its ability to discretise complex geometries using different element types. Once you see how shape functions extend from 1D to 2D and 3D, the whole framework clicks into place — the element type changes, but the process does not.
The pattern that holds throughout: All FEM formulations approximate the solution using shape functions defined over elements, derive element equations from the weak form or variational principle, and assemble into a global system. The complexity increases with dimensionality but the process remains the same.
The simplest finite elements are one-dimensional:
Linear 2-Node Bar Element:
\[u(\xi) = N_1(\xi)u_1 + N_2(\xi)u_2\]
Shape functions in natural coordinates \(\xi \in [-1, 1]\):
\[N_1(\xi) = \frac{1 - \xi}{2}, \quad N_2(\xi) = \frac{1 + \xi}{2}\]
Quadratic 3-Node Element:
\[N_1 = \frac{\xi(\xi-1)}{2}, \quad N_2 = 1 - \xi^2, \quad N_3 = \frac{\xi(\xi+1)}{2}\]
The middle node provides curvature, improving accuracy.
Two-dimensional elements can be triangular or quadrilateral:
3-Node Triangle (Linear):
Natural coordinates: \(L_1, L_2, L_3\) (area coordinates)
\[N_i = L_i, \quad L_1 + L_2 + L_3 = 1\]
Displacement field: \(u = N_1 u_1 + N_2 u_2 + N_3 u_3\)
4-Node Quadrilateral (Bilinear):
Natural coordinates: \(\xi, \eta \in [-1, 1]\)
\[N_1 = \frac{(1-\xi)(1-\eta)}{4}, \quad N_2 = \frac{(1+\xi)(1-\eta)}{4}\]
\[N_3 = \frac{(1+\xi)(1+\eta)}{4}, \quad N_4 = \frac{(1-\xi)(1+\eta)}{4}\]
The isoparametric formulation uses the same shape functions for both the element geometry and the field variables.
Coordinate Mapping:
\[x(\xi, \eta) = \sum_{i=1}^n N_i(\xi, \eta) x_i\]
\[y(\xi, \eta) = \sum_{i=1}^n N_i(\xi, \eta) y_i\]
Field Variable Approximation:
\[u(\xi, \eta) = \sum_{i=1}^n N_i(\xi, \eta) u_i\]
This allows curved elements to fit complex geometries accurately.
Derivatives in physical space require the Jacobian matrix:
\[[J] = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{bmatrix}\]
Derivatives transform as:
\[\begin{bmatrix} \frac{\partial N}{\partial x} \\ \frac{\partial N}{\partial y} \end{bmatrix} = [J]^{-1} \begin{bmatrix} \frac{\partial N}{\partial \xi} \\ \frac{\partial N}{\partial \eta} \end{bmatrix}\]
Integration requires: \(dx\,dy = \det(J)\,d\xi\,d\eta\)
FEM typically uses the weak (or variational) form of PDEs:
Strong Form (Heat equation):
\[-\kappa\nabla^2 T = Q \quad \text{in } \Omega\]
where \(\kappa\) is the thermal conductivity, \(T\) the temperature, and \(Q\) the heat source.
Weak Form (after integration by parts):
\[\int_\Omega \kappa \nabla w \cdot \nabla T\,d\Omega = \int_\Omega wQ\,d\Omega + \int_{\partial\Omega} wq\,d\Gamma\]
where \(w\) is the weighting function (test function).
Galerkin Method: Choose \(w = N_i\) (shape functions)
This leads to the element stiffness matrix:
\[[K^e]_{ij} = \int_{\Omega^e} \kappa \nabla N_i \cdot \nabla N_j\,d\Omega\]
Gauss Quadrature: Integrals are evaluated numerically using Gauss points.
For 1D: \(\int_{-1}^{1} f(\xi)\,d\xi \approx \sum_{i=1}^{n_g} w_i f(\xi_i)\)
For 2D: \(\int_{-1}^{1}\int_{-1}^{1} f(\xi,\eta)\,d\xi\,d\eta \approx \sum_{i=1}^{n_g}\sum_{j=1}^{n_g} w_i w_j f(\xi_i, \eta_j)\)
Common schemes: 2×2 for bilinear quads, 3×3 for biquadratic quads.
Element quality affects accuracy and convergence:
Requirements for Convergence:
Convergence Rate: For polynomial order p, error decreases as \(O(h^{p+1})\) where h is element size.
Explore different element types, mesh refinement, and shape functions. Click on elements to see their shape functions evaluated.
Number of Elements: 4
Number of Nodes: 9
Degrees of Freedom: 18
Element Order: Linear
The Jacobian determinant's sign depends on the ordering of element nodes. For a valid element:
There are two standard approaches to resolve negative Jacobian determinants:
1. Correct the Node Ordering: Reorder element nodes to restore counter-clockwise orientation in physical space. For the cylinder, this would mean reversing certain node pairs.
2. Use Absolute Value (with caution): Taking |det(J)| allows the calculation to proceed, but this:
Important: The magnitude of det(J) does NOT directly indicate element quality!
Common misconceptions about the Jacobian determinant:
What det(J) represents:
Element quality is determined by:
Problems caused by negative or zero Jacobian determinants:
Non-uniform Jacobian determinants (even if all positive) indicate:
Professional FEM software will detect and reject meshes with negative Jacobian elements, requiring the user to fix the mesh before analysis can proceed.
Advanced FEM formulations enable analysis of complex engineering systems:
| Concept | Equation | Description |
|---|---|---|
| Linear 1D Shape Function | \(N_1 = \frac{1-\xi}{2}, \quad N_2 = \frac{1+\xi}{2}\) | 2-node bar element |
| Quadratic 1D Shape Function | \(N_1 = \frac{\xi(\xi-1)}{2}, \quad N_2 = 1-\xi^2\) | 3-node element with mid-node |
| Q4 Shape Function | \(N_i = \frac{(1+\xi_i\xi)(1+\eta_i\eta)}{4}\) | Bilinear quadrilateral |
| Isoparametric Mapping | \(x = \sum N_i x_i, \quad y = \sum N_i y_i\) | Geometry from same shape functions |
| Jacobian Matrix | \([J] = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{bmatrix}\) | Coordinate transformation |
| Stiffness Matrix | \([K^e]_{ij} = \int_{\Omega^e} [B_i]^T[D][B_j]\det(J)\,d\xi\,d\eta\) | Element stiffness (general): \([B_i]\) = strain–displacement matrix, \([D]\) = material constitutive matrix |
| Gauss Quadrature | \(\int_{-1}^{1} f(\xi)d\xi \approx \sum w_i f(\xi_i)\) | Numerical integration |
For a 4-node quadrilateral element with nodes at (0,0), (2,0), (2,1), (0,1), evaluate: (a) all four shape functions at the element centre (ξ=0, η=0), (b) the physical coordinates at this point, and (c) verify partition of unity.
Show Solution
(a) Shape functions at centre (ξ=0, η=0):
\[N_1 = \frac{(1-0)(1-0)}{4} = \frac{1}{4}\]
\[N_2 = \frac{(1+0)(1-0)}{4} = \frac{1}{4}\]
\[N_3 = \frac{(1+0)(1+0)}{4} = \frac{1}{4}\]
\[N_4 = \frac{(1-0)(1+0)}{4} = \frac{1}{4}\]
(b) Physical coordinates:
\[x = N_1(0) + N_2(2) + N_3(2) + N_4(0) = \frac{1}{4}(0 + 2 + 2 + 0) = 1\]
\[y = N_1(0) + N_2(0) + N_3(1) + N_4(1) = \frac{1}{4}(0 + 0 + 1 + 1) = 0.5\]
Result: (1, 0.5) is indeed the centre of the rectangle.
(c) Partition of unity:
\[\sum_{i=1}^4 N_i = \frac{1}{4} + \frac{1}{4} + \frac{1}{4} + \frac{1}{4} = 1 \quad ✓\]
Answer: (a) \(N_1 = N_2 = N_3 = N_4 = \tfrac{1}{4}\); (b) physical point \((x, y) = (1,\ 0.5)\); (c) \(\sum N_i = 1\), so partition of unity holds.
For the element in Problem 1, calculate the Jacobian matrix at (ξ=0, η=0) and find its determinant. What does this value represent physically?
Show Solution
Shape function derivatives in natural coordinates:
\[\frac{\partial N_1}{\partial \xi} = -\frac{(1-\eta)}{4}, \quad \frac{\partial N_1}{\partial \eta} = -\frac{(1-\xi)}{4}\]
At (ξ=0, η=0):
\[\frac{\partial N_1}{\partial \xi} = -\frac{1}{4}, \quad \frac{\partial N_1}{\partial \eta} = -\frac{1}{4}\]
Jacobian components:
\[\frac{\partial x}{\partial \xi} = \sum \frac{\partial N_i}{\partial \xi}x_i = -\frac{1}{4}(0) + \frac{1}{4}(2) + \frac{1}{4}(2) + (-\frac{1}{4})(0) = 1\]
\[\frac{\partial x}{\partial \eta} = \sum \frac{\partial N_i}{\partial \eta}x_i = -\frac{1}{4}(0) + (-\frac{1}{4})(2) + \frac{1}{4}(2) + \frac{1}{4}(0) = 0\]
\[\frac{\partial y}{\partial \xi} = 0, \quad \frac{\partial y}{\partial \eta} = 0.5\]
Jacobian matrix:
\[[J] = \begin{bmatrix} 1 & 0 \\ 0 & 0.5 \end{bmatrix}\]
\[\det(J) = 1 \times 0.5 - 0 \times 0 = 0.5\]
Physical meaning: det(J) = 0.5 is the local scaling factor between natural and physical coordinates — the ratio of physical area to natural area, not a fraction of the physical area itself. Here physical area = 2×1 = 2 and natural area = 2×2 = 4, so det(J) = 2/4 = 0.5. For any quadrilateral the natural area is always 4, so det(J) rescales that fixed reference area onto the physical element.
Answer: \([J] = \begin{bmatrix} 1 & 0 \\ 0 & 0.5 \end{bmatrix}\), \(\det(J) = 0.5\) (the dimensionless natural-to-physical area-scaling factor).
Use 2-point Gauss quadrature to evaluate \(\int_{-1}^{1} (3x^2 + 2x + 1)dx\). Compare with the exact analytical result. The Gauss points are ξ₁ = -0.577, ξ₂ = 0.577, both with weight w = 1.
Show Solution
Gauss quadrature approximation:
\[I \approx w_1 f(\xi_1) + w_2 f(\xi_2)\]
\[f(\xi_1) = 3(-0.577)^2 + 2(-0.577) + 1 = 0.9988 - 1.154 + 1 = 0.8448\]
\[f(\xi_2) = 3(0.577)^2 + 2(0.577) + 1 = 0.9988 + 1.154 + 1 = 3.1528\]
\[I \approx 1 \times 0.8448 + 1 \times 3.1528 = 3.9976 \approx 4.0\]
Exact analytical result:
\[I = \int_{-1}^{1} (3x^2 + 2x + 1)dx = [x^3 + x^2 + x]_{-1}^{1}\]
\[I = (1 + 1 + 1) - (-1 + 1 - 1) = 3 - (-1) = 4\]
Error: |4 - 3.9976| / 4 ≈ 0.06% - essentially exact!
Explanation: 2-point Gauss quadrature is exact for polynomials up to degree 3. Our polynomial is degree 2, so the result is exact (the tiny residual error comes only from using the rounded Gauss points ±0.577 instead of the exact ±√(1/3) = ±0.57735...).
Answer: \(I \approx 3.9976\) versus the exact value \(4\); error \(\approx 0.06\%\).