"""An undrained clay slope on a foundation layer of different strength.

Example 4 of Griffiths & Lane (1999), "Slope stability analysis by finite elements",
Geotechnique 49(3), 387-403: a 2:1 slope of height H in clay of undrained strength
cu1, on a foundation layer of thickness H and strength cu2, with cu1 / (gamma H) = 0.25.
The factor of safety is computed by strength reduction for a range of cu2 / cu1.

Two classical limits bracket the answer (Taylor 1937):

    cu2  = cu1   deep-seated base failure through the foundation     FoS = 1.47
    cu2 >> cu1   toe failure confined to the slope                   FoS = 2.10

and the paper reports the change from one mechanism to the other at cu2/cu1 ~ 1.5.

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

H = 10.0                    # slope height = foundation thickness [m]
GAMMA = 20.0                # unit weight [kN/m3]
CU1 = 0.25 * GAMMA * H      # undrained strength of the slope clay [kPa]
RATIOS = (0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 3.0)


def slope(ratio, mesh_size=1.5):
    prj = katai.Project(f"Griffiths & Lane (1999) Example 4, cu2/cu1 = {ratio:g}",
                        mesh_size=mesh_size, auto_refine=False)
    # Undrained, total stress: phi = 0 and the cohesion is the undrained strength.
    upper = prj.materials.mohr_coulomb("Slope clay", E=1.0e5, nu=0.3,
                                       c=CU1, phi=0.0, gamma=GAMMA)
    lower = prj.materials.mohr_coulomb("Foundation clay", E=1.0e5, nu=0.3,
                                       c=ratio * CU1, phi=0.0, gamma=GAMMA)

    # Foundation: 6H wide, H thick. Fixed base, rollers on the two sides.
    prj.geometry.polygon([(0, 0), (6 * H, 0), (6 * H, H), (0, H)],
                         material=lower, name="Foundation",
                         fix=["full", "horizontal", "free", "horizontal"])
    # Slope: crest 2H wide, then a 2:1 face down to the toe at x = 4H.
    prj.geometry.polygon([(0, H), (4 * H, H), (2 * H, 2 * H), (0, 2 * H)],
                         material=upper, name="Slope",
                         fix=["free", "free", "free", "horizontal"])

    prj.initial(procedure="safety")     # strength reduction on the gravity state
    return prj


if __name__ == "__main__":
    print("cu2/cu1      FoS")
    for ratio in RATIOS:
        prj = slope(ratio)
        prj.save(f"clay-slope-cu-ratio-{ratio:g}.k2d")
        fos = prj.run().results()[-1].fos
        print(f"  {ratio:5.2f}   {fos:7.4f}", flush=True)