QUAD4 bilinear shape functions on the reference square

Renders the four bilinear (4-node) shape functions \(N_i(\xi, \eta)\) on the reference quadrilateral \(\hat\Omega = [-1, +1]^2\):

\[N_i(\xi, \eta) = \tfrac{1}{4}(1 + \xi_i\xi)(1 + \eta_i\eta), \qquad i = 1, \dots, 4,\]

with corner-node signs \((\xi_i, \eta_i) \in \{-1, +1\}^2\). Each function equals 1 at its own node and 0 at every other node — the Kronecker-delta property at the heart of standard Lagrange interpolation. The four basis functions sum to 1 everywhere on the square (partition of unity), so any constant or linear field is reproduced exactly.

References

• Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002) Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, §6.3 (4-node quad), Table 6.3-1.

• Hughes, T. J. R. (2000) The Finite Element Method — Linear Static and Dynamic Finite Element Analysis, Dover, §3.5 + §3.6.

• Zienkiewicz, O. C., Taylor, R. L. (2013) The Finite Element Method: Its Basis and Fundamentals, 7th ed., §6.4.

Implementation: femorph_solver.elements.quad4_plane.Quad4Plane.

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np

Corner signs in VTK_QUAD order: (-1,-1), (+1,-1), (+1,+1), (-1,+1).

corner_signs = np.array(
    [[-1, -1], [+1, -1], [+1, +1], [-1, +1]],
    dtype=float,
)


def n_corner(xi, eta, signs):
    """Bilinear corner shape function for a given (xi_i, eta_i)."""
    xi_i, eta_i = signs
    return 0.25 * (1 + xi_i * xi) * (1 + eta_i * eta)

Sample the reference square on a (xi, eta) grid

n = 81
xi_grid, eta_grid = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n), indexing="ij")

Render — 2×2 grid of contour plots, one per corner.

fig, axes = plt.subplots(2, 2, figsize=(9, 7), constrained_layout=True)
fig.suptitle(
    "QUAD4 bilinear shape functions on $\\hat\\Omega = [-1, +1]^2$",
    fontsize=13,
)

for idx, signs in enumerate(corner_signs):
    row, col = divmod(idx, 2)
    ax = axes[row, col]
    N = n_corner(xi_grid, eta_grid, signs)
    img = ax.contourf(xi_grid, eta_grid, N, levels=21, cmap="plasma")
    ax.contour(
        xi_grid,
        eta_grid,
        N,
        levels=[0.0, 0.25, 0.5, 0.75, 1.0],
        colors="black",
        linewidths=0.4,
    )
    sx, sy = signs.astype(int)
    ax.scatter(corner_signs[:, 0], corner_signs[:, 1], s=30, color="black", zorder=5)
    ax.scatter([sx], [sy], s=110, edgecolor="black", facecolor="white", zorder=6)
    ax.set_aspect("equal")
    ax.set_title(f"$N_{{{idx + 1}}}$  $(\\xi_i, \\eta_i) = ({sx:+d}, {sy:+d})$", fontsize=11)
    ax.set_xlabel(r"$\xi$")
    ax.set_ylabel(r"$\eta$")
    fig.colorbar(img, ax=ax, shrink=0.85)

fig.show()

Verify partition of unity at every reference-element node.

total = sum(n_corner(corner_signs[:, 0], corner_signs[:, 1], s) for s in corner_signs)
np.testing.assert_allclose(total, 1.0, atol=1e-14)
print("OK — QUAD4 partition-of-unity holds at every corner node.")

OK — QUAD4 partition-of-unity holds at every corner node.

Verify Kronecker-delta interpolation at every reference-element node.

basis = np.zeros((4, 4))
for k, (xi_k, eta_k) in enumerate(corner_signs):
    for i, signs in enumerate(corner_signs):
        basis[k, i] = n_corner(xi_k, eta_k, signs)
np.testing.assert_allclose(basis, np.eye(4), atol=1e-14)
print("OK — QUAD4 Kronecker-delta interpolation verified at every node.")

OK — QUAD4 Kronecker-delta interpolation verified at every node.

Verify partition of unity on a dense interior grid as well.

dense_total = sum(n_corner(xi_grid, eta_grid, s) for s in corner_signs)
np.testing.assert_allclose(dense_total, 1.0, atol=1e-14)
print("OK — QUAD4 partition-of-unity holds across the full reference square.")

OK — QUAD4 partition-of-unity holds across the full reference square.