3  Galerkin Discretization: From Weak Form to Element Matrices

This notebook demonstrates the Galerkin discretization process, where we substitute finite element shape functions into the weak form to derive the element stiffness matrix and load vector. This is where abstract weak forms become concrete linear algebra.

import sympy as sp
from IPython.display import Math, display
from symbolic_fem_workbench import (
    make_field_1d, LinearElement1D,
    local_trial_expansion, local_test_expansion,
    substitute_fe, extract_coefficient_matrix, extract_coefficient_vector,
    build_bar_1d_local_problem,
)

3.1 Step 1: Recall the weak form

From the previous notebook, the weak form for the 1D bar is: \[a(u, v) = \int_0^L EA \frac{du}{dx} \frac{dv}{dx} \, dx\] \[F(v) = -\int_0^L q v \, dx + P \cdot v(L)\]

In the Galerkin method, we substitute \(u \approx u_h = \sum_i N_i d_i\) and \(v = N_j\) to get a finite-dimensional problem.

problem = build_bar_1d_local_problem()
x, L = problem["x"], problem["L"]
u, v = problem["u"], problem["v"]
E, A = problem["E"], problem["A"]

print("Bilinear form a(u,v):")
display(Math(sp.latex(problem["weak_bilinear"])))
print("\nLinear form F(v):")
display(Math(sp.latex(problem["weak_linear"])))

Bilinear form a(u,v):

\(\displaystyle A E \int\limits_{0}^{L} \frac{d}{d x} u{\left(x \right)} \frac{d}{d x} v{\left(x \right)}\, dx\)

Linear form F(v):

\(\displaystyle q \int\limits_{0}^{L} v{\left(x \right)}\, dx\)

3.2 Step 2: Define shape functions on an element

For a linear 1D element on the interval \([0, L]\), we use two linear shape functions (P1 element): \[N_1(x) = 1 - \frac{x}{L}, \quad N_2(x) = \frac{x}{L}\]

These form a partition of unity: \(N_1 + N_2 = 1\)

element = LinearElement1D(x=x, x0=0, x1=L)
N = element.shape_functions

print(f"Element length: {element.length}")
print(f"Shape function N1 = {N[0]}")
print(f"Shape function N2 = {N[1]}")
print(f"Sum (partition of unity): {sp.simplify(N[0] + N[1])}")

Element length: L
Shape function N1 = (L - x)/L
Shape function N2 = x/L
Sum (partition of unity): 1

3.3 Step 3: Build trial and test expansions

The Galerkin discretization replaces the infinite-dimensional spaces with finite-dimensional ones: \[u_h(x) = N_1(x) d_0 + N_2(x) d_1\] \[v_h(x) = N_1(x) w_0 + N_2(x) w_1\]

where \(d_0, d_1\) are nodal values and \(w_0, w_1\) are test function coefficients.

d0, d1 = sp.symbols("d0 d1", real=True)
w0, w1 = sp.symbols("w0 w1", real=True)

uh = local_trial_expansion(N, [d0, d1])
vh = local_test_expansion(N, [w0, w1])

print("u_h =", uh)
print("v_h =", vh)

u_h = d0 - d0*x/L + d1*x/L
v_h = w0 - w0*x/L + w1*x/L

3.4 Step 4: Substitute into the weak form

The Galerkin discretization consists of substituting \(u \to u_h\) and \(v \to v_h\) into the weak form integrals. This produces a bilinear form that depends on coefficients \(d_i\) and \(w_j\).

a_disc = substitute_fe(problem["weak_bilinear"], {u: uh, v: vh}, x)
l_disc = substitute_fe(problem["weak_linear"], {v: vh}, x)

print("Discretized bilinear form:")
display(Math(sp.latex(a_disc)))
print("\nDiscretized linear form:")
display(Math(sp.latex(l_disc)))

Discretized bilinear form:

\(\displaystyle \frac{A E d_{0} w_{0}}{L} - \frac{A E d_{0} w_{1}}{L} - \frac{A E d_{1} w_{0}}{L} + \frac{A E d_{1} w_{1}}{L}\)

Discretized linear form:

\(\displaystyle \frac{L q w_{0}}{2} + \frac{L q w_{1}}{2}\)

3.5 Step 5: Extract element matrices

The key insight: The discretized bilinear form is quadratic in coefficients. By comparing coefficients of \(d_i w_j\), we extract the element stiffness matrix \(K_e\). Similarly, coefficients of \(w_j\) give the element load vector \(f_e\).

Ke = extract_coefficient_matrix(a_disc, [d0, d1], [w0, w1])
fe = extract_coefficient_vector(l_disc, [w0, w1])

print("Element stiffness matrix Ke:")
display(Math(r"K_e = " + sp.latex(sp.simplify(Ke))))
print("\nElement load vector fe:")
display(Math(r"f_e = " + sp.latex(sp.simplify(fe))))

Element stiffness matrix Ke:

\(\displaystyle K_e = \left[\begin{matrix}\frac{A E}{L} & - \frac{A E}{L}\\- \frac{A E}{L} & \frac{A E}{L}\end{matrix}\right]\)

Element load vector fe:

\(\displaystyle f_e = \left[\begin{matrix}\frac{L q}{2}\\\frac{L q}{2}\end{matrix}\right]\)

3.6 Verify: Use the pre-built workflow

The symbolic_fem_workbench provides pre-built workflows that perform all these steps at once. Let’s verify our manual derivation matches the library.

print("Pre-built workflow gives the same result:")
display(Math(r"K_e = " + sp.latex(problem["Ke"])))
display(Math(r"f_e = " + sp.latex(problem["fe"])))

Pre-built workflow gives the same result:

\(\displaystyle K_e = \left[\begin{matrix}\frac{A E}{L} & - \frac{A E}{L}\\- \frac{A E}{L} & \frac{A E}{L}\end{matrix}\right]\)

\(\displaystyle f_e = \left[\begin{matrix}\frac{L q}{2}\\\frac{L q}{2}\end{matrix}\right]\)

3.7 Exercise

Question: What happens if you use a 3-node (P2) element instead? Define shape functions for a quadratic element on \([0, L]\): \[N_1 = (1 - \xi)(1 - 2\xi), \quad N_2 = 4\xi(1-\xi), \quad N_3 = \xi(2\xi - 1)\] where \(\xi = x/L\).

Substitute and extract the \(3 \times 3\) element stiffness matrix. How do the entries compare to the P1 case?

# Exercise: What happens if you use a 3-node (P2) element instead?
# Define shape functions for a quadratic element on [0, L]:
#   N1 = (1 - xi)(1 - 2*xi),  N2 = 4*xi*(1-xi),  N3 = xi*(2*xi - 1)
# where xi = x/L.
# Substitute and extract the 3x3 element stiffness matrix.