10  2D Boundary Conditions: Neumann Edges and Dirichlet Lifting

This notebook extends the 1D boundary condition concepts to 2D triangular meshes, showing how Neumann conditions are applied on edge boundaries and how Dirichlet conditions are enforced across a full problem.

import numpy as np
import sympy as sp
from IPython.display import Math, display

from symbolic_fem_workbench import (
    ReferenceTriangleP1,
    AffineTriangleMap2D,
    integrate_reference_triangle_exact,
    neumann_load_vector_triangle_edge,
)
from symbolic_fem_workbench.assembly import (
    assemble_dense_matrix,
    assemble_dense_vector,
    apply_dirichlet_by_row_substitution,
    apply_dirichlet_by_lifting,
    expand_reduced_solution,
)

print("Imports successful. Ready for 2D FEM with boundary conditions.")

Imports successful. Ready for 2D FEM with boundary conditions.

10.1 Neumann BCs on Triangle Edges

In 2D, the boundary integral becomes a line integral along each edge:

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

For a straight edge of length |e| in a P1 triangle, if the flux t-star is constant, the contribution to each endpoint node is:

\[f_i^{\text{edge}} = t^* \cdot \frac{|e|}{2}\]

This comes from integrating the constant flux against the linear shape function over the edge.

# Define reference coordinates first
xi, eta = sp.symbols('xi eta', real=True, positive=True)

# Create a reference P1 triangle
ref_tri = ReferenceTriangleP1(xi, eta)
print(f"Reference triangle shape functions:")
for idx, N in enumerate(ref_tri.shape_functions):
    print(f"  N_{idx}: {N}")

# Define a constant Neumann flux and edge length
t_star = sp.Symbol('t_star', real=True)  # Prescribed flux
edge_len = sp.Symbol('edge_len', positive=True)

# Compute Neumann load contribution for edge 0 (between nodes 0 and 1)
# Nodes 0 and 1 are at (0,0) and (1,0) on the reference triangle
f_edge = neumann_load_vector_triangle_edge(
    ref_tri,
    edge_local_nodes=(0, 1),
    prescribed_flux=t_star,
    edge_length=edge_len,
    xi=xi,
    eta=eta,
)

print(f"\nNeumann load vector for edge 0-1:")
print(f"f_edge = {f_edge.T}")
print()
print(f"Each endpoint gets t_star times edge_len divided by 2")
print(f"The opposite node (node 2) gets 0")

Reference triangle shape functions:
  N_0: -eta - xi + 1
  N_1: xi
  N_2: eta

Neumann load vector for edge 0-1:
f_edge = Matrix([[edge_len*t_star/2, edge_len*t_star/2, 0]])

Each endpoint gets t_star times edge_len divided by 2
The opposite node (node 2) gets 0

10.2 2D Problem Setup: 4-Triangle Mesh with Mixed BCs

Consider a simple domain made of 4 triangles in a 1x1 square:

3 -------- 4 -------- 5
|\\\\ tri2 /|\\\\ tri3 /|
| \\\\ /  |  \\\\ /  |
| tri0 \\/ tri1 |  (interior edges omitted)
|    /  \\  |    /  \\  |
|   /    \\ |   /    \\ |
|  /      \\|  /      \\|
0 -------- 1 -------- 2

Boundary Conditions: - Bottom (nodes 0, 1, 2): Dirichlet u = 0 (fixed) - Top (nodes 3, 4, 5): Neumann t-star = 1.0 (prescribed flux) - Left & Right: Natural (homogeneous Neumann, no load)

# Simple 4-triangle mesh for a 1x1 domain
# Nodes
nodes = np.array([
    [0.0, 0.0],  # 0: bottom-left
    [0.5, 0.0],  # 1: bottom-center
    [1.0, 0.0],  # 2: bottom-right
    [0.0, 1.0],  # 3: top-left
    [0.5, 1.0],  # 4: top-center
    [1.0, 1.0],  # 5: top-right
])
n_nodes = len(nodes)

# Connectivity (triangles)
# Each row: [node_a, node_b, node_c]
connectivity = np.array([
    [0, 1, 3],  # triangle 0 (bottom-left)
    [1, 4, 3],  # triangle 1 (top-left)
    [1, 2, 4],  # triangle 2 (bottom-right)
    [2, 5, 4],  # triangle 3 (top-right)
])
n_elem = len(connectivity)

print(f"Mesh: {n_nodes} nodes, {n_elem} elements")
print(f"\nNode coordinates:")
for i, (x, y) in enumerate(nodes):
    print(f"  Node {i}: ({x:.1f}, {y:.1f})")

print(f"\nElement connectivity:")
for e, (a, b, c) in enumerate(connectivity):
    print(f"  Element {e}: nodes {a}, {b}, {c}")

Mesh: 6 nodes, 4 elements

Node coordinates:
  Node 0: (0.0, 0.0)
  Node 1: (0.5, 0.0)
  Node 2: (1.0, 0.0)
  Node 3: (0.0, 1.0)
  Node 4: (0.5, 1.0)
  Node 5: (1.0, 1.0)

Element connectivity:
  Element 0: nodes 0, 1, 3
  Element 1: nodes 1, 4, 3
  Element 2: nodes 1, 2, 4
  Element 3: nodes 2, 5, 4

# For a teaching example, we use a uniform Poisson problem:
#   -nabla u = 0 in Omega
#   u = 0 on Gamma_D (bottom)
#   du/dn = t-star on Gamma_N (top)

# Local stiffness matrix for P1 triangle (standard Poisson)
# For a uniform reference triangle with area 0.5:
K_ref = np.array([
    [1, -1/2, -1/2],
    [-1/2, 1, -1/2],
    [-1/2, -1/2, 1],
]) / 2

# Assemble global system
K_global = np.zeros((n_nodes, n_nodes))
F_global = np.zeros(n_nodes)

# Assemble element matrices
for e, (a, b, c) in enumerate(connectivity):
    # Get physical coordinates
    x_a, y_a = nodes[a]
    x_b, y_b = nodes[b]
    x_c, y_c = nodes[c]
    
    # Compute physical element area (0.5 times |det(J)|)
    J = np.array([
        [x_b - x_a, x_c - x_a],
        [y_b - y_a, y_c - y_a],
    ])
    area = 0.5 * abs(np.linalg.det(J))
    
    # Scale local stiffness matrix by area
    K_e = K_ref / area
    
    # Assemble
    for i, global_i in enumerate([a, b, c]):
        for j, global_j in enumerate([a, b, c]):
            K_global[global_i, global_j] += K_e[i, j]

print(f"Assembled global K ({n_nodes}x{n_nodes}):")
print(np.round(K_global, 3))
print(f"\nF_global (before boundary conditions): {F_global}")

Assembled global K (6x6):
[[ 2. -1.  0. -1.  0.  0.]
 [-1.  6. -1. -2. -2.  0.]
 [ 0. -1.  4.  0. -2. -1.]
 [-1. -2.  0.  4. -1.  0.]
 [ 0. -2. -2. -1.  6. -1.]
 [ 0.  0. -1.  0. -1.  2.]]

F_global (before boundary conditions): [0. 0. 0. 0. 0. 0.]

10.3 Apply Neumann BC on Top Edge

The top boundary consists of edges between nodes (3,4) and (4,5). Each has length 0.5.

# Apply Neumann BC: t-star = 1.0 on top edges
t_star_val = 1.0
edge_len_val = 0.5  # Length of each top edge

# Top edge 1: nodes 3 and 4
f_edge_34 = t_star_val * edge_len_val / 2
F_global[3] += f_edge_34
F_global[4] += f_edge_34

# Top edge 2: nodes 4 and 5
f_edge_45 = t_star_val * edge_len_val / 2
F_global[4] += f_edge_45  # Node 4 appears in both edges!
F_global[5] += f_edge_45

print(f"After adding Neumann flux (t-star={t_star_val}) on top edges:")
print(f"F_global: {np.round(F_global, 4)}")
print()
print(f"Note: Node 4 (top-center) receives contributions from both edges.")

After adding Neumann flux (t-star=1.0) on top edges:
F_global: [0.   0.   0.   0.25 0.5  0.25]

Note: Node 4 (top-center) receives contributions from both edges.

10.4 Enforce Dirichlet BC on Bottom (Row Substitution)

# Enforce Dirichlet: u = 0 on bottom (nodes 0, 1, 2)
K_mod, F_mod = apply_dirichlet_by_row_substitution(
    K_global.copy(), F_global.copy(),
    constrained_dofs=[0, 1, 2],
    prescribed_values=[0.0, 0.0, 0.0]
)

print("After enforcing u=0 on bottom (row substitution):")
print(f"Modified K:")
print(np.round(K_mod, 3))
print(f"\nModified F: {np.round(F_mod, 4)}")

After enforcing u=0 on bottom (row substitution):
Modified K:
[[ 1.  0.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.]
 [ 0.  0.  0.  4. -1.  0.]
 [ 0.  0.  0. -1.  6. -1.]
 [ 0.  0.  0.  0. -1.  2.]]

Modified F: [0.   0.   0.   0.25 0.5  0.25]

# Solve the system
u_row = np.linalg.solve(K_mod, F_mod)

print("Solution (row substitution):")
for i in range(n_nodes):
    x, y = nodes[i]
    print(f"  u[{i}] = {u_row[i]:.6f}  (at ({x:.1f}, {y:.1f}))")

Solution (row substitution):
  u[0] = 0.000000  (at (0.0, 0.0))
  u[1] = 0.000000  (at (0.5, 0.0))
  u[2] = 0.000000  (at (1.0, 0.0))
  u[3] = 0.095238  (at (0.0, 1.0))
  u[4] = 0.130952  (at (0.5, 1.0))
  u[5] = 0.190476  (at (1.0, 1.0))

10.5 Enforce Same BCs Using Lifting

# Re-assemble (lifting does not modify in place)
K_global2 = np.zeros((n_nodes, n_nodes))
F_global2 = np.zeros(n_nodes)

for e, (a, b, c) in enumerate(connectivity):
    x_a, y_a = nodes[a]
    x_b, y_b = nodes[b]
    x_c, y_c = nodes[c]
    
    J = np.array([
        [x_b - x_a, x_c - x_a],
        [y_b - y_a, y_c - y_a],
    ])
    area = 0.5 * abs(np.linalg.det(J))
    K_e = K_ref / area
    
    for i, global_i in enumerate([a, b, c]):
        for j, global_j in enumerate([a, b, c]):
            K_global2[global_i, global_j] += K_e[i, j]

# Add Neumann contributions
F_global2[3] += f_edge_34
F_global2[4] += f_edge_34 + f_edge_45
F_global2[5] += f_edge_45

print("Re-assembled global system for lifting approach")

Re-assembled global system for lifting approach

# Apply lifting: decompose u = u_0 + u_D
K_ff, F_f, free_dofs, u_lift = apply_dirichlet_by_lifting(
    K_global2, F_global2,
    constrained_dofs=[0, 1, 2],
    prescribed_values=[0.0, 0.0, 0.0]
)

print(f"Lifting: u = u_0 + u_D where")
print(f"  u_D (lift): {u_lift}")
print(f"  Free DOFs: {free_dofs}")
print(f"\nReduced system K_ff ({len(free_dofs)}x{len(free_dofs)}):")
print(np.round(K_ff, 3))
print(f"\nModified RHS F_f: {np.round(F_f, 4)}")

Lifting: u = u_0 + u_D where
  u_D (lift): [0. 0. 0. 0. 0. 0.]
  Free DOFs: [3, 4, 5]

Reduced system K_ff (3x3):
[[ 4. -1.  0.]
 [-1.  6. -1.]
 [ 0. -1.  2.]]

Modified RHS F_f: [0.25 0.5  0.25]

# Solve for the free DOFs
u0_free = np.linalg.solve(K_ff, F_f)

# Expand back to full solution
u0_full = expand_reduced_solution(u0_free, n_nodes, free_dofs, [0, 1, 2], [0.0, 0.0, 0.0])

# Total solution
u_lifting = u0_full + u_lift

print("Solution (lifting):")
for i in range(n_nodes):
    x, y = nodes[i]
    print(f"  u[{i}] = {u_lifting[i]:.6f}  (at ({x:.1f}, {y:.1f}))")

Solution (lifting):
  u[0] = 0.000000  (at (0.0, 0.0))
  u[1] = 0.000000  (at (0.5, 0.0))
  u[2] = 0.000000  (at (1.0, 0.0))
  u[3] = 0.095238  (at (0.0, 1.0))
  u[4] = 0.130952  (at (0.5, 1.0))
  u[5] = 0.190476  (at (1.0, 1.0))

10.6 Comparison of Both Approaches

print("Row substitution:", np.round(u_row, 6))
print("Lifting:         ", np.round(u_lifting, 6))
print()
print(f"Solutions match: {np.allclose(u_row, u_lifting)}")
print(f"Max difference: {np.max(np.abs(u_row - u_lifting)):.2e}")

Row substitution: [0.       0.       0.       0.095238 0.130952 0.190476]
Lifting:          [0.       0.       0.       0.095238 0.130952 0.190476]

Solutions match: True
Max difference: 0.00e+00

10.7 Exercise

Modify the problem setup above:

1. Change the Neumann flux value from 1.0 to 2.0 and re-solve using both approaches.
2. Apply a non-zero Dirichlet condition on the bottom: u(bottom) = 1.0 instead of 0.0
3. Change the boundary condition on the top edge to Dirichlet (u=0 on top, keep u=0 on bottom) and compare the solution.
4. Visualize the solution by plotting nodal values (e.g., a simple scatter plot at each node with the value as color).

For each case, verify that both methods still give the same result.

# Your exercise code here