14) FEM Assembly

14) FEM Assembly#

Last time:

Introduction to the Finite Element Method
Variational formulations
Galerkin Spectral Element Method (SEM): \(L^2\) projection

Today:

Linear Finite Elements
Galerkin method for 1D Poisson equation
Abstractions

1. Linear Finite Elements#

The domain \(\Omega = [0,1]\) is discretized with elements (subintervals \(\Omega_h\))

\[ 0 = x_0 < x_1 < \cdots < x_N = 1 \]

such that \(x_i = i h \), \(i=0,1,\ldots, n\), where

\[ h_i = x_i - x_{i-1} \]

is the length of the subinterval \((x_{i-1}, x_i)\), \(i=1, 2, \ldots, n\).

Linear (Hat) Basis Functions#

\[\begin{split} \phi_i(x) = \begin{cases} \frac{x - x_{i-1}}{h_i}, & x \in [x_{i-1}, x_i] \\[8pt] \frac{x_{i+1} - x}{h_{i+1}}, & x \in [x_i, x_{i+1}] \\[8pt] 0, & \text{otherwise} \end{cases} \end{split}\]

Nodal Property#

\[ \phi_i(x_j) = \delta_{ij} \]

where \(\delta_{ij}\) denotes the Kronecker delta. This just says that \(\phi_i(x)\) has zeros at all the \(x_j\)’s (with \(j \neq i\)) except \(x_i\), where it is identically equal to \(1\).

Compact support#

\[ \text{supp}(\phi_i) = [x_{i-1}, x_{i+1}] \]

Derivatives#

\[\begin{split} \phi_i'(x) = \begin{cases} \frac{1}{h_i}, & x \in [x_{i-1}, x_i] \\[8pt] -\frac{1}{h_{i+1}}, & x \in [x_i, x_{i+1}] \\[8pt] 0, & \text{otherwise} \end{cases} \end{split}\]

Uniform Mesh Case#

When a mesh is partitioned in uniform (all identically equal) elements, we have that \(h_i = h\) and therefore, we can rewrite

\[\begin{split} \phi_i(x) = \begin{cases} \frac{x - x_{i-1}}{h}, & x \in [x_{i-1}, x_i] \\[8pt] \frac{x_{i+1} - x}{h}, & x \in [x_i, x_{i+1}] \\[8pt] 0, & \text{otherwise} \end{cases} \end{split}\]

\[\begin{split} \phi_i'(x) = \begin{cases} \frac{1}{h}, & x \in [x_{i-1}, x_i] \\[8pt] -\frac{1}{h}, & x \in [x_i, x_{i+1}] \end{cases} \end{split}\]

using Plots
using LaTeXStrings
default(linewidth=4, legendfontsize=12)

# Piecewise linear hat function centered at node xc,
# with left neighbor xl and right neighbor xr
function hat(x, xl, xc, xr)
    if x < xl || x > xr
        0.0
    elseif x <= xc
        (x - xl) / (xc - xl)
    else
        (xr - x) / (xr - xc)
    end
end

xim1 = -1.0
xi   =  0.0
xip1 =  1.0

x = range(xim1 - 0.5, xip1 + 0.5, length=500)
phi_i = [hat(xx, xim1, xi, xip1) for xx in x]

plot(x, phi_i,
    lw=3,
    label=L"\phi_i(x)",
    xticks=([xim1, xi, xip1], [L"x_{i-1}", L"x_i", L"x_{i+1}"]),
    yticks=([0, 1], ["0", "1"]),
    ylims=(-0.05, 1.1),
    framestyle=:origin
)

scatter!([xim1, xi, xip1], [0, 1, 0], label="")

# Example node locations
xim2 = -2.0
xim1 = -1.0
xi   =  0.0
xip1 =  1.0
xip2 =  2.0

# Plot grid
x = range(xim2, xip2, length=1000)

# Three hat functions
phi_im1 = [hat(xx, xim2, xim1, xi)   for xx in x]
phi_i   = [hat(xx, xim1, xi,   xip1) for xx in x]
phi_ip1 = [hat(xx, xi,   xip1, xip2) for xx in x]

# Base plot
plot(x, phi_im1,
    lw=3,
    label=L"\phi_{i-1}(x)",
    xlabel=L"x",
    ylabel="",
    legend=:top,
    ylims=(-0.05, 1.2),
    xlims=(xim2 - 0.4, xip2 + 0.5),
    xticks=([xim2, xim1, xi, xip1, xip2],
            [L"x_{i-2}", L"x_{i-1}", L"x_i", L"x_{i+1}", L"x_{i+2}"]),
    yticks=[],
    framestyle=:origin
)

plot!(x, phi_i,   lw=3, label=L"\phi_i(x)")
plot!(x, phi_ip1, lw=3, label=L"\phi_{i+1}(x)")

# Mark the nodes
scatter!([xim2, xim1, xi, xip1, xip2], zeros(5),
    label="",
    markersize=5
)

# Optional vertical guide lines at the peaks
plot!([xim1, xim1], [0, 1], lw=1, ls=:dash, color=:black, label="")
plot!([xi,   xi],   [0, 1], lw=1, ls=:dash, color=:black, label="")
plot!([xip1, xip1], [0, 1], lw=1, ls=:dash, color=:black, label="", legend =:topright)

2. Galerkin method for 1D Poisson equation#

For the 1D two-point BVP for the Poisson equation:

\[ - u_{xx}(x) = f(x), \qquad 0 \leq x \leq 1, \qquad u(0)=0, \; u(1) = 0 \, . \]

We have seen that we start by constructing a weak form or variational formulation of the problem

Step 1: Multiply both sides of the PDE by a test function \(v(x)\) with compact support (i.e., \( v(0)=0, \; v(1) = 0 \)).

\[ - u'' v = fv \]

Step 2: Integrate both sides of the equation across the domain \(\Omega = [0,1]\) and apply integrations by parts. On the left-hand-side we have:

\[ \int_0^1 (-u'' v) dx = \left[ -u' v \right]^1_0 + \int_0^1 u' v' dx = \int_0^1 u' v' dx \]

We note that the boundary term vanishes because \(v\) has compact support.

Hence the weak form of the PDE is:

\[ \int_0^1 u' v' dx = \int_0^1 fv dx \]

Step 3: Construct the set of linear (hat) basis function as above

\[\begin{split} \phi_i(x) = \begin{cases} \frac{x - x_{i-1}}{h}, & x \in [x_{i-1}, x_i] \\[8pt] \frac{x_{i+1} - x}{h}, & x \in [x_i, x_{i+1}] \\[8pt] 0, & \text{otherwise} \end{cases} \end{split}\]

\[\begin{split} \phi_i'(x) = \begin{cases} \frac{1}{h}, & x \in [x_{i-1}, x_i] \\[8pt] -\frac{1}{h}, & x \in [x_i, x_{i+1}] \end{cases} \end{split}\]

Step 4: Expand all functions in terms of basis functions, such as

\[\begin{split} \begin{align} u_h(\mathbf x) &= \sum_{j=1}^n u_j \, \phi_j(\mathbf x) ,\\ u'_h(\mathbf x) &= \sum_{j=1}^n u_j \, \phi'_j(\mathbf x) \, . \end{align} \end{split}\]

where the coefficients \(\{u_j\}\) are the unknown nodal values to be determined.

On assuming the hat basis functions, obviously \(u_h(x)\) is also a piecewise linear function, although this is not usually the case for the true solution \(u(x)\).

We then derive a linear system of equations for the coefficients by substituting the approximate solution \(u_h(x)\) for the exact solution \(u(x)\) in the weak form

\[ \int_0^1 u' v' dx = \int_0^1 fv dx\]

Hence,

\[ \int_0^1 \sum_{i=1}^{n} u_i \phi'_iv' dx = \sum_{i=1}^{n} u_i \int_0^1 \phi_i'v' dx = \int_0^1 fv dx \, . \]

Hence, we get a system of equations

\[\begin{split} \begin{align*} \left( \int_0^1 \phi_1' \phi_1' dx \right) u_1 + \ldots + \left( \int_0^1 \phi_1' \phi_n' dx \right) u_n &= \int_0^1 f \phi_1 dx \\ \left( \int_0^1 \phi_2' \phi_1' dx \right) u_1 + \ldots + \left( \int_0^1 \phi_2' \phi_n' dx \right) u_n &= \int_0^1 f \phi_2 dx \\ \ldots & \\ \left( \int_0^1 \phi_i' \phi_1' dx \right) u_1 + \ldots + \left( \int_0^1 \phi_i' \phi_n' dx \right) u_n &= \int_0^1 f \phi_i dx \\ \ldots & \\ \left( \int_0^1 \phi_{n-1}' \phi_1' dx \right) u_1 + \ldots + \left( \int_0^1 \phi_{n-1}' \phi_n' dx \right) u_n &= \int_0^1 f\phi_{n-1} dx \\ \left( \int_0^1 \phi_n' \phi_1' dx \right) u_1 + \ldots + \left( \int_0^1 \phi_n' \phi_n' dx \right) u_n &= \int_0^1 f \phi_n dx \end{align*} \end{split}\]

Or in matrix form:

\[\begin{split} \begin{align*} \left[ \begin{array}{cccc} a(\phi_1, \phi_1) & a(\phi_1, \phi_2) & \ldots & a(\phi_1, \phi_n) \\ a(\phi_2, \phi_1) & a(\phi_2, \phi_2) & \ldots & a(\phi_2, \phi_n) \\ \vdots & \vdots & \vdots & \vdots \\ a(\phi_{n-1}, \phi_1) & a(\phi_{n-1}, \phi_2) & \ldots & a(\phi_{n-1}, \phi_n)\\ a(\phi_n, \phi_1) & a(\phi_n, \phi_2) & \ldots & a(\phi_n, \phi_n) \end{array} \right] \left[ \begin{array}{c} u_1 \\ u_2\\ \vdots\\ u_{n-1} \\ u_n \end{array} \right] = \left[ \begin{array}{c} \langle f,\phi_1\rangle_h\\ \langle f,\phi_2\rangle_h\\ \vdots \\ \langle f,\phi_{n-1}\rangle_h\\ \langle f,\phi_n\rangle_h \end{array} \right] \,. \end{align*} \end{split}\]

Where we have used the bilinear form

\[\begin{split} a(\phi_i,\phi_j) = \int_{0}^1 \nabla \phi_i \cdot \nabla \phi_j ~\mathrm{d} x = \begin{cases} \frac{2}{h} &, i = j \\ \frac{-1}{h} &, |i-j| = 1 \\ 0 &, |i-j| \geq 2 \end{cases} \end{split}\]

and the inner product

\[ \langle f, \phi_i \rangle = \int_0^1 f \phi_i ~\mathrm{d} x \,. \]

The matrix problem above can be written more compactly as

\[ K \mathbf{u} = F \]

where the matrix \(K\) is called the Stiffness Matrix.

Therefore, for this case we get the matrix

\[\begin{split} \begin{align*} \frac{1}{h} \left[ \begin{array}{cccccc} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & & \ddots & \ddots & \ddots & \\ & & & -1 & 2 & -1 \\ & & & & -1 & 2 \end{array} \right] \left[ \begin{array}{c} u_1 \\ u_2\\ \vdots\\ u_{n-1} \\ u_n \end{array} \right] = \left[ \begin{array}{c} \langle f,\phi_1\rangle_h\\ \langle f,\phi_2\rangle_h\\ \vdots \\ \langle f,\phi_{n-1}\rangle_h\\ \langle f,\phi_n\rangle_h \end{array} \right] \,. \end{align*} \end{split}\]

The right-hand side can be explicitly written as

\[ \langle f,\phi_{i}\rangle_h = \sum_{j=0}^n \frac{1}{2} h \left[ f(x_j) \phi_i (x_j) + f(x_{j+1}) \phi_{i} (x_{j+1})\right] = f_i h \]

where we have used the simple trapezoidal rule for the \(L^2\) inner product integral.

Hence, we can rewrite our matrix problem as

\[\begin{split} \begin{align*} \frac{1}{h} \left[ \begin{array}{cccccc} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & & \ddots & \ddots & \ddots & \\ & & & -1 & 2 & -1 \\ & & & & -1 & 2 \end{array} \right] \left[ \begin{array}{c} u_1 \\ u_2\\ \vdots\\ u_{n-1} \\ u_n \end{array} \right] = h \left[ \begin{array}{c} f_1\\ f_2\\ \vdots \\ f_{n-1}\\ f_n \end{array} \right] \,. \tag{1} \end{align*} \end{split}\]

Step 5: Solve the linear system of equations for the coefficients \( \{u_j\}\) and hence obtain the approximate solution
\[ u_h(\mathbf x) = \sum_{j=1}^n u_j \, \phi_j(\mathbf x) \]

Observations#

Does the matrix in (1) look familiar??
We note that for solving \(-u'' = f(x)\), \(u(0) = u(1) = 0\), the matrix form of the \(P^1\) linear Finite Element Method with trapezoidal quadrature is precisely the second order centered Finite Difference Scheme.

3. Abstractions#

In general, In the FEM/SEM fomrulations, we always try to keep the order of derivatives of \(u\) and \(v\) as small as possible. We reach a general form, where we can express as \(f_0\) all terms in the weak form that multiply the test function \(v\) and \(f_1\) all terms that multiply its gradient \(\nabla v\):

General form:

\[ \int_\Omega v\cdot f_0(u, \nabla u) + \nabla v\cdot f_1(u, \nabla u) = 0, \quad \forall v\]

Many FEM/SEM libraries (like libCEED) discretize this as

\[ \sum_e \mathcal E_e^T \Big( B^T W \left\lvert \frac{\partial x}{\partial X} \right\rvert f_0(\tilde u, \nabla \tilde u) + D^T \left(\frac{\partial X}{\partial x}\right)^{T} W \left\lvert \frac{\partial x}{\partial X} \right\rvert f_1(\tilde u, \nabla\tilde u) \Big) = 0 \]

where \(\tilde u = B \mathcal E_e u\) and \(\nabla \tilde u = \frac{\partial X}{\partial x} D \mathcal E_e u\) are the values and gradients evaluated at quadrature points.

libCEED’s operator decomposition is summarized in the following graphic.

libCEED operator decomposition

Isoparametric mapping#

Given the reference coordinates \(X \in \hat K \subset R^n\) and physical coordinates \(x(X)\), an integral on the physical element can be written

\[ \int_{K = x(\hat K)} f(x) dx = \int_K \underbrace{\left\lvert \frac{\partial x}{\partial X} \right\rvert}_{\text{Jacobian determinant}} f(x(X)) dX .\]

Notation	Meaning
\(x\)	physical coordinates
\(X\)	reference coordinates
\(\mathcal E_e\)	restriction from global vector to element \(e\)
\(B\)	values of nodal basis functions at quadrature ponits on reference element
\(D\)	gradients of nodal basis functions at quadrature points on reference element
\(W\)	diagonal matrix of quadrature weights on reference element
\(\frac{\partial x}{\partial X} = D \mathcal E_e x \)	gradient of physical coordinates with respect to reference coordinates
\(\left\lvert \frac{\partial x}{\partial X}\right\rvert\)	determinant of coordinate transformation at each quadrature point
\(\frac{\partial X}{\partial x} = \left(\frac{\partial x}{\partial X}\right)^{-1}\)	derivative of reference coordinates with respect to physical coordinates

Finite element mesh and restriction#

using LinearAlgebra
using SparseArrays
using FastGaussQuadrature

function my_spy(A)
    cmax = norm(vec(A), Inf)
    s = max(1, ceil(120 / size(A, 1)))
    spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

function vander_legendre_deriv(x, k=nothing)
    if isnothing(k)
        k = length(x) # Square by default
    end
    m = length(x)
    Q = ones(m, k)
    dQ = zeros(m, k)
    Q[:, 2] = x
    dQ[:, 2] .= 1
    for n in 1:k-2
        Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
        dQ[:, n+2] = (2*n + 1) * Q[:,n+1] + dQ[:,n]
    end
    Q, dQ
end

function febasis(P, Q)
    x, _ = gausslobatto(P)
    q, w = gausslegendre(Q)
    Pk, _ = vander_legendre_deriv(x)
    Bp, Dp = vander_legendre_deriv(q, P)
    B = Bp / Pk
    D = Dp / Pk
    x, q, w, B, D
end

function L2_galerkin(P, Q, f)
    x, q, w, B, _ = febasis(P, Q)
    A = B' * diagm(w) * B
    rhs = B' * diagm(w) * f.(q)
    u = A \ rhs
    x, u
end

L2_galerkin (generic function with 1 method)

function fe1_mesh(P, nelem)
    x = LinRange(-1, 1, nelem+1)
    rows = Int[]
    cols = Int[]
    for i in 1:nelem
        append!(rows, (i-1)*P+1:i*P)
        append!(cols, (i-1)*(P-1)+1:i*(P-1)+1)
    end
    x, sparse(cols, rows, ones(nelem*P))'
end
P, nelem = 4, 3
x, E = fe1_mesh(P, nelem)
E

12×10 Adjoint{Float64, SparseMatrixCSC{Float64, Int64}} with 12 stored entries:
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0

function xnodal(x, P)
    xn = Float64[]
    xref, _ = gausslobatto(P)
    for i in 1:length(x)-1
        xL, xR = x[i:i+1]
        append!(xn, (xL+xR)/2 .+ (xR-xL)/2 * xref[1+(i>1):end])
    end
    xn
end
xn = xnodal(x, 10)
@show xn
scatter(xn, zero, marker=:square)

xn = [-1.0, -0.9731779693888196, -0.912924621701835, -0.8259749832701482, -0.7217596525554624, -0.6115736807778711, -0.5073583500631853, -0.42040871163149846, -0.36015536394451386, -0.3333333333333335, -0.3065113027221531, -0.24625795503516845, -0.1593083166034816, -0.05509298588879577, 0.055092985888795604, 0.15930831660348144, 0.24625795503516829, 0.3065113027221529, 0.33333333333333326, 0.36015536394451364, 0.4204087116314983, 0.5073583500631851, 0.6115736807778709, 0.7217596525554624, 0.8259749832701482, 0.912924621701835, 0.9731779693888196, 1.0]

Finite Element Method (FEM) building blocks in Julia#

struct FESpace
    P::Int
    Q::Int
    nelem::Int
    x::Vector
    xn::Vector
    Et::SparseMatrixCSC{Float64, Int64}
    q::Vector
    w::Vector
    B::Matrix
    D::Matrix
    function FESpace(P, Q, nelem)
        x, E = fe1_mesh(P, nelem)
        xn = xnodal(x, P)
        _, q, w, B, D = febasis(P, Q)
        new(P, Q, nelem, x, xn, E', q, w, B, D)
    end
end

# Extract out what we need for element e
function fe_element(fe, e)
    xL, xR = fe.x[e:e+1]
    q = (xL+xR)/2 .+ (xR-xL)/2*fe.q
    w = (xR - xL)/2 * fe.w
    E = fe.Et[:, (e-1)*fe.P+1:e*fe.P]'
    dXdx = ones(fe.Q) * 2 / (xR - xL)
    q, w, E, dXdx
end

fe = FESpace(3, 3, 5)
q, w, E, dXdx = fe_element(fe, 1);
@show q
@show sum(w)
E

q = [-0.9549193338482966, -0.8, -0.6450806661517035]
sum(w) = 0.3999999999999999

3×11 Adjoint{Float64, SparseMatrixCSC{Float64, Int64}} with 3 stored entries:
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 

Finite element residual assembly#

\[ v^T F(u) \sim \int_\Omega v\cdot f_0(u, \nabla u) + \nabla v\cdot f_1(u, \nabla u) = 0, \quad \forall v\]