Geomechanics: Drucker-Prager strip footing

This example combines the two existing example-only geomechanics building blocks into one nonlinear boundary-value problem: the local Drucker-Prager triaxial constitutive driver and the elastic strip-footing setup. A rectangular soil block is loaded by a centered strip footing, the footing pressure is ramped in load steps, and the per-quadrature Drucker-Prager history is committed after each converged step.

It is deliberately example-only. It does not add a public geomechanics API; the constitutive code is copied locally (Torch only) so the script is self-contained. TensorMesh keeps the internal solid-mechanics convention stress tension-positive. For geomechanics reporting, settlement is shown positive downward,

\[s = -u_y.\]

This is a compact educational example, not a foundation-design method.

Problem

The model is a two-dimensional plane-strain-style soil block with unit out-of-plane thickness. The constitutive model is small-strain, associated Drucker-Prager plasticity with linear isotropic hardening, written internally with the tension-positive yield function

\[f(\sigma, \alpha) = q + \eta I_1 - (k + H\alpha) \le 0, \qquad I_1 = \mathrm{tr}(\sigma), \qquad q = \sqrt{\tfrac{3}{2}\,s:s}.\]

Boundary conditions

The example reuses the elastic footing roller setup:

• bottom boundary: vertical displacement fixed, u_y = 0;

• left and right boundaries: horizontal displacement fixed, u_x = 0;

• top boundary: free except for the loaded footing patch.

The footing pressure is lumped over the top-surface nodes inside the footing patch.

Solver

Because the problem is now path-dependent, it is solved by nonlinear energy minimization with load stepping. At each load step the total potential energy internal - external is minimized over the free displacement DOFs with L-BFGS. The example follows the same per-quadrature history lifecycle as the drucker_prager_triaxial example:

• per-quadrature history variables are stored in self.history[etype];

• previous-step eps_p and alpha are passed through element_data;

• update_state(u) is called after each converged load step under torch.no_grad().

Sanity checks

The script reports the final footing settlement, the maximum settlement, the maximum and mean committed plastic history, and the plastic centroid.

The associated test checks that the footing settles downward, that the committed plastic history and settlement grow monotonically with load, that the plastic zone localizes under the footing, that a low-load / high-cohesion case stays essentially elastic and matches the elastic footing settlement, and that a higher-cohesion soil develops less plasticity.

Fig. 47 Output of drucker_prager_footing.py. The left panel shows the deformed soil mesh colored by settlement -u_y. The middle panel shows the committed Drucker-Prager plastic history variable, which develops beneath the footing. The right panel shows the load-settlement curve, including nonlinear growth as plasticity accumulates. Deformations are exaggerated for visibility.

Running it

cd examples/solid/geomechanics/drucker_prager_footing
python drucker_prager_footing.py

For a fast numerical-only run without writing the plot:

python drucker_prager_footing.py --no-plot --steps 6 --chara-length 0.60

Core implementation

The load-stepped nonlinear solver is the heart of the example. It reuses the local Drucker-Prager assembler’s energy, element_data_from_history and update_state methods.

def solve_drucker_prager_footing(
    problem: FootingProblem,
    params: DruckerPragerParameters,
    n_steps: int = 10,
    dtype: torch.dtype = torch.float64,
    device: str | torch.device = "cpu",
    lbfgs_max_iter: int = 50,
) -> Dict[str, Any]:
    """Solve the footing problem with Drucker-Prager plasticity and load stepping.

    The total potential energy ``internal - external`` is minimized over the free
    displacement DOFs with L-BFGS at each load step.  After each converged step
    the per-quadrature Drucker-Prager history is committed with ``update_state``.
    """
    device = torch.device(device)
    mesh = build_mesh(problem, dtype=dtype, device=device)
    dp = DruckerPragerPlasticity.from_mesh(mesh, params=params)

    dim = mesh.dim
    n_points = mesh.n_points
    n_dof = n_points * dim

    fixed = make_boundary_mask(mesh, problem)
    rhs_full, loaded_nodes, total_vertical_load = make_load_vector(mesh, problem)
    rhs_full_flat = rhs_full.flatten()

    free_flat = (~fixed).flatten()
    free_idx = torch.where(free_flat)[0]
    base_zeros = torch.zeros(n_dof, dtype=dtype, device=device)

    free_u = torch.zeros(free_idx.numel(), dtype=dtype, device=device, requires_grad=True)

    def recover_full(free_values: torch.Tensor) -> torch.Tensor:
        """Scatter the free DOFs into a full [n_points, dim] displacement field."""
        u_flat = base_zeros.index_add(0, free_idx, free_values)
        return u_flat.reshape(n_points, dim)

    load_factors: List[float] = []
    pressures_kpa: List[float] = []
    footing_settlements: List[float] = []
    max_settlements: List[float] = []
    max_alphas: List[float] = []
    mean_alphas: List[float] = []

    for step in range(1, n_steps + 1):
        load_factor = step / n_steps
        f_ext_flat = load_factor * rhs_full_flat

        optimizer = torch.optim.LBFGS(
            [free_u],
            max_iter=lbfgs_max_iter,
            tolerance_grad=1.0e-10,
            tolerance_change=1.0e-14,
            history_size=50,
            line_search_fn="strong_wolfe",
        )

        def closure() -> torch.Tensor:
            optimizer.zero_grad()
            u_full = recover_full(free_u)
            internal = dp.energy(
                point_data={"displacement": u_full},
                element_data=dp.element_data_from_history(),
            )
            external = torch.dot(f_ext_flat, u_full.flatten())
            loss = internal - external
            loss.backward()
            return loss

        optimizer.step(closure)

        with torch.no_grad():
            u_full = recover_full(free_u)
            dp.update_state(u_full)

            settlement = -u_full[:, 1]
            load_factors.append(float(load_factor))
            pressures_kpa.append(float(load_factor * problem.footing_pressure / 1.0e3))
            footing_settlements.append(float(settlement[loaded_nodes].mean()))
            max_settlements.append(float(settlement.max()))
            max_alphas.append(float(dp.max_alpha()))
            mean_alphas.append(float(dp.mean_alpha()))

    with torch.no_grad():
        u_final = recover_full(free_u).detach()

    return {
        "mesh": mesh,
        "dp": dp,
        "u": u_final,
        "fixed": fixed,
        "loaded_nodes": loaded_nodes,
        "total_vertical_load_N_per_m": total_vertical_load,
        "load_factors": load_factors,
        "pressures_kpa": pressures_kpa,
        "footing_settlements_m": footing_settlements,
        "max_settlements_m": max_settlements,
        "max_alphas": max_alphas,
        "mean_alphas": mean_alphas,
        "n_nodes": int(n_points),
        "n_steps": int(n_steps),
    }

What’s next

This example stays intentionally example-only. A natural follow-up would be to promote a stabilized geomechanics assembler into tensormesh/assemble/ once the model surface and public API direction are agreed.