Using awatif-soil

1Scope

awatif-soil is a 2D plane-strain geotechnical finite element solver: six-noded triangles, drained Mohr–Coulomb with a tension cut-off, layered ground with a water table and a K₀ initial state, staged excavation with the water that goes with it, and retaining walls with Coulomb soil–wall interfaces and struts, written in C++ and driven from Python. Meshing is not part of it — the caller supplies the mesh, under the conventions of §4. This document is how to use it, ending in a complete reference of every interface (§11). The theory and the benchmark record — including the cross-validation against OpenGeoSys — are in paper.html, which ships alongside it.

The evaluation build is limited to 500 elements — building a bigger model raises ValueError. Everything else is the shipping solver, same API and numerics, and the limit is lifted by the unrestricted build; version() (§11) always tells you which build you have.

2Install

From the unzipped archive directory:

python -m pip install --upgrade awatif-soil==0.4.2 --find-links wheels

pip picks the wheel matching your platform and Python version out of wheels/ itself, and refuses with a clear message if none matches — no filename to choose by hand. Python 3.10 to 3.13; the only runtime dependency is numpy, which comes from PyPI as usual (on a machine with no internet access, install numpy first and add --no-index). For reference, what wheels/ contains:

cp310 to cp313 in the filename are Python 3.10 to 3.13. On Windows, use the standard 64-bit Python from python.org — including on ARM machines, where the x64 build runs under Windows’ built-in emulation; an ARM64 build of Python has no matching wheel here, and pip will say so.

Then run python selftest.py: it runs five numerical checks against closed-form answers and checks the evaluation cap at exactly 500 and 501 elements. It should end with All checks passed.

3First solve

import numpy as np
from awatif_soil import Model, material

clay = material(E=10e3, nu=0.3, c=10.0, phi=25.0, psi=5.0)

model = Model(nodes, elements, clay)      # nodes (n, 2), elements (m, 6)
model.fix(base_nodes, dx=0.0, dy=0.0)     # supports
model.fix(side_nodes, dx=0.0)
model.traction(top_edges, ty=-50.0)       # load

result = model.solve()
if not result.ok:
    raise RuntimeError(result.message)
print(result.displacements)               # (node, 2)
print(result.nodal_stress())              # (node, 4): xx, yy, zz, xy

examples/ has six runnable scripts: a footing driven to collapse (01_footing.py), layered ground under self-weight (02_layered_ground.py), two chained construction phases (03_staged_phases.py), layered ground under a water table started from the K₀ procedure and then dewatered (04_ground_conditions.py), a pit dug in two lifts below the water table, flooded and then pumped dry (05_staged_excavation.py), and a braced excavation — walls, a strut, a surcharge, four layers and a dewatered pit, eight phases (06_braced_excavation.py). A guided braced excavation runs in the browser at awatif.co/api/soil/demo, with construction stages, displacements, wall forces and a Python export. It is an illustrative example; a verified same-input Plaxis comparison is still outstanding.

4Conventions

These are the things that most often cost a day if they are assumed rather than read.

• Tension is positive, throughout the API and the internals, matching Plaxis. There is no conversion layer. Compressive stresses come back negative.
• Stresses are effective. All strength and stiffness parameters are effective parameters.
• Angles are in degrees at this boundary — phi and psi in material() — and converted to radians once, inside.
• Units are yours, used consistently. With kN and m, the Plaxis defaults, stresses come out in kPa.
• gamma and gamma_sat are unit weights: force per unit volume, acting downwards in −y, above and below the phreatic line. Gravity is already folded in, exactly as Plaxis states gamma_unsat and gamma_sat, so no 9.81 appears anywhere. A typical soil is 18 to 20.
• Water is a phreatic line, hydrostatic below it and no suction above, with the pore pressure negative in the water (tension positive) and gamma_w defaulting to Plaxis’s 10. The solver carries effective stress; the water is a load. A cluster of elements may carry a line of its own.
• Water standing on a boundary is applied automatically. On every exposed edge — the outer boundary and the faces an excavation opens — the water on the far side presses with its own pressure. A phreatic line above the ground surface is standing water, and a submerged support reports the effective reaction (§8).
• A phase declares everything that acts in it — supports, loads, self-weight, water, which elements exist, which walls and struts exist — and inherits stress from the previous phase. The solver ramps only what differs (§6, §8, §9).
• A wall’s sides are the left and right of its declared node order. Looking from the first node toward the last, left is the left hand: a wall declared top-down has its left on +x. Result.interface(wall, side) is named that way, and so is nothing else.
• Element nodes are in gmsh / Abaqus CPE6 order: corners 0, 1, 2 counter-clockwise, then mid-sides 3 = (0,1), 4 = (1,2), 5 = (2,0). Mid-side nodes are expected at the edge midpoints.
• Meshing is not part of this package. The solver takes the nodes and the connectivity you give it.

5Reading the result

solve() returns a Result whether or not the run converged — displacements, reactions, stresses, the yield map and the per-step history. Every member, with its exact shape, is specified in §11; what needs saying here is how to read a run that stops short.

A stalled solve is not an error and does not raise. Perfect plasticity has a genuine limit load, and past it there is no solution to find. A run that stops short still carries the achieved load_factor, the converged part of the step history, and the displacements that go with it — which is exactly what a load-settlement curve up to failure is made of. Check ok before treating the answer as complete, but do show the partial result: it is the useful part.

To reach a collapse load at all, prescribe displacement rather than force, and pass monitor= to solve() to get the reaction back per step. Force control has no solution past the limit load, by definition. A displacement-controlled collapse run:

clay  = material(E=10e3, nu=0.3, c=10, phi=25, psi=5)
model = Model(nodes, elements, clay)
model.fix(base_nodes, dx=0.0, dy=0.0)
model.fix(side_nodes, dx=0.0)
model.fix(footing_nodes, dy=-0.15)  # displacement, not force
result = model.solve(monitor=footing_nodes)

result.load_factor  # < 1 = incomplete; inspect result.message
result.steps        # load_factor, monitor_force per step: the curve

6Chaining construction phases

A phase declares everything that acts in it and inherits the previous phase’s stress through initial_stress:

phase2 = Model(nodes, elements, clay)
phase2.fix(base_nodes, dx=0.0, dy=0.0)        # everything that acts in this phase:
phase2.traction(top_edges, ty=-50.0)          #   the load phase 1 applied, still there
phase2.traction(new_edges, tx=-10.0)          #   and the one this phase adds
phase2.phreatic(-2.0)                         #   the water table (unchanged, or not)
phase2.initial_stress(phase1_result.state)    # where phase 1 left the soil

The solver reads what the inherited stress already carries straight off the stress and ramps only the difference between that and what the phase declares. So a load present in both phases is held, not re-applied; a load that this phase leaves out is unloaded, proportionally over the phase; a load that this phase adds is applied; and a water table that moved, or self-weight switched on, is ramped like any other load. This is how Plaxis defines a phase, and it is what lets self-weight and water be carried across phases at all.

The consequence to remember: a load left out is removed. Omit phase 1’s surcharge from phase 2 and it is quietly unloaded — the solve still converges and still reports Ok, and plasticity is path dependent, so the answer is simply wrong. examples/03_staged_phases.py shows the difference side by side. Declare, in every phase, everything that acts in it.

Displacements are reported per phase; Result.total_displacements adds them up across the phases chained through the state. Result.state also carries the load that was acting at the end of the phase, which is what keeps reactions exact while the next phase is under way; a plain list of Gauss-point states is accepted too, in which case reactions are exact at the end of a phase and approximate before, and the total restarts.

7Layers, water and initial stress

Layers are materials with material_ids. Water is a phreatic line per phase. The initial state comes either from gravity loading — a phase with gamma set and no inherited state, which switches self-weight on from zero stress and gives the elastic lateral ratio ν/(1−ν) — or from the K₀ procedure, which prescribes the lateral stress instead:

sand = material(E=30e3, nu=0.25, c=0.0, phi=34.0, psi=4.0,
                gamma=17.0, gamma_sat=20.0, k0=0.45)
clay = material(E=8e3, nu=0.35, c=15.0, phi=22.0, psi=0.0,
                gamma=16.0, gamma_sat=18.0)           # k0 defaults to 1 - sin phi

model = Model(nodes, elements, [sand, clay], material_ids)
model.fix(base_nodes, dx=0.0, dy=0.0)
model.fix(side_nodes, dx=0.0)
model.phreatic(-1.0)                          # water table 1 m down; gamma_w = 10

model.initial_stress(model.k0_state())        # K0 procedure: the initial state
result = model.solve()                        # nil-step: nothing moves

phase2 = Model(nodes, elements, [sand, clay], material_ids)
phase2.fix(base_nodes, dx=0.0, dy=0.0)
phase2.fix(side_nodes, dx=0.0)
phase2.phreatic(-4.0)                         # lower the water: ramped as a load
phase2.initial_stress(result.state)
dewatered = phase2.solve()

• phreatic(level) takes an elevation for a horizontal water table, or (x, y) points with increasing x for a general line — linear between them, flat beyond the ends. Below the line the pore pressure is p_w = −γ_wh and gamma_sat applies; above it p_w = 0 and gamma applies. Both are evaluated per Gauss point, so put water tables and layer boundaries on element boundaries: there the integration is exact, and an element the line cuts through is integrated approximately. pore_water() shows what the solver will use.
• k0_state() generates the initial stress from the profile: total vertical stress from the overburden, pore pressure off it, and σ′_h = k0 σ′_v with each material’s k0 (default 1 − sin φ′). Like Plaxis’s K₀ procedure it needs a horizontal ground surface and horizontal layers; anything else raises ValueError and wants gravity loading. Pass the state to a phase that declares the same ground and the same water: it starts in equilibrium and does not move.
• A k0 outside the active–passive range puts the state outside the yield surface. k0_state() reports how many Gauss points that is (violations, with a warning), and the phase that inherits the state corrects them as a plastic nil-step — the lateral stress lands on the active or passive value, the vertical stress stays.
• Lowering the water table is a phase that states the lower line and inherits the previous state. The effective stress rises where the water went, by γ_w − (γ_sat − γ) per metre dewatered above the point, and nothing else changes. examples/04_ground_conditions.py prints the whole profile against the closed form.

8Staged excavation

Excavation is element activation on top of the phase rule of §6. A phase declares which elements exist; an element left out has no stiffness, no weight, no water and no stress, and a phase that inherits stress from one in which it was active removes, ramped over the phase, the force that soil exerted on what remains. Nothing else is needed to dig:

lift1 = Model(nodes, elements, clay)
lift1.fix(base_nodes, dx=0.0, dy=0.0)
lift1.fix(side_nodes, dx=0.0)
lift1.phreatic(-1.0)                          # the retained water table
lift1.deactivate(pit_elements_above(-2.0))    # dig the first lift: those elements are gone
lift1.initial_stress(initial.state)           # the ramp removes the force they exerted
result1 = lift1.solve()

lift2 = Model(nodes, elements, clay)
lift2.fix(base_nodes, dx=0.0, dy=0.0)
lift2.fix(side_nodes, dx=0.0)
lift2.phreatic(-1.0)
lift2.deactivate(pit_elements_above(-4.0))    # both lifts: soil that is dug stays dug
lift2.phreatic(-4.0, elements=pit_column)     # pit and ground under it: pumped to the floor
lift2.initial_stress(result1.state)
result2 = lift2.solve()

result2.active                # (element,) bool
result2.total_displacements   # accumulated over every phase
lift2.surface_water()         # the water standing on the boundary, (node, 2)

• deactivate(elements) leaves the elements out of this phase; activate(elements) puts them back (every element starts active). Soil that is dug stays dug only if every later phase leaves it out — a phase declares everything. Soil put back — backfill — starts stress-free, whatever it carried before, and takes its weight through the ramp.
• Nodes no active element touches drop out: zero displacement, zero reaction, NaN in nodal_stress(). A point load or traction on such a node is refused (Status.InvalidMesh, naming the node) rather than dropped, because a surcharge left on excavated ground is a wrong declaration; a prescribed displacement there is ignored. Excavated elements read zero in gauss_stress().
• Water on the boundary. Every element, dug or not, carries a phreatic line — the global one unless phreatic(level, elements=...) gives a cluster its own — and the water standing on an exposed edge is applied automatically from the line on the far side of it. An excavated pit left on the global line is therefore flooded, its floor and lower face carrying the water column; a pit whose elements and the ground under them are given a line at the floor is pumped dry. surface_water() shows the load; pore_water() shows the pore pressure per element, dug elements included.
• The unbalanced water load. In a dewatered pit the pore pressure jumps across the pit’s outline: the retained water presses on the face between the table and the floor with nothing against it, and along the line under the pit’s edge the two pore pressures differ by γ_w times the drawdown. This is the load a retaining wall carries. Without a wall it has no equilibrium under a tension cut-off: the total stress on a free face is zero, so the effective stress there would have to be tensile by the pore pressure behind it, and the cut-off refuses it. The solve then stops short (NotConverged, with the load factor it reached, and the tension cut-off active along the face), and that is the solver reporting a real limit state. examples/05_staged_excavation.py shows a flooded pit converging and the same pit pumped dry stopping.
• total_displacements accumulates across the phases chained through the state, counted from when each element’s soil was placed: zero at removed nodes, and fill starts undeformed. Result.active says which elements took part.
• The K₀ procedure works on the active elements — the surface that must be horizontal is the exposed one — and refuses a cluster line: an initial state has one water table.

9Walls, interfaces and struts

A wall is declared on a line of mesh nodes the mesh conforms to — every consecutive (corner, mid, corner) triple an edge with soil on both sides — and the solver does the rest: the soil on each face gets its own copy of the nodes along the line, the wall keeps the originals, and a zero-thickness interface joins each pair of faces. Struts are bars between two nodes. Both are declared in every phase, active or not, exactly like elements, so the node numbering is the same throughout:

chain = nodes_on_line(x=3.5, top=18.5, toe=9.5)   # mesh nodes, top first, mid-sides included

initial = Model(nodes, elements, layers, material_ids)
initial.fix(base_nodes, dx=0.0, dy=0.0); initial.fix(side_nodes, dx=0.0)
initial.phreatic(17.26)
initial.wall(chain, EA=1.8e6, EI=2.3e4, R_inter=0.67, active=False)   # declared, not yet built
initial.strut(wall_node, axis_node, EA=2.7e6, spacing=3.0, active=False)
state = initial.k0_state()                    # ground alone: structures must be inactive

install = Model(nodes, elements, layers, material_ids)
install.fix(base_nodes, dx=0.0, dy=0.0); install.fix(side_nodes, dx=0.0)
install.phreatic(17.26)
install.wall(chain, EA=1.8e6, EI=2.3e4, R_inter=0.67)   # active: the tie is released, the
install.strut(wall_node, axis_node, EA=2.7e6, spacing=3.0, active=False)  #   earth pressure ramps onto it
install.initial_stress(state)
result = install.solve()

wall = result.wall(0)                         # .N .Q .M per metre at the nodes, .displacements, .rotations
face = result.interface(0, "right")           # declared top-down, right is -x: .normal .shear .mode
install.wall_water()                          # the water on the wall, (node, 2)

dig = Model(nodes, elements, layers, material_ids)
# ... supports, water, wall(active), then:
dig.strut(wall_node, axis_node, EA=2.7e6, spacing=3.0)   # installed: carries what moves from here
dig.deactivate(pit_elements_above(13.5))
dig.phreatic(13.5, elements=pit_column)       # pumped dry; the wall takes the unbalanced water
dig.initial_stress(result.state)
result = dig.solve()
result.strut(0)                               # kN per strut, tension positive

• The wall is Plaxis’s plate: a Mindlin beam with EA and EI per metre of wall, weight w per metre, and nu for its shear stiffness (κEA/(2(1+ν)), κ = 5/6). Each wall node gets a rotation (Result.rotations; fix_rotation holds it). Section forces come back per node through Result.wall(i): N, Q, M per metre, N tension positive, Q and M in the element’s frame (tangent along the declared order, normal to its left, M = EI dθ/ds), extrapolated from the two Gauss points.
• A wall toe inside the mesh is a free plate end in 0.4.2. With toe="free" (the default), the soil under the toe has one separate node shared by both sides, and the wall tip interacts with it through the last interface pair on each face. toe="tied" restores the previous wall-tip connection to the soil node. The free toe changes node_count and node_origin; rebuild every phase with the same wall declarations and toe choice rather than reusing old node maps or states.
• The interface is Coulomb’s, in effective stress. Elastic with kn and ks (stress per unit relative displacement); it slips at R_inter (c' − σ'n tan φ') of the adjacent soil and opens where its effective normal stress would exceed R_inter σ_t — the opening is a gap that must close before compression returns, and an open pair carries no shear. Integration is at the node pairs, so Result.interface(i, side) is nodal: normal, shear, gap, slip, mode along the wall. Left out, kn and ks follow Plaxis’s rule — a virtual thickness of 0.1 times the average element size, G_i = R_inter² G_soil, ν_i = 0.45 — from the soil on each side. That rule scales with R_inter squared: a wall made smooth by a tiny R_inter must be given kn and ks, or the interface is softer than the soil.
• Water reaches the wall. The soil face along an interface carries its own water like any exposed boundary (§8), so the interface sees effective stress and a retained face stays in compression; the wall receives the equal and opposite traction from each side (wall_water()). Two sides under different water levels leave the wall carrying the unbalanced water load; the dewatered pit of §8 that had no equilibrium has one with a wall. The interface transmits no pore pressure itself.
• An inactive wall is continuous soil. Its node pairs are tied, so the ground before installation is modelled on the same mesh; activating the wall in a later phase releases the tie and ramps the earth pressure onto the interfaces — the same rule that excavates soil. interface=False keeps the tie with the wall active: a wall bonded rigidly to the soil. Where the soil on one side is excavated, that side’s interface is off, and the water standing in the pit presses on the wall instead.
• Struts carry what moves after they exist. Force EA / (L spacing) times the change in distance between the two nodes, per metre of model; Result.strut(i) is the force per strut, tension positive — a strut jacked against the walls reads negative. prestress imposes an axial force during the phase (stiffness off, the change of force applied to the walls and ramped) and the strut leaves the phase holding it; later phases hold that force and add to it. A fixed-end anchor is a strut to a fixed node.
• Prescribing a wall’s displacement is fix(nodes, dx=, dy=, sides=False): the wall alone, the soil faces free against it. A plain fix on a wall node holds the wall and both faces, which is what a boundary support means.
• Node counts. With a wall declared, Model.node_count exceeds the mesh’s node count by the copies, and every per-node array — displacements, reactions, applied, the state’s — has that many rows; Model.nodes and Model.node_origin say where each row is and which mesh node it came from. The first mesh_node_count rows are the mesh’s own nodes, the wall nodes among them. Supports, loads and tractions are still given in mesh node ids; a traction on an edge that ends at a wall loads the side the edge belongs to.
• The initial state is ground alone. k0_state() refuses an active wall, interface or strut, as Plaxis does; declare them inactive there and activate them in a construction phase. examples/06_braced_excavation.py runs the whole sequence.

10Known limitations in this stage

Stated plainly, because they matter for how far you can take the integration:

• Edge tractions are uniform per edge — one constant vector for the whole edge. A surcharge that varies along the boundary can only be approximated as a staircase of edges, and one that follows a curved or inclined boundary cannot be represented. Water pressure does not depend on this: the water standing on any exposed boundary, curved or not, is applied automatically from the phreatic lines (§8).
• No arc-length control, so post-peak (softening) branches are not traced. Displacement control up to the limit load is supported and is what a curve needs.
• Collapse loads converge from above. Like every displacement-based FE code, a T6 mesh overestimates a limit load — about 5% when well-resolved and considerably more on a coarse mesh, so a coarse-mesh capacity is unconservative. Refine until the capacity settles before relying on it. The benchmark paper measures this convergence.
• Volumetric locking is marginal for fully incompressible plastic flow (phi' = 0); B-bar is not implemented.
• No seepage. A cluster line is a prescribed level, not a drawdown curve: between two levels the pore pressure jumps at the cluster boundary rather than falling through the ground. That is a drained analysis without a groundwater-flow calculation, as in Plaxis; an inclined global line is accepted where you want to prescribe one.
• A dewatered pit needs a wall on every face. Under a retained water table an unsupported face has no equilibrium with a tension cut-off (§8); expect NotConverged there. A half model with the pit centre as a symmetry axis needs one wall.
• Soil in tension is where the solver is weakest. Wherever the tension cut-off is active the exact tangent has no stiffness along the flow direction, and walls put soil in tension routinely. A load step the Newton iteration cannot converge is retried with the elastic soil matrix before the increment is bisected (fallback_iterations, fallback_tolerance), which carries the usual cases, R_inter = 1 and a wall pulled away from a c' = 0 sand included, at the cost of many more iterations. One does not always get through and stops with NotConverged: gravity loading from zero with an active wall glued to the soil by interface cohesion, which the initial phase without structures that Plaxis requires anyway avoids.
• Interfaces are on walls only, with Coulomb slip and a gap; no interface dilatancy or hardening, no interface between two soils.
• The K₀ procedure is for horizontal ground, like Plaxis’s; a sloping surface or non-horizontal layers take gravity loading. An inclined phreatic line is accepted but is not in equilibrium with a shear-free initial state, and the first phase corrects it.
• An element the water table cuts through is integrated approximately: put water tables and layer boundaries on element boundaries.

11API reference

Everything importable is listed here. The public surface is deliberately small: material, Model, Result, State, K0State, WallResult, InterfaceResult, Status, Region, InterfaceMode and version, plus a few module constants. State(states, applied=None, displacements=None, walls=None, interfaces=None, struts=None) can also be built by hand. Anything else reachable under awatif_soil._soil is internal and carries no compatibility promise. Arrays go in and come out as numpy; shapes are written (rows, columns).

material()

material(E, nu, c, phi, psi, sigma_t=0.0, gamma=0.0,
         gamma_sat=None, k0=None) -> MaterialMC

A drained Mohr–Coulomb material. All stiffness and strength parameters are effective.

Any violated bound raises ValueError. psi = phi is associated flow — accepted, but for soil it grossly overpredicts both dilation and collapse load. Without a cut-off (sigma_t=math.inf) the bare Mohr–Coulomb pyramid carries isotropic tension all the way to its apex at c′ cot φ′. Self-weight given through gamma is a load like any other: with no inherited state it is ramped from zero, so a solve with gamma set is gravity being switched on; with an inherited state that already carries it, it is held (§6).

On the returned MaterialMC object, phi and psi are stored in radians; build materials with material() rather than mutating the object’s attributes.

Model()

Model(nodes, elements, materials, material_ids=None)

nodes is (n, 2) float coordinates. elements is (m, 6) integer connectivity in the gmsh / Abaqus CPE6 node order of §4. materials is a single MaterialMC or a sequence of them; with more than one, material_ids is required — one index into materials per element, in element order (examples/02_layered_ground.py shows a two-material column).

A wrong shape, a missing material_ids, or a model over the evaluation element cap raises ValueError — the cap fires here, at construction, not after a solve you have already waited for. Model.node_count is the number of nodes, Model.element_count the number of elements. The module constants NODES_PER_ELEMENT = 6 and POINTS_PER_ELEMENT = 3 (Gauss points per element) name the two layout numbers the shapes below are built from; GAMMA_WATER = 10.0 is the default unit weight of water.

Supports and loads

model.fix(nodes, dx=None, dy=None, sides=True)
model.fix_rotation(nodes)
model.point_load(node, fx=0.0, fy=0.0)
model.traction(edges, tx=0.0, ty=0.0)

Each returns the model, so calls chain. Loads accumulate across calls. Together with self-weight and water they make up what the phase declares; the solver ramps whatever differs from the inherited state (§6) — from zero when there is none. Prescribed displacements are ramped from zero in every phase. A point load or a traction on a node no active element touches makes solve() return Status.InvalidMesh (§8); a prescribed displacement there is ignored.

sides matters only on a node a wall runs through: True (the default) holds the wall and both soil faces — a boundary — and False holds the wall alone, which is how a wall is given a prescribed displacement. fix_rotation holds the wall rotation at the given nodes at zero, Plaxis’s rotation fixity; ignored where no active wall reaches the node. A traction on an edge that ends at a wall loads the side the edge belongs to; one on the wall line itself raises ValueError.

• fix prescribes dx and/or dy on the given nodes; giving neither raises ValueError. 0.0 is a support; a non-zero value is displacement control (§5).
• point_load adds a nodal force (value at load factor 1).
• traction applies a uniform traction, force per unit length, to T6 edges given as (corner, mid-side, corner) node triples. The resulting nodal forces are the consistent Simpson (1, 4, 1)·L/6 ones, not an equal three-way split. One constant vector per edge — the limitation in §8.

Ground water

model.phreatic(level, gamma_w=10.0, elements=None)
model.pore_pressure(values)
model.pore_water() -> (ndarray (element, 3), ndarray (element, 3))
model.surface_water() -> ndarray (node, 2)
model.wall_water() -> ndarray (node, 2)

• phreatic sets this phase’s phreatic line: a number for a horizontal water table, or (k, 2) points with strictly increasing x for a general line, linear between the points and flat beyond the ends; None clears it. Unsorted points or a negative gamma_w raise ValueError. Below the line the pore pressure is −γ_wh and gamma_sat applies; above it the pore pressure is zero and gamma applies (§7). With elements the line belongs to that cluster alone (None makes it dry); a later call for some of the same elements overrides the earlier one, and inactive elements take part (§8). gamma_w is one value for the model.
• pore_pressure prescribes the pore pressure per Gauss point directly — 3 × n_elements values, flattened in element order, tension positive — and overrides the phreatic lines point by point; the choice of gamma against gamma_sat still follows the lines. A pore pressure on a free boundary is water standing on it.
• pore_water returns what the solver will use: the pore pressure and the unit weight at every Gauss point, each (element, 3), inactive elements included.
• surface_water returns the water standing on the boundary as nodal forces: on every exposed edge of the active mesh, the pore pressure on the far side times the outward normal, integrated along the edge. Applied automatically as part of the phase’s load; this shows it.

wall_water() is the water of the soil on each side of a wall pressing on the wall, as nodal forces on the wall nodes — minus each face’s share of surface_water(). Applied automatically; this shows it.

Staged excavation

model.deactivate(elements)
model.activate(elements)
model.active -> ndarray (element,) bool

deactivate leaves the given elements (indices into the connectivity) out of this phase and activate puts them back; every element starts active, and an index out of range raises ValueError. What that means for stress, nodes, loads and water is §8. Both return the model.

Walls, interfaces and struts

model.wall(nodes, EA, EI, w=0.0, nu=0.0, R_inter=1.0, kn=None, ks=None,
           active=True, interface=True, toe="free") -> int
model.strut(a, b, EA, spacing=1.0, prestress=None, active=True) -> int
model.nodes -> ndarray (node, 2)      model.node_origin -> ndarray (node,)
model.node_count, model.mesh_node_count, model.walls, model.struts -> int
model.elements -> ndarray (element, 6)

wall declares a wall along nodes, the chain of mesh nodes from one end to the other, mid-side nodes included (an odd count of at least three), and returns its index. Each consecutive (corner, mid, corner) triple must be an edge with soil on both sides, walls may not share nodes, and the parameters are checked (EA > 0, EI ≥ 0, w ≥ 0, 0 ≤ nu < 0.5, 0 < R_inter ≤ 1); anything else raises ValueError. What it builds and what the arguments mean is §9. The split happens as the wall is declared: node_count, nodes, node_origin and elements reflect it at once, and mesh_node_count stays the mesh’s. strut declares a bar between two mesh nodes with the axial stiffness of one strut and its out-of-plane spacing, returning its index; prestress is the axial force to impose in this phase, tension positive.

Result.wall(), Result.interface(), Result.strut()

result.wall(index=0) -> WallResult
result.interface(wall=0, side="left") -> InterfaceResult
result.strut(index=0) -> float        result.strut_forces -> ndarray (strut,)
result.rotations -> ndarray (node,)

Construction phases and initial stress

model.k0_state() -> K0State
model.initial_stress(state)
model.applied_force() -> ndarray (node, 2)

k0_state runs the K₀ procedure on this model’s ground profile and water (§7), over the active elements, and returns a K0State: .states (opaque), .applied (the self-weight and water it balances), .violations (Gauss points outside the yield surface; a warning is issued when it is not zero) and .stress(), (element, 3, 4). It raises ValueError for a non-horizontal exposed surface, non-horizontal layers, or a cluster phreatic line; water above the surface is standing water and is accepted. Point loads and tractions on the model play no part.

initial_stress takes the previous phase’s Result.state or a K0State — or a plain sequence of Gauss-point states, in which case the reactions of this phase are exact at its end and approximate before, and the displacement total restarts. Anything but 3 × n_elements states raises ValueError, and so does a state whose walls or struts do not match the ones declared — declare the same walls and struts, in the same order, in every phase. applied_force() returns this model’s own external load at load factor 1, (node, 2) — point loads, tractions, self-weight, water, wall weight and strut prestress assembled to nodes — and is informational: the solver works it out for itself. initial_force, from release 0.1.0, is gone; calling it raises TypeError naming the replacement.

Model.solve()

model.solve(monitor=None, **settings) -> Result

monitor is a set of node indices whose vertical reaction is summed and reported per step as StepRecord.monitor_force. This is how a load-displacement curve is obtained under displacement control: the total reaction over the whole boundary is zero by equilibrium, so a curve needs the reaction on a named subset — the loaded footing, say.

**settings overrides the solver settings below; an unknown name raises TypeError. The defaults are what every benchmark in the paper ran on, and exist to be read, not tuned.

Result

Returned by solve() whether or not the run converged (§5). Attributes first, then the methods that derive fields from the Gauss-point state.

StepRecord

One entry per attempted load step, all read-only. converged is False on a bisected attempt; failed attempts are kept on purpose — they are how the load path near a limit load is diagnosed — so a curve must filter on it rather than take every record. Note the name collision: StepRecord.plastic_points is a count, Result.plastic_points() is the boolean map.

Status, Region and InterfaceMode

Region names where the return map left each Gauss point’s stress, with principal stresses ordered σ₁ ≥ σ₂ ≥ σ₃, tension positive:

InterfaceMode is the regime of an interface pair in a phase:

version()

>>> import awatif_soil
>>> awatif_soil.version()
'awatif-soil a1b2c3d | eigen fetched 3.4.0 @ 3147391d | Release | evaluation build, 500 element limit'

Build provenance: the git commit the binary was built from, the exact Eigen it was linked against, the build type, and — on an evaluation build — the element limit (§13 asks for this line in every bug report).

12Licence and notices

LICENSE.txt governs your use of this software. THIRD_PARTY.md carries the notices for the two open-source libraries this binary is built with, Eigen and pybind11, both used unmodified.

13Support

Questions and bug reports: [email protected]. Please include the full version() line, your OS and Python version, and the output of selftest.py.