"""Rankine active and passive earth pressure, by finite elements.

A smooth vertical wall of height H retains a horizontal soil surface. The wall is
moved away from the soil (active) or into it (passive) and the force it carries is
read from the reactions. Because the wall and the base are smooth, the whole block
strains uniformly and reaches the Rankine state everywhere, so the closed forms are
exact for this boundary value problem:

    Pa = 1/2 Ka gamma H^2 - 2 c sqrt(Ka) H        Ka = tan^2(45 - phi/2)
    Pp = 1/2 Kp gamma H^2 + 2 c sqrt(Kp) H        Kp = tan^2(45 + phi/2)

and, when the soil cannot carry tension, the active force is that of the soil
below the tension crack of depth zc = 2 c / (gamma sqrt(Ka)):

    Pa = 1/2 Ka gamma (H - zc)^2

Units: kN, m.   Run:  python rankine_earth_pressure.py
Written for KATAI 2D 0.9.0 (https://kataigeo.com). Each case is also saved as a
.k2d project that opens in the desktop application or runs with `katai solve`.
"""
import math

import numpy as np

import katai

H, L = 6.0, 12.0          # wall height and length of the retained block [m]
GAMMA = 18.0              # unit weight [kN/m3]
PHI = 30.0                # friction angle [deg]
E, NU = 4.0e4, 0.3        # elastic constants [kN/m2, -]

KA = math.tan(math.radians(45.0 - PHI / 2.0)) ** 2
KP = math.tan(math.radians(45.0 + PHI / 2.0)) ** 2
K0 = 1.0 - math.sin(math.radians(PHI))


def wall_force(ux, c, tension_cutoff=True, save=None):
    """Move the wall by ux [m] (negative = away from the soil); return the force [kN/m]."""
    prj = katai.Project("Rankine wall", mesh_size=1.0, auto_refine=False)
    soil = prj.materials.mohr_coulomb("Soil", E=E, nu=NU, c=c, phi=PHI, psi=PHI,
                                      gamma=GAMMA, tension_cutoff=tension_cutoff)
    # Edges in order: base, far side, surface, wall.
    # "vertical" holds u_y only (a smooth base); "horizontal" holds u_x only (a roller).
    prj.geometry.polygon([(0, 0), (L, 0), (L, H), (0, H)], material=soil,
                         fix=["vertical", "horizontal", "free", "horizontal"])
    wall = prj.displacements.line((0, 0), (0, H), ux=ux, name="Wall movement")

    prj.initial(procedure="k0", exclude=[wall])          # at rest: K0 = 1 - sin(phi)
    prj.phases.plastic("Move the wall", activate=[wall])
    if save:
        prj.save(save)

    at_rest, moved = prj.run().results()
    on_wall = np.abs(np.asarray(moved.node_x)) < 1e-9
    force = lambda r: abs(np.asarray(r.reaction).reshape(-1, 2)[on_wall, 0].sum())
    return force(at_rest), force(moved)


def pct(a, b):
    return 100.0 * (a - b) / b


print(f"Ka = {KA:.6f}   Kp = {KP:.6f}   K0 = {K0:.6f}\n")

# --- 1. sand: c = 0 ----------------------------------------------------------
half = 0.5 * GAMMA * H * H
p0, pa = wall_force(-0.005, 0.0, save="rankine-sand-active.k2d")
_, pp = wall_force(+0.300, 0.0, save="rankine-sand-passive.k2d")
print("1. Sand (c = 0)")
print(f"   at rest   P0 = {p0:10.4f} kN/m   closed form {K0 * half:10.4f}   {pct(p0, K0 * half):+.4f} %")
print(f"   active    Pa = {pa:10.4f} kN/m   closed form {KA * half:10.4f}   {pct(pa, KA * half):+.4f} %")
print(f"   passive   Pp = {pp:10.4f} kN/m   closed form {KP * half:10.4f}   {pct(pp, KP * half):+.4f} %\n")

# --- 2. how much movement each state needs -----------------------------------
print("2. Wall force against wall movement (sand)")
print("     u [mm]    u/H [%]     P [kN/m]    K = P / (gamma H^2 / 2)")
for u in (-5.0, -2.0, -1.0, -0.5, 5.0, 20.0, 50.0, 100.0, 300.0):
    _, p = wall_force(u / 1000.0, 0.0)
    print(f"   {u:8.1f}   {100 * u / 1000 / H:8.3f}   {p:10.3f}   {p / half:8.4f}")
print()

# --- 3. a soil with cohesion: c = 10 kPa -------------------------------------
C = 10.0
zc = 2.0 * C / (GAMMA * math.sqrt(KA))
pa_full = KA * half - 2.0 * C * math.sqrt(KA) * H
pa_crack = 0.5 * KA * GAMMA * (H - zc) ** 2
pp_c = KP * half + 2.0 * C * math.sqrt(KP) * H
_, pa_cut = wall_force(-0.020, C, tension_cutoff=True, save="rankine-clayey-active.k2d")
_, pa_ten = wall_force(-0.020, C, tension_cutoff=False)
_, pp = wall_force(+0.300, C, save="rankine-clayey-passive.k2d")
print(f"3. Soil with cohesion (c = {C:g} kPa); tension crack depth zc = {zc:.4f} m")
print(f"   passive                 Pp = {pp:10.4f}   closed form {pp_c:10.4f}   {pct(pp, pp_c):+.4f} %")
print(f"   active, tension allowed Pa = {pa_ten:10.4f}   closed form {pa_full:10.4f}   {pct(pa_ten, pa_full):+.4f} %")
print(f"   active, no tension      Pa = {pa_cut:10.4f}   closed form {pa_crack:10.4f}   {pct(pa_cut, pa_crack):+.4f} %")