 codex/center-of-pressure-nonuniqueness.py

"""Original pressure profiles with identical resultant force and centroid."""
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def profiles(u):
    first = 1 + 0.6 * (2 * u - 1)
    return first, first + 0.3 * (6 * u**2 - 6 * u + 1)


def main():
    # Gauss quadrature integrates these low-degree polynomials exactly.
    nodes, weights = np.polynomial.legendre.leggauss(8)
    nodes, weights = (nodes + 1) / 2, weights / 2
    length, breadth, force = 0.4, 0.3, 800.0
    for pressure in profiles(nodes):
        total = np.dot(weights, pressure)
        first_moment = np.dot(weights, nodes * pressure)
        assert np.isclose(total, 1, rtol=0, atol=1e-14)
        assert np.isclose(first_moment, 0.6, rtol=0, atol=1e-14)
        assert np.isclose(length * first_moment / total, 0.24)
        assert np.isclose(breadth / 2, 0.15)
    p0 = force / (length * breadth)
    assert np.isclose(p0 * length * breadth, 800)
    u = np.linspace(0, 1, 1201)
    first, second = profiles(u)
    assert np.min(first) > 0 and np.min(second) > 0
    assert np.isclose(first.max(), 1.6) and np.isclose(second.max(), 1.9)

    fig, ax = plt.subplots(figsize=(8.6, 4.2), dpi=100, layout="constrained")
    fig.set_facecolor("white")
    ax.plot(u, first, color="#286f99", linewidth=2,
            label="Profile 1: linear pressure variation")
    ax.plot(u, second, color="#9f4b27", linewidth=2,
            label="Profile 2: an added quadratic variation")
    ax.axvline(0.6, color="#52762a", linestyle=":", linewidth=1.8)
    ax.annotate("Shared pressure-weighted position: x = 0.60L",
                xy=(0.6, 1.35), xytext=(0.42, 2.13), fontsize=10,
                arrowprops={"arrowstyle": "->", "color": "#52762a"})
    ax.text(0.025, 0.16, "Both: total normal force 800 N\nCenter of pressure (0.240, 0.150) m",
            fontsize=10, bbox={"facecolor": "white", "edgecolor": "#dddddd"})
    ax.set(xlim=(0, 1), ylim=(0, 2.4), xlabel="Position u = x / L",
           ylabel="Normal pressure / p₀",
           title="Same center of pressure; different local pressure loads")
    ax.legend(loc="upper left", fontsize=9)
    ax.grid(alpha=0.15)
    fig.savefig(Path(__file__).with_suffix(".png").name, facecolor="white")
    plt.close(fig)


if __name__ == "__main__":
    main()