"""The vertical cut in undrained clay: how high can it stand?

A vertical face of height H in a purely cohesive (Tresca) soil of undrained strength
cu and unit weight gamma collapses when the stability number

    Ns = gamma H / cu

reaches a critical value. That value has no closed form, but limit analysis has
bracketed it between 3.7752 (lower bound) and 3.7776 (upper bound), so the factor
of safety of a cut is known to four figures:  FoS = 3.776 cu / (gamma H).

This script computes the factor of safety by finite element strength reduction on
several meshes, with the soil able to carry tension (the case the bounds are for)
and with the default tension cut-off (a soil that cracks).

Units: kN, m.   Run:  python vertical_cut_undrained_clay.py
Written for KATAI 2D 0.9.0 (https://kataigeo.com). The cases are saved as .k2d
projects that open in the desktop application or run with `katai solve`.
"""
import katai

H = 5.0           # height of the cut [m]
GAMMA = 20.0      # unit weight [kN/m3]
CU = 40.0         # undrained shear strength [kPa]
NS_EXACT = 3.776  # between the bounds 3.7752 and 3.7776


def cut(element, mesh_size, tension_cutoff):
    prj = katai.Project("Vertical cut in undrained clay", element=element,
                        mesh_size=mesh_size, auto_refine=False)
    clay = prj.materials.mohr_coulomb("Clay", E=1.0e5, nu=0.3, c=CU, phi=0.0,
                                      gamma=GAMMA, tension_cutoff=tension_cutoff)
    # L-shaped section, counter-clockwise from the bottom-left corner: retained
    # ground 3H wide and 2H deep, the cut face at x = 3H, the floor 2H wide.
    prj.geometry.polygon(
        [(0, 0), (5 * H, 0), (5 * H, H), (3 * H, H), (3 * H, 2 * H), (0, 2 * H)],
        material=clay,
        fix=["full", "horizontal", "free", "free", "free", "horizontal"])
    prj.initial(procedure="safety")
    return prj


if __name__ == "__main__":
    target = NS_EXACT * CU / (GAMMA * H)
    print(f"limit analysis: Ns = {NS_EXACT}, so FoS = {target:.4f}\n")
    print("element  mesh [m]  tension     FoS      Ns = FoS gamma H / cu   vs 3.776")
    for element, mesh, cutoff in [("tri6", 2.0, False), ("tri6", 1.0, False),
                                  ("tri6", 0.5, False), ("tri15", 2.0, False),
                                  ("tri15", 1.0, False), ("tri15", 2.0, True)]:
        prj = cut(element, mesh, cutoff)
        prj.save(f"vertical-cut-{element}-{mesh:g}m-{'cracking' if cutoff else 'tension'}.k2d")
        fos = prj.run().results()[-1].fos
        ns = fos * GAMMA * H / CU
        print(f"{element:6s}   {mesh:5.2f}    {'cut off' if cutoff else 'carried'}   {fos:7.4f}   {ns:8.4f}"
              f"            {100 * (ns - NS_EXACT) / NS_EXACT:+6.2f} %", flush=True)