Geomechanics: elastic strip footing

This example adds a small geomechanics boundary-value problem to the solid mechanics example family. A linear-elastic soil block is loaded by a centered strip footing, then solved with TensorMesh’s standard direct linear-elasticity pipeline.

The example is deliberately modest: it demonstrates geotechnical boundary conditions, a footing load patch, settlement reporting, and reaction/load balance without adding a public geomechanics API.

Problem

The model is a two-dimensional plane-strain-style soil block with unit out-of-plane thickness. The displacement field satisfies small-strain linear elasticity,

\[\nabla \cdot \boldsymbol{\sigma} = \mathbf{0}, \qquad \boldsymbol{\sigma} = \mathbb{C} : \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} = \tfrac12(\nabla\mathbf{u} + \nabla\mathbf{u}^T).\]

The domain is loaded by a compression-positive footing pressure on a limited top-surface patch. Internally, the load is applied in the negative y direction. For geomechanics reporting, settlement is shown as positive downward,

\[s = -u_y.\]

Boundary conditions

The example uses a simple 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. This keeps the example close to the existing solid-mechanics direct solve examples, while the sanity check verifies that the total reaction balances the total applied load.

Sanity checks

The script reports the total applied vertical load, the vertical reaction at fixed vertical degrees of freedom, the reaction/load relative error, and the maximum settlement.

The associated test checks that:

  • the soil settles downward under the footing;

  • the vertical reaction balances the applied vertical load;

  • doubling the footing pressure approximately doubles the settlement.

Elastic strip-footing settlement contour and surface settlement profile

Fig. 46 Output of elastic_footing.py. The left panel shows the deformed soil mesh colored by settlement -u_y; red arrows mark the footing load, blue markers show vertical supports on the base, and purple markers show horizontal roller constraints on the side boundaries. The right panel shows the settlement profile along the ground surface. The largest settlement occurs beneath the centered footing, and the deformation is exaggerated for visibility.

Running it

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

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

python elastic_footing.py --no-plot --chara-length 0.5

Core implementation

The full driver is in examples/solid/geomechanics/elastic_footing/elastic_footing.py. It reuses the same direct-solve components as the existing linear-elastic solid examples: LinearElasticityElementAssembler, Condenser, and SparseMatrix.solve.

    problem: FootingProblem,
    dtype: torch.dtype = torch.float64,
    device: str | torch.device = "cpu",
) -> Dict[str, Any]:
    """Assemble and solve the elastic footing boundary-value problem."""
    mesh = build_mesh(problem, dtype=dtype, device=device)

    assembler = LinearElasticityElementAssembler.from_mesh(
        mesh,
        E=problem.E,
        nu=problem.nu,
    )
    K = assembler()

    fixed = make_boundary_mask(mesh, problem)
    rhs, loaded_nodes, total_vertical_load = make_load_vector(mesh, problem)
    rhs_flat = rhs.flatten()

    condenser = Condenser(fixed.flatten())
    K_cond, rhs_cond = condenser(K, rhs_flat)

    u_cond = K_cond.solve(rhs_cond)
    u_flat = condenser.recover(u_cond).to(dtype=rhs_flat.dtype)
    u = u_flat.reshape(mesh.n_points, mesh.dim)

    # Reactions at constrained DOFs: r = K u - f.
    # SparseMatrix inherits @ from torch-sla SparseTensor.
    residual = K @ u_flat - rhs_flat
    reactions = residual.reshape(mesh.n_points, mesh.dim)

    vertical_fixed = fixed[:, 1]
    vertical_reaction = reactions[vertical_fixed, 1].sum()

    uy = u[:, 1]
    settlement = -uy
    footing_settlement = settlement[loaded_nodes].mean()
    max_settlement = settlement.max()
    min_uy = uy.min()

    load_balance_error = torch.abs(
        vertical_reaction + torch.as_tensor(
            total_vertical_load,
            dtype=vertical_reaction.dtype,
            device=vertical_reaction.device,
        )
    ) / max(abs(total_vertical_load), 1.0)

    result = {
        "mesh": mesh,
        "K": K,
        "rhs": rhs,
        "u": u,
        "fixed": fixed,
        "loaded_nodes": loaded_nodes,
        "total_vertical_load_N_per_m": total_vertical_load,
        "vertical_reaction_N_per_m": float(vertical_reaction.detach().cpu()),
        "load_balance_relative_error": float(load_balance_error.detach().cpu()),
        "max_settlement_m": float(max_settlement.detach().cpu()),
        "footing_settlement_m": float(footing_settlement.detach().cpu()),
        "min_vertical_displacement_m": float(min_uy.detach().cpu()),
        "n_nodes": int(mesh.n_points),
        "n_loaded_nodes": int(loaded_nodes.numel()),
        "n_fixed_dofs": int(fixed.sum().detach().cpu()),
        "n_free_dofs": int(fixed.numel() - fixed.sum().detach().cpu()),
    }
    return result

What’s next

This example is elastic and intentionally small. A later follow-up could reuse the same footing geometry with the example-local Drucker-Prager material model, or promote a stabilized geomechanics assembler into tensormesh/assemble/ if the API direction becomes clear.