r"""
Simply-supported beam with off-center point load
================================================

A prismatic beam of length :math:`L` is simply supported at both
ends and loaded by a single concentrated transverse force
:math:`P` at :math:`x = a` (with :math:`b = L - a`).  The
Euler-Bernoulli closed form (Timoshenko 1955 §28; Roark Table 8
case 2) gives reactions

.. math::

    R_{\mathrm{left}} = \frac{P\, b}{L},
    \qquad
    R_{\mathrm{right}} = \frac{P\, a}{L},

a peak bending moment under the load

.. math::

    M_{\mathrm{load}} = \frac{P\, a\, b}{L},

and a deflection profile that splits at :math:`x = a`:

.. math::
    :label: ss-beam-off-center-deflection

    v(x) = \begin{cases}
      \dfrac{P\, b\, x\,(L^{2} - b^{2} - x^{2})}{6\, L\, E\, I},
        & 0 \le x \le a \\[8pt]
      \dfrac{P\, a\,(L - x)\,(2 L\, x - a^{2} - x^{2})}{6\, L\, E\, I},
        & a \le x \le L.
    \end{cases}

So the deflection at the load and the absolute maximum deflection
follow as

.. math::

    v(a)
    \;=\;
    \frac{P\, a^{2}\, b^{2}}{3\, E\, I\, L},
    \qquad
    v_{\max}
    \;=\;
    \frac{P\, b\, (L^{2} - b^{2})^{3/2}}{9\sqrt{3}\, E\, I\, L},

with the maximum occurring in the **longer** span at
:math:`x_{\max} = \sqrt{(L^{2} - b^{2}) / 3}` from the closer
support (Roark Table 8 case 2; both formulas valid for
:math:`a \le L/2`).

For the centred case :math:`a = L/2` the formula collapses to
:math:`v(L/2) = P L^{3} / (48\, E\, I)`, recovering the
:ref:`sphx_glr_gallery_verification_example_verify_ss_beam_central_load.py`
result.  Off-centring the load reduces both :math:`v_{\max}` and
:math:`M_{\max}` because the moment-arm distribution changes — at
:math:`a = L / 3`, the peak moment drops to
:math:`2 P L / 9` (vs :math:`P L / 4` for the centred case) and
the peak deflection falls to roughly 88 % of the centred value.

Implementation
--------------

A 60-element BEAM2 (Hermite-cubic Bernoulli) line spans the beam.
The two ends are simply supported (only :math:`u_y = 0`); the load
is applied at the closest mesh node to :math:`x = a`.  Hermite-cubic
shape functions interpolate Bernoulli kinematics exactly on the
piecewise-linear moment field this point load creates, so the FE
deflection at the loaded node and at the analytical peak position
match the closed form to machine precision.

References
----------

* Timoshenko, S. P. (1955) *Strength of Materials, Part I:
  Elementary Theory and Problems*, 3rd ed., Van Nostrand, §28
  (simply-supported beam, concentrated load at any point).
* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for
  Stress and Strain*, 6th ed., McGraw-Hill, Table 8 case 2.
* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.
  (2002) *Concepts and Applications of Finite Element
  Analysis*, 4th ed., Wiley, §2.5 (Hermite-cubic shape
  functions).
"""

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS

# %%
# Problem data — point load at L / 3
# ----------------------------------

E = 2.0e11
NU = 0.30
RHO = 7850.0
WIDTH = 0.05
HEIGHT = 0.05
A_section = WIDTH * HEIGHT
I_z = WIDTH * HEIGHT**3 / 12.0
I_y = HEIGHT * WIDTH**3 / 12.0
J = (1.0 / 3.0) * min(WIDTH, HEIGHT) ** 3 * max(WIDTH, HEIGHT)

L = 1.0  # span [m]
a_load = L / 3.0  # load position from the left support
b_load = L - a_load
P = 1.0e3  # transverse load [N, downward = -y]

EI = E * I_z

# Closed-form references (Roark Table 8 case 2).
R_left_pub = P * b_load / L
R_right_pub = P * a_load / L
M_load_pub = P * a_load * b_load / L
v_load_pub = -P * a_load**2 * b_load**2 / (3.0 * EI * L)

# Peak deflection: in the longer span (from a_load to L since a < L/2).
# x_peak measured FROM the right support; convert to absolute x.
x_peak_from_right = np.sqrt((L**2 - b_load**2) / 3.0)
x_peak_abs = L - x_peak_from_right
v_peak_pub = -P * b_load * (L**2 - b_load**2) ** 1.5 / (9.0 * np.sqrt(3.0) * EI * L)

print("Simply-supported beam, off-center point load")
print(f"  L = {L} m, a = L/3 = {a_load:.4f} m, b = 2L/3 = {b_load:.4f} m")
print(f"  P = {P} N (downward), E I = {EI:.3e} N m²")
print()
print("Closed-form references (Roark Table 8 case 2):")
print(f"  R_left  = {R_left_pub:+.4e} N      (= P b / L = 2P/3)")
print(f"  R_right = {R_right_pub:+.4e} N      (= P a / L = P/3)")
print(f"  M_load  = {M_load_pub:+.4e} N m    (= P a b / L = 2 P L / 9)")
print(f"  v(a)    = {v_load_pub * 1e3:+.4e} mm   (= P a² b² / (3 E I L))")
print(f"  v_max   = {v_peak_pub * 1e3:+.4e} mm   at x = {x_peak_abs:.4f} m")

# %%
# Build a 60-element BEAM2 line
# -----------------------------

N_ELEM = 60
xs = np.linspace(0.0, L, N_ELEM + 1)
pts = np.column_stack((xs, np.zeros_like(xs), np.zeros_like(xs)))

cells = np.empty((N_ELEM, 3), dtype=np.int64)
cells[:, 0] = 2
cells[:, 1] = np.arange(N_ELEM)
cells[:, 2] = np.arange(1, N_ELEM + 1)
grid = pv.UnstructuredGrid(
    cells.ravel(),
    np.full(N_ELEM, 3, dtype=np.int64),
    pts,
)

m = femorph_solver.Model.from_grid(grid)
m.assign(
    ELEMENTS.BEAM2,
    material={"EX": E, "PRXY": NU, "DENS": RHO},
    real=(A_section, I_z, I_y, J),
)

# %%
# Boundary conditions
# -------------------
#
# Simply supported at both ends — pin UY only.  Suppress out-of-
# plane motion + axial rigid body across the line; pin UX at the
# left support to remove the axial rigid-body mode.

m.fix(nodes=1, dof="UY")
m.fix(nodes=int(N_ELEM + 1), dof="UY")
m.fix(nodes=1, dof="UX")

for i in range(N_ELEM + 1):
    m.fix(nodes=int(i + 1), dof="UZ")
    m.fix(nodes=int(i + 1), dof="ROTX")
    m.fix(nodes=int(i + 1), dof="ROTY")

# %%
# Apply the off-center point load
# -------------------------------

i_load = int(np.argmin(np.abs(xs - a_load)))
m.apply_force(int(i_load + 1), fy=-P)

# %%
# Solve + extract
# ---------------

res = m.solve_static()
dof_map = m.dof_map()
disp = np.asarray(res.displacement, dtype=np.float64)


def _node_dof(node_id: int, dof_idx: int) -> float:
    """0=UX, 1=UY, 2=UZ, 3=ROTX, 4=ROTY, 5=ROTZ."""
    rows = np.where(dof_map[:, 0] == node_id)[0]
    for r in rows:
        if int(dof_map[r, 1]) == dof_idx:
            return float(disp[r])
    return 0.0


v_load_fe = _node_dof(int(i_load + 1), 1)


def v_closed_form(x: float) -> float:
    """Piecewise-cubic deflection profile (Roark Table 8 case 2)."""
    if x <= a_load:
        return -P * b_load * x * (L**2 - b_load**2 - x**2) / (6.0 * L * EI)
    return -P * a_load * (L - x) * (2.0 * L * x - a_load**2 - x**2) / (6.0 * L * EI)


# The Bernoulli-Hermite-cubic FE response is exact at every node
# for this problem.  The analytical absolute peak at
# ``x_peak_abs`` falls between mesh nodes — comparing FE vs that
# would conflate FE accuracy with mesh-sampling resolution.
# Cleaner check: at the FE peak node, the FE value should match
# the closed form *evaluated at the same x* to machine precision.

fe_uy_all = np.array([_node_dof(int(i + 1), 1) for i in range(N_ELEM + 1)])
i_peak = int(np.argmax(np.abs(fe_uy_all)))
v_peak_fe = float(fe_uy_all[i_peak])
v_at_peak_node_pub = v_closed_form(float(xs[i_peak]))

err_load = (abs(v_load_fe) - abs(v_load_pub)) / abs(v_load_pub) * 100.0
err_peak_node = (abs(v_peak_fe) - abs(v_at_peak_node_pub)) / abs(v_at_peak_node_pub) * 100.0

print()
print(f"{'quantity':<28}  {'FE':>14}  {'published':>14}  {'err':>9}")
print(f"{'-' * 28:<28}  {'-' * 14:>14}  {'-' * 14:>14}  {'-' * 9:>9}")
print(
    f"{'v(a) at load':<28}  {v_load_fe * 1e3:>10.4f} mm  "
    f"{v_load_pub * 1e3:>10.4f} mm  {err_load:>+8.4f}%"
)
print(
    f"{'v at FE peak node':<28}  {v_peak_fe * 1e3:>10.4f} mm  "
    f"{v_at_peak_node_pub * 1e3:>10.4f} mm  {err_peak_node:>+8.4f}%"
)
print(f"  FE peak-node x      = {xs[i_peak]:.4f} m  vs  analytical peak x = {x_peak_abs:.4f} m")
print(
    f"  analytical |v_max|  = {abs(v_peak_pub) * 1e3:.4f} mm  "
    f"(off-grid; FE-peak-node value above is on the same cubic curve)"
)

assert abs(err_load) < 0.1, f"v(a) deviation {err_load:.4f} % exceeds 0.1 %"
assert abs(err_peak_node) < 0.1, f"v at FE peak node deviation {err_peak_node:.4f} % exceeds 0.1 %"

# %%
# Plot the deflection profile against the closed form
# ---------------------------------------------------

fe_uy = fe_uy_all

import matplotlib.pyplot as plt  # noqa: E402

fig, ax = plt.subplots(1, 1, figsize=(7.5, 4.0))
xs_dense = np.linspace(0.0, L, 401)
ax.plot(
    xs_dense,
    [v_closed_form(x) * 1e3 for x in xs_dense],
    "-",
    color="#1f77b4",
    lw=2.0,
    label="closed form (Roark Table 8 case 2)",
)
ax.plot(xs, fe_uy * 1e3, "o", color="#d62728", markersize=4, label="BEAM2 FE")
ax.axvline(a_load, color="grey", ls="--", lw=1.0, label=f"load at x = a = {a_load:.4f} m")
ax.axvline(x_peak_abs, color="black", ls=":", lw=1.0, label=f"v_max at x = {x_peak_abs:.4f} m")
ax.axhline(0.0, color="black", lw=0.5)
ax.set_xlabel("x [m]")
ax.set_ylabel("v(x)  [mm]   (downward negative)")
ax.set_title("SS beam, off-center point load — FE vs closed form", fontsize=11)
ax.legend(loc="lower left", fontsize=9)
ax.grid(True, ls=":", alpha=0.5)
fig.tight_layout()
fig.show()

# %%
# Render the deformed beam
# ------------------------

g = m.grid.copy()
disp_xyz = np.zeros((g.n_points, 3))
for ni in range(g.n_points):
    disp_xyz[ni, 1] = _node_dof(int(ni + 1), 1)
g.point_data["displacement"] = disp_xyz
g.point_data["UY"] = disp_xyz[:, 1]

scale = 1.0 / max(abs(v_peak_fe), 1e-12) * 0.05
warped = g.copy()
warped.points = np.asarray(g.points) + scale * disp_xyz

plotter = pv.Plotter(off_screen=True, window_size=(900, 320))
plotter.add_mesh(g, color="grey", opacity=0.4, line_width=2, label="undeformed")
plotter.add_mesh(warped, scalars="UY", cmap="coolwarm", line_width=4)
support_pts = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]])
plotter.add_points(
    support_pts,
    render_points_as_spheres=True,
    point_size=18,
    color="black",
    label="simple supports",
)
plotter.add_points(
    np.array([[a_load, scale * v_load_fe, 0.0]]),
    render_points_as_spheres=True,
    point_size=14,
    color="#d62728",
    label=f"load (P = {P:.0f} N at x = {a_load:.3f} m)",
)
plotter.add_points(
    np.array([[xs[i_peak], scale * v_peak_fe, 0.0]]),
    render_points_as_spheres=True,
    point_size=12,
    color="#9467bd",
    label=f"v_max FE = {v_peak_fe * 1e3:.4f} mm",
)
plotter.view_xy()
plotter.add_legend()
plotter.show()

# %%
# Take-aways
# ----------

print()
print("Take-aways:")
print(f"  • v(a) at the load within {abs(err_load):.4f} % of P a² b² / (3 E I L).")
print(
    f"  • v at the FE peak node within {abs(err_peak_node):.4f} % of the closed-form "
    "deflection at the same x — Hermite-cubic FE is exact at every node for this load."
)
print(
    "  • Off-centering the load (a = L/3) reduces both peak moment "
    "(2PL/9 vs PL/4 centred) and peak deflection (~88 % of centred) — "
    "the moment-arm asymmetry dominates."
)