r"""
L-shaped frame under tip load — Castigliano on a two-member portal
==================================================================

Two prismatic Bernoulli beams welded at a rigid corner form an
L-shaped portal frame.  A vertical column of length :math:`L_{v}`
runs from the clamp at the origin to the corner; a horizontal beam
of length :math:`L_{h}` cantilevers from the corner.  A
concentrated transverse load :math:`P` acts at the free tip
:math:`(L_{h}, L_{v})` in the :math:`-y` direction (downward).

Reading off the moments by superposition:

* **Column** (vertical, :math:`0 \le y \le L_{v}`): axial force
  :math:`N = -P` (compression); the tip load applied at horizontal
  offset :math:`L_{h}` produces a constant moment about :math:`z`,

  .. math::
     M_{\mathrm{col}}(y) \;=\; -\, P\, L_{h}.

* **Horizontal beam** (:math:`0 \le x \le L_{h}`): standard
  cantilever-with-tip-load moment,

  .. math::
     M_{\mathrm{beam}}(x) \;=\; -\, P\, (L_{h} - x).

Castigliano's theorem gives the tip deflection in the load
direction (Timoshenko & Young 1968 §80; Roark Table 9 case 6):

.. math::
    :label: l-frame-tip-deflection

    \delta_{\mathrm{tip}}
    \;=\;
    \frac{1}{E\, I}
    \Biggl[
        \int_{0}^{L_{v}} M_{\mathrm{col}}^{2}\, \mathrm{d}y
        +
        \int_{0}^{L_{h}} M_{\mathrm{beam}}^{2}\, \mathrm{d}x
    \Biggr]_{ \!\!\!/ P}
    \;=\;
    \frac{P\, L_{h}^{2}\, L_{v}}{E\, I}
    \;+\;
    \frac{P\, L_{h}^{3}}{3\, E\, I}
    \;+\;
    \frac{P\, L_{v}}{E\, A},

with the axial-strain term :math:`P L_{v} / (E A)` from the
column.  For a slender frame :math:`(I / A \ll L^{2})` the axial
contribution is negligible — at the default geometry below it
sits 6 000× below the bending term.

Two limits collapse the formula:

* :math:`L_{v} \to 0` — the column vanishes and the structure
  reduces to a horizontal cantilever:
  :math:`\delta = P L_{h}^{3} / (3 E I)`.  Recovers the
  :ref:`sphx_glr_gallery_verification_example_verify_cantilever_eb.py`
  result.
* :math:`L_{h} \to 0` — only axial compression of the column
  survives: :math:`\delta = P L_{v} / (E A)`.  Recovers the
  :ref:`sphx_glr_gallery_verification_example_verify_single_hex_uniaxial.py`
  / Hooke's-law axial response.

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

40-element BEAM2 (Hermite-cubic Bernoulli) line for the column,
40 more for the horizontal beam — sharing a common node at the
rigid corner so all moment / shear / axial transfer is implicit.
The clamp at :math:`(0, 0, 0)` pins all six DOFs; out-of-plane
DOFs are pinned across the line so the response stays strictly
2D in the :math:`x`-:math:`y` plane.

References
----------

* Timoshenko, S. P. and Young, D. H. (1968) *Elements of
  Strength of Materials*, 5th ed., Van Nostrand, §80
  (Castigliano on a curved-beam / frame).
* Roark, R. J. and Young, W. C. (1989) *Roark's Formulas for
  Stress and Strain*, 6th ed., McGraw-Hill, Table 9 case 6
  (right-angle bend, end load).
* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J.
  (2002) *Concepts and Applications of Finite Element
  Analysis*, 4th ed., Wiley, §16.5 (multi-element frame
  assembly).
"""

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver
from femorph_solver import ELEMENTS

# %%
# Problem data — equal-leg L-frame
# --------------------------------

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_v = 1.0  # column length [m]
L_h = 1.0  # beam length [m]
P = 1.0e3  # tip load [N, downward = -y]

EI = E * I_z

# Closed-form tip deflection (Castigliano, Roark Table 9 case 6).
delta_bending = P * L_h**2 * L_v / EI + P * L_h**3 / (3.0 * EI)
delta_axial = P * L_v / (E * A_section)
delta_tip_pub = -(delta_bending + delta_axial)  # downward in -y

print("L-shaped frame under tip load")
print(f"  column L_v = {L_v} m, beam L_h = {L_h} m, P = {P} N (-y)")
print(f"  E = {E:.2e} Pa, I = {I_z:.3e} m^4, A = {A_section:.3e} m^2")
print()
print("Closed-form references (Castigliano / Roark Table 9 case 6):")
print(f"  bending δ           = {-delta_bending * 1e3:+.4e} mm")
print("    = P L_h² L_v / E I + P L_h³ / (3 E I)")
print(f"  axial δ (column)    = {-delta_axial * 1e3:+.4e} mm   (= P L_v / E A)")
print(f"  total δ_tip         = {delta_tip_pub * 1e3:+.4e} mm")
print(
    f"  axial / bending ratio = {delta_axial / delta_bending:.2e}  "
    "(slender-frame regime — axial negligible)"
)

# %%
# Build the two-segment BEAM2 frame
# ---------------------------------

N_PER_SEG = 40

# Column: 0 → L_v along +y axis.
col_y = np.linspace(0.0, L_v, N_PER_SEG + 1)
col_pts = np.column_stack((np.zeros_like(col_y), col_y, np.zeros_like(col_y)))

# Beam: 0 → L_h along +x axis at y = L_v.  Skip the corner node (it's
# already in the column) and concatenate the rest.
beam_x = np.linspace(0.0, L_h, N_PER_SEG + 1)[1:]
beam_pts = np.column_stack((beam_x, np.full_like(beam_x, L_v), np.zeros_like(beam_x)))

pts = np.vstack((col_pts, beam_pts))

# Cells — column elements then beam elements.  Both segments share
# node ``N_PER_SEG`` (the corner), so beam-element indices have an
# implicit offset of ``N_PER_SEG`` for their first node.
cells = []
for i in range(N_PER_SEG):
    cells.append([2, i, i + 1])
# Beam segment: connect corner (node N_PER_SEG) to first beam point
# (which is the point right after corner = N_PER_SEG + 1).
cells.append([2, N_PER_SEG, N_PER_SEG + 1])
for i in range(1, N_PER_SEG):
    cells.append([2, N_PER_SEG + i, N_PER_SEG + i + 1])
cells_arr = np.array(cells, dtype=np.int64)

n_cells = cells_arr.shape[0]
grid = pv.UnstructuredGrid(
    cells_arr.ravel(),
    np.full(n_cells, 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),
)

n_nodes = pts.shape[0]
print(f"\nL-frame mesh: {n_nodes} nodes, {n_cells} BEAM2 cells")

# %%
# Boundary conditions
# -------------------
#
# Clamp at (0, 0, 0) — node 1.  Out-of-plane motion suppressed at
# every node so the response stays strictly 2D in x-y.

m.fix(nodes=1, dof="ALL")

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

# %%
# Tip load at (L_h, L_v, 0)
# -------------------------

i_tip = n_nodes - 1  # last node added (the beam tip)
m.apply_force(int(i_tip + 1), fy=-P)

# %%
# Solve + extract tip deflection
# ------------------------------

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_tip_fe = _node_dof(int(i_tip + 1), 1)
err = (abs(v_tip_fe) - abs(delta_tip_pub)) / abs(delta_tip_pub) * 100.0

print()
print(f"{'quantity':<22}  {'FE':>14}  {'published':>14}  {'err':>9}")
print(f"{'-' * 22:<22}  {'-' * 14:>14}  {'-' * 14:>14}  {'-' * 9:>9}")
print(
    f"{'v_tip (UY)':<22}  {v_tip_fe * 1e3:>10.4f} mm  "
    f"{delta_tip_pub * 1e3:>10.4f} mm  {err:>+8.4f}%"
)

assert abs(err) < 0.1, f"v_tip deviation {err:.4f} % exceeds 0.1 %"

# %%
# Render the deformed frame
# -------------------------

g = m.grid.copy()
disp_xyz = np.zeros((g.n_points, 3))
for ni in range(g.n_points):
    disp_xyz[ni, 0] = _node_dof(int(ni + 1), 0)
    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_tip_fe), 1e-12) * 0.1
warped = g.copy()
warped.points = np.asarray(g.points) + scale * disp_xyz

plotter = pv.Plotter(off_screen=True, window_size=(600, 600))
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)
plotter.add_points(
    np.array([[0.0, 0.0, 0.0]]),
    render_points_as_spheres=True,
    point_size=18,
    color="black",
    label="clamp",
)
plotter.add_points(
    np.array([[0.0, L_v, 0.0]]),
    render_points_as_spheres=True,
    point_size=14,
    color="#888888",
    label="rigid corner",
)
plotter.add_points(
    np.array([[L_h, L_v + scale * v_tip_fe, 0.0]]),
    render_points_as_spheres=True,
    point_size=14,
    color="#d62728",
    label=f"tip — v_tip = {v_tip_fe * 1e3:.4f} mm",
)
plotter.view_xy()
plotter.add_legend()
plotter.show()

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

print()
print("Take-aways:")
print(f"  • v_tip within {abs(err):.4f} % of P L_h² L_v / E I + P L_h³ / (3 E I) + P L_v / (E A).")
print(
    "  • Bending dominates the slender-frame response; the column's axial "
    f"compression contributes only {abs(delta_axial / delta_bending) * 100:.4f} % "
    "of the tip deflection at this geometry."
)
print(
    "  • Limit cases collapse cleanly: L_v → 0 ⇒ horizontal cantilever P L³/(3 E I); "
    "L_h → 0 ⇒ axial column P L / (E A)."
)