9  Boundary Conditions in FEM: Theory and Practice

Boundary conditions are where the physics meets the algebra. This notebook covers how Neumann and Dirichlet conditions enter the FEM formulation.

9.1 How Boundary Conditions Arise

When we apply integration by parts to a flux-divergence term, we get a boundary integral:

\[\int_\Omega \frac{d(\text{flux})}{dx} \cdot v \, dx = -\int_\Omega \text{flux} \cdot \frac{dv}{dx} \, dx + \left[\text{flux} \cdot v\right]_{\text{boundary}}\]

This boundary term is the mathematical origin of both boundary condition types:

  • Neumann (natural) BCs: The flux itself is prescribed on part of the boundary
  • Dirichlet (essential) BCs: The solution value is prescribed on part of the boundary

The key insight: Neumann conditions enter the weak form naturally through the boundary integral, while Dirichlet conditions must be enforced by constraining the solution space.

9.2 Neumann (Natural) Boundary Conditions

Why “natural”? When flux multiplied by n is prescribed on Gamma_N, the boundary integral becomes a known contribution:

\[\int_{\Gamma_N} t^* \, v \, ds\]

This term is simply added to the load vector. No special handling of the stiffness matrix is needed—the condition is automatically incorporated when we assemble the system.

import sympy as sp
from IPython.display import Math, display
from symbolic_fem_workbench import (
    make_field_1d,
    integrate_divergence_1d,
    drop_dirichlet_boundary,
    apply_neumann_flux,
)

# Set up a simple 1D problem
x, L = sp.symbols("x L", positive=True)
E, A = sp.symbols("E A", positive=True)
P = sp.Symbol("P", real=True)  # External load
u = make_field_1d("u", x)
v = make_field_1d("v", x)
domain = (x, 0, L)

# Define the flux (stress times area)
flux = E * A * sp.diff(u, x)

# Apply integration by parts
domain_term, (left_bc, right_bc) = integrate_divergence_1d(flux, v, x, domain)

print("Boundary terms from integration by parts:")
print(f"Left boundary (x=0): {left_bc.evaluate()}")
print(f"Right boundary (x=L): {right_bc.evaluate()}")
print()
print("These boundary terms will be handled by Neumann and Dirichlet BCs.")
Boundary terms from integration by parts:
Left boundary (x=0): -A*E*v(0)*Subs(Derivative(u(x), x), x, 0)
Right boundary (x=L): -A*E*v(L)*Derivative(u(L), L)

These boundary terms will be handled by Neumann and Dirichlet BCs.
# Apply Neumann BC at x=L: prescribed force P
print("Applying Neumann BC: flux = P at x=L")
neumann_contrib = apply_neumann_flux(right_bc, flux_expr=flux, prescribed_flux=P)
print(f"Neumann contribution to load vector: {neumann_contrib}")
print()
print("Key insight: This term enters F(v) directly—no matrix modification needed.")
print("This is why Neumann BCs are called 'natural': they are built into the weak form.")
Applying Neumann BC: flux = P at x=L
Neumann contribution to load vector: -P*v(L)

Key insight: This term enters F(v) directly—no matrix modification needed.
This is why Neumann BCs are called 'natural': they are built into the weak form.

9.3 Dirichlet (Essential) Boundary Conditions

Why “essential”? Dirichlet BCs constrain solution values directly. Unlike Neumann BCs, they don’t enter through the weak form naturally. Instead, they must be enforced on the solution space and then on the algebraic system.

When u = g on Gamma_D, we restrict the test space to functions vanishing on Gamma_D:

\[V_D = \{v : v = 0 \text{ on } \Gamma_D\}\]

After discretization and assembly, we have a system K u = F with all DOFs. To enforce the Dirichlet condition u_i = g_i, we modify the system algebraically. Two main approaches:

  1. Row Substitution: Textbook approach, simple and explicit
  2. Lifting: Modern approach (FEniCS, deal.II), reveals mathematical structure

9.3.1 Approach 1: Row Substitution

The algorithm: For each constrained DOF i with prescribed value g_i:

  1. Subtract K_{:,i} times g_i from the RHS: F becomes F - K_{:,i} times g_i
  2. Zero out row i and column i of K
  3. Set K_{ii} = 1
  4. Set F_i = g_i

After these modifications, solving K u = F automatically gives u_i = g_i.

Pros: Simple, preserves matrix size, explicit
Cons: Destroys symmetry of K (temporarily), modifies matrix structure

import numpy as np
from symbolic_fem_workbench import (
    build_bar_1d_local_problem,
)
from symbolic_fem_workbench.assembly import (
    assemble_dense_matrix,
    assemble_dense_vector,
    apply_dirichlet_by_row_substitution,
)

# Build a 1D bar problem with 3 elements
problem = build_bar_1d_local_problem()
ke_fn = sp.lambdify((problem["L"], problem["E"], problem["A"]), problem["Ke"], "numpy")
fe_fn = sp.lambdify((problem["L"], problem["q"], problem["P"]), problem["fe"], "numpy")

n_elem, L_val, E_val, A_val, q_val = 3, 1.0, 1.0, 1.0, 1.0
h = L_val / n_elem
ndofs = n_elem + 1
K = np.zeros((ndofs, ndofs))
F = np.zeros(ndofs)

# Assemble the system
for e in range(n_elem):
    Ke = np.array(ke_fn(h, E_val, A_val), dtype=float)
    fe = np.array(fe_fn(h, q_val, 0.0), dtype=float).flatten()
    conn = [e, e + 1]
    assemble_dense_matrix(K, Ke, conn)
    assemble_dense_vector(F, fe, conn)

print("Assembled K (before BC enforcement):")
print(np.round(K, 4))
print()
print("Assembled F (before BC enforcement):", np.round(F, 4))
Assembled K (before BC enforcement):
[[ 3. -3.  0.  0.]
 [-3.  6. -3.  0.]
 [ 0. -3.  6. -3.]
 [ 0.  0. -3.  3.]]

Assembled F (before BC enforcement): [0.1667 0.3333 0.3333 0.1667]
# Apply Dirichlet BCs: u(0) = 0 and u(L) = 0 (fixed at both ends)
K_mod, F_mod = apply_dirichlet_by_row_substitution(
    K.copy(), F.copy(), constrained_dofs=[0, ndofs-1], prescribed_values=[0.0, 0.0]
)

print("After row substitution (u(0)=0, u(L)=0):")
print()
print("K_mod:")
print(np.round(K_mod, 4))
print()
print("F_mod:", np.round(F_mod, 4))
print()
print("Note: Rows/columns 0 and 3 are now identity rows")
print("      RHS entries 0 and 3 are now the prescribed values (0).")
After row substitution (u(0)=0, u(L)=0):

K_mod:
[[ 1.  0.  0.  0.]
 [ 0.  6. -3.  0.]
 [ 0. -3.  6.  0.]
 [ 0.  0.  0.  1.]]

F_mod: [0.     0.3333 0.3333 0.    ]

Note: Rows/columns 0 and 3 are now identity rows
      RHS entries 0 and 3 are now the prescribed values (0).
# Solve the modified system
u = np.linalg.solve(K_mod, F_mod)
print("Solution (row substitution):")
print(np.round(u, 6))
print()
print(f"Verification: u[0] = {u[0]:.6f}, u[3] = {u[3]:.6f}")
print("Both are 0.0 as prescribed.")
Solution (row substitution):
[0.       0.111111 0.111111 0.      ]

Verification: u[0] = 0.000000, u[3] = 0.000000
Both are 0.0 as prescribed.

9.3.2 Approach 2: Lifting (Modern FEniCS Approach)

The mathematical idea:

Decompose the solution as: u = u_0 + u_D

where: - u_D is a “lift” that satisfies the Dirichlet BCs: u_D = g on Gamma_D, zero elsewhere - u_0 is the homogeneous part: u_0 = 0 on Gamma_D

Substituting into the weak form a(u, v) = F(v): a(u_0 + u_D, v) = F(v) a(u_0, v) = F(v) - a(u_D, v)

In matrix form: K_{ff} u_{0,f} = F_f - K_{fc} g_c

Why is this better? 1. The test space has homogeneous BCs—no special handling needed 2. The Dirichlet data enters as a correction to the RHS, not as matrix surgery 3. Generalises naturally to non-homogeneous BCs and higher dimensions 4. This is how modern FEM frameworks (FEniCS, deal.II) handle Dirichlet conditions internally

from symbolic_fem_workbench.assembly import apply_dirichlet_by_lifting, expand_reduced_solution

# Re-assemble (since row substitution modified in place)
K2 = np.zeros((ndofs, ndofs))
F2 = np.zeros(ndofs)
for e in range(n_elem):
    Ke = np.array(ke_fn(h, E_val, A_val), dtype=float)
    fe = np.array(fe_fn(h, q_val, 0.0), dtype=float).flatten()
    assemble_dense_matrix(K2, Ke, [e, e+1])
    assemble_dense_vector(F2, fe, [e, e+1])

# Apply lifting
K_ff, F_f, free_dofs, u_lift = apply_dirichlet_by_lifting(
    K2, F2, constrained_dofs=[0, ndofs-1], prescribed_values=[0.0, 0.0]
)

print("Lifting decomposition:")
print(f"Lift vector u_D: {u_lift}")
print(f"Free DOFs: {free_dofs}")
print()
print("Modified RHS F_f:")
print(np.round(F_f, 4))
print()
print("Reduced stiffness K_ff (free-free block):")
print(np.round(K_ff, 4))
Lifting decomposition:
Lift vector u_D: [0. 0. 0. 0.]
Free DOFs: [1, 2]

Modified RHS F_f:
[0.3333 0.3333]

Reduced stiffness K_ff (free-free block):
[[ 6. -3.]
 [-3.  6.]]
# Solve for the homogeneous part u_0
u0_free = np.linalg.solve(K_ff, F_f)
print(f"Homogeneous solution u_0 (at free DOFs): {np.round(u0_free, 6)}")
print()

# Expand back to full size
u0_full = expand_reduced_solution(u0_free, ndofs, free_dofs, [0, ndofs-1], [0.0, 0.0])
print(f"u_0 (full size): {np.round(u0_full, 6)}")
print()

# Total solution = u_0 + u_D
u_total = u0_full + u_lift
print(f"u = u_0 + u_D: {np.round(u_total, 6)}")
Homogeneous solution u_0 (at free DOFs): [0.111111 0.111111]

u_0 (full size): [0.       0.111111 0.111111 0.      ]

u = u_0 + u_D: [0.       0.111111 0.111111 0.      ]

9.4 Comparison: Both Approaches Give the Same Answer

The mathematical principle behind both methods is the same—they are just different algebraic implementations. Let us verify:

print("Row substitution result:", np.round(u, 6))
print("Lifting result:         ", np.round(u_total, 6))
print()
print("Are they equal?", np.allclose(u, u_total))
print("Max difference:", np.max(np.abs(u - u_total)))
Row substitution result: [0.       0.111111 0.111111 0.      ]
Lifting result:          [0.       0.111111 0.111111 0.      ]

Are they equal? True
Max difference: 0.0

9.5 Non-Homogeneous Dirichlet Example

Now let us try a non-homogeneous case: u(0) = 0 and u(L) = 1

# Re-assemble
K3 = np.zeros((ndofs, ndofs))
F3 = np.zeros(ndofs)
for e in range(n_elem):
    Ke = np.array(ke_fn(h, E_val, A_val), dtype=float)
    fe = np.array(fe_fn(h, q_val, 0.0), dtype=float).flatten()
    assemble_dense_matrix(K3, Ke, [e, e+1])
    assemble_dense_vector(F3, fe, [e, e+1])

print("Non-homogeneous Dirichlet: u(0)=0, u(L)=1")
print()

# Row substitution
K3_mod, F3_mod = apply_dirichlet_by_row_substitution(
    K3.copy(), F3.copy(), constrained_dofs=[0, ndofs-1], prescribed_values=[0.0, 1.0]
)
u_row = np.linalg.solve(K3_mod, F3_mod)
print("Row substitution:", np.round(u_row, 6))

# Lifting
K3_ff, F3_f, free_dofs, u_lift3 = apply_dirichlet_by_lifting(
    K3, F3, constrained_dofs=[0, ndofs-1], prescribed_values=[0.0, 1.0]
)
u0_free = np.linalg.solve(K3_ff, F3_f)
# IMPORTANT: use zeros here because u_lift already carries the prescribed values
u0_full = expand_reduced_solution(u0_free, ndofs, free_dofs, [0, ndofs-1], [0.0, 0.0])
u_lift_total = u0_full + u_lift3
print("Lifting:         ", np.round(u_lift_total, 6))
print()
print("Match:", np.allclose(u_row, u_lift_total))
Non-homogeneous Dirichlet: u(0)=0, u(L)=1

Row substitution: [0.       0.444444 0.777778 1.      ]
Lifting:          [0.       0.444444 0.777778 1.      ]

Match: True

9.6 Exercise

Modify the problem above to apply a different boundary condition at the right end: - Try u(L) = 2 - Try u(0) = 1, u(L) = 0 - Try only fixing the left end: u(0) = 0 (what happens at the right?)

For each case, verify that both approaches still give the same answer.

# Your exercise code here