Stage 6 deformation analysis.

2.1 Plane strain

The deformation route is formulated in plane strain. Out-of-plane strain is zero, while the out-of-plane stress is recovered through the constitutive law.

This is appropriate for long strip-like loading, long embankment sections, or sections where out-of-plane variation is secondary. It is not a substitute for a true three-dimensional footing or excavation model.

2.2 Small strain

The kinematics are small-strain. Geometry is not updated during the solve, and the result is interpreted on the original mesh. This is appropriate for serviceability screening and early plastic-zone development, but not for large deformation or post-failure runout.

2.3 Effective-stress interpretation

The route is fundamentally a drained effective-stress screen. Initial pore pressure comes from either the seepage workspace or a hydrostatic phreatic-line reconstruction. The mechanical load step itself does not generate excess pore pressure.

2.4 Engineering interpretation of the current route

• Use it for long-term drained settlement and stress-screening questions.
• Use it for plastic-zone visualization and displacement trends.
• Do not read it as an undrained short-term clay model.
• Do not read Stage 2 plastic bands as a fully validated FEM failure surface in the PLAXIS sense.

3.1 Shared geometry

The deformation workspace reuses the Stage 6 section geometry rather than introducing a second modelling environment. The active terrain polyline, model base, lateral boundaries, soil polygons, retaining-wall geometry, seepage field, and measurement line all live in the same shared state.

• If custom polygons are disabled, the interpreted CPT layer column is extended laterally across the section.
• If custom polygons are enabled, those polygons define the deformation regions.
• The same mesh backbone is shared conceptually with seepage, even if the mechanical and hydraulic analyses are solved separately.

3.2 Material parameters carried into deformation

3.3 Load model

The user defines a terrain interval and applies either a direct pressure or an equivalent pressure derived from total load and out-of-plane length:

The current public route applies this as a vertical traction on the terrain interval. It does not yet introduce a rigid structural slab, contact law, interface element, or a structural plate stiffness in the deformation solve itself.

3.4 Support conditions and domain sizing

The standard screening supports are used on a truncated half-space:

This is acceptable only if the domain is sufficiently wide and deep. The technical guidance used in the spec is to keep the lateral and vertical extents on the order of several load widths, with roughly five load widths as a practical first-pass target.

3.5 Mesh model

The current route uses a constrained triangular mesh and supports two element families, dispatched per analysis through a unified element kernel:

• T3 — three-node constant-strain triangle, 1 Gauss point. Default for first-pass screening; fast and easy to audit.
• T6 — six-node quadratic triangle, 3 Gauss points. Used where bending response or near-incompressible behaviour matters; combines naturally with the B-bar near-incompressibility projection (§5.4).

Automatic meshing starts relatively coarse for first-pass user runs and can be manually refined where needed. The element type is selected through the analysis options; T6 increases the free-DOF count by roughly 4× over a comparable T3 mesh and the solver warns when the resulting DOF count is large enough to slow the run materially.

• Local load-edge refinement is important because stresses and strains concentrate there.
• Slope crests and geometry kinks can legitimately create high gradients and local plasticity.
• For T3, mesh convergence should be checked for sensitive bending-dominated problems.

4.1 Small-strain kinematics

The plane-strain engineering strain vector used in the element formulation is:

4.2 Full 3D stress or strain state inside the constitutive core

Although the global FE kinematics are plane strain, the constitutive logic is organized around full six-component stress or strain vectors, because Mohr-Coulomb yielding and plastic flow are naturally formulated in principal stress space.

This distinction is important: the deformation route uses engineering shear strain in Voigt notation. Therefore stress-like, strain-like, and gradient-like quantities are not interchangeable under a naive mapping, and work-conjugate handling of the shear terms must remain explicit.

4.3 Element strain-displacement and internal force

For a constant-strain triangle:

The global residual equation is:

4.4 Tangent stiffness

The consistent linearization of the global solve is expressed through the material tangent returned by the constitutive update:

For linear elasticity, D_tan = D^e. For Stage 1, the tangent is either elastic or reduced-shear. For Stage 2, the constitutive update returns the current elastoplastic tangent from the exact shear return; the default global solve uses its symmetrized form for robustness unless the optional unsymmetric path is enabled.

5.1 Isotropic elastic constants

These constants control both the linear elastic constitutive branch and the elastic part of the plastic constitutive update.

5.2 Full 3D elastic stress law

The full six-by-six isotropic elastic matrix is used inside the constitutive logic so that the out-of-plane normal stress is available for principal-stress evaluation and Mohr-Coulomb checking.

5.3 Plane-strain stress closure

Because ε_zz = 0 in plane strain, the elastic constitutive law generates an out-of-plane normal stress increment:

The incremental form matters. The absolute zero-state shortcut σ′_zz = ν(σ′_xx + σ′_yy) is not valid once an initial geostatic stress field exists and pore pressure has been subtracted.

5.4 Near-incompressibility handling

Plane-strain analysis can drive the volumetric response toward incompressibility, either through high Poisson ratio in elasticity or through low dilation angle (ψ → 0) under non-associated Mohr-Coulomb plastic flow. The current route treats this in two element-specific ways:

• T3 — single-Gauss-point constant-strain element. Volumetric locking is mitigated by capping Poisson’s ratio for numerical stability. The cap is a numerical safeguard, not a material law, and only matters as the ratio approaches 0.5.
• T6 — three-Gauss-point quadratic element. Volumetric locking is removed by the B-bar projection: the volumetric part of the strain-displacement matrix B is replaced by its element-average <B> while the deviatoric part is left untouched. The 2D plane-strain B-bar form is applied per Gauss point at element-stiffness assembly, so the same projection is consistent between K_tan and the internal-force integration.

6.1 Linear gravity predictor

The initial route always begins with a linear elastic gravity predictor on the same mesh and the same support set used by the subsequent deformation analysis. This predictor provides a geometry-aware vertical stress and initial shear field and forms the starting point for both shipped initial-stress modes.

This gives a geometry-driven total stress state including non-zero initial shear stress on slopes and near geometry breaks. The production Stage 2 route does not use that field blindly: it rebuilds the normal confinement from K_0,nc, keeps the recovered shear only as far as the exact Mohr-Coulomb and tension-cutoff envelope allow, and then solves the remaining self-weight equilibrium correction.

6.2 Effective stress reconstruction and K_0,nc predictor state

After the gravity step, pore pressure is subtracted from the normal total stress components only:

Storing the full six-component effective stress state is a crucial correction in the current solver architecture. Earlier shortcuts that stored only σ′_xx, σ′_yy, and τ_xy were not sufficient for a reliable plane-strain Mohr-Coulomb route.

The shipped predictor then reconstructs the in-situ normal confinement with the interpreted at-rest ratio rather than letting elastic Poisson coupling define the initial lateral stress state:

Here λ is chosen per integration point by admissibility clipping: λ = 1 when the full recovered shear is inside the exact Mohr-Coulomb and tension envelope; otherwise a bisection search selects the largest admissible value. This preserves the useful slope-equilibrium shear without seeding an inadmissible stress state.

6.3 Plastic self-weight equilibration

The predictor state above is not accepted as the final production initial condition. It becomes the starting point for a full Stage 2 self-weight equilibrium correction under exact Mohr-Coulomb elastoplasticity. Service loading and c-phi reduction require this self-weight phase to converge.

This phase is a predictor-correction problem, not a fresh loading from a stress-free state and not a replay of the gravity-step displacement field. The constitutive update sees the incremental strain operator B Δu around the seeded predictor stress state. That distinction is required to prevent double-counting of the gravity strain field.

• If the correction converges, the equilibrated stress and plastic state become the constitutive reference state for the service phase.
• If the correction does not converge, the service phase is not started and the reported state remains the last accepted self-weight state.

6.4 K_0,nc controls initial confinement, not ν

The present public route makes a strict conceptual distinction:

• K_0,nc controls the in-situ effective confinement state.
• ν controls elastic stiffness during the load solve.

This is important because allowing elastic ν to control initial confinement can make weak soils appear to “fail” everywhere under light load simply because the lateral stress state starts too low.

6.5 Displacement reset between the initial phase and service phase

If the plastic self-weight equilibration phase converges, the later service-load phase keeps the equilibrated stress state and plastic history but resets the displacement baseline for reported contour quantities. This matches the standard staged-analysis engineering reading used in commercial FEM workflows.

So the default deformation contours and settlement values still represent the service-load increment, not the sum of self-weight settlement and service settlement. The total and initial displacement fields are retained internally for diagnostics and reconstruction.

6.6 Flat K₀ fallback

If the gravity step fails numerically, the solver falls back to a flat-ground K₀ reconstruction. That fallback still constructs the total stress state first and only then subtracts pore pressure, so the out-of-plane effective stress is not corrupted by direct subtraction on an already reduced 2D state.

6.7 Effective stress and load increment

The FE solve returns an incremental stress field in the solver sign convention. For public engineering interpretation, the increment is converted back to geotechnical compression-positive convention and interpreted as an effective stress increment:

If a seepage result exists, its pore pressure field provides u₀. Otherwise the phreatic line provides a hydrostatic reconstruction.

7.1 Classical shear yield function

For compression-positive principal effective stresses σ′₁ ≥ σ′₂ ≥ σ′₃, the Mohr-Coulomb shear yield function is written as:

• f_s < 0: elastic admissibility.
• f_s = 0: on the yield surface.
• f_s > 0: elastic trial stress exceeds Mohr-Coulomb strength.

7.2 Principal stress evaluation

The constitutive algorithm works with the full compression-positive effective stress tensor in plane strain:

The full three-principal representation is required because the plane-strain constitutive state still carries an out-of-plane stress component and the exact Mohr-Coulomb active-set return is formulated in principal effective stress space.

7.3 Utilization ratio

This is a local reserve indicator, not a global factor of safety. Its interpretation is:

• η_MC < 1: inside the envelope.
• η_MC ≈ 1: near yield.
• η_MC > 1: trial stress exceeds Mohr-Coulomb strength.

7.4 Tension cut-off

Classical Mohr-Coulomb can imply an unrealistic effective tensile capacity. The exact Stage 2 route therefore includes the optional tension cut-off:

The exact constitutive implementation uses one physical tension surface only:

In Stage 1 this condition remains diagnostic-only. In Stage 2 it is now part of the exact active-set elastoplastic return. Tension-governed accepted states are reported with active yield surface TENSION; the app does not reinterpret those zones as ordinary finite shear-utilization states.

8.1 Additive strain split

Only Stage 2 computes true plastic strain. Stage 1 is not true plasticity and must not be described as such.

8.2 Non-associated plastic potential

Soils are typically modelled with non-associated flow:

The shipped Stage 2 route uses this exact principal-stress logic for the Mohr-Coulomb shear return, with face and edge handling in the local active-set solve.

8.3 Flow rule and Kuhn-Tucker conditions

These conditions define elastic loading, plastic loading, and unloading or reloading around a plastic state.

8.4 Backward-Euler return mapping

The Stage 2 constitutive path follows the standard elastic predictor/plastic corrector structure:

The shipped Stage 2 route uses an exact nonsmooth multisurface return in principal effective stress space. The accepted branch may be a shear face, a repeated-eigenvalue shear edge, the formal shear apex, a tension face, a mixed shear-tension branch, or the triple tension point.

8.5 Exact Stage 2 branch structure

The present implementation exposes the following exact Stage 2 branch families:

• MC face: single active shear surface.
• MC edge: two active shear surfaces with repeated principal stresses.
• MC apex: formal triple-shear branch where admissible.
• Tension face: single active cut-off surface T3.
• Mixed shear-tension: active set {F13, T3}.
• Tension apex: hydrostatic cut-off point with σ′₁ = σ′₂ = σ′₃ = −σ_t.

Repeated-eigenvalue branches are solved with representative principal-space normals and flow directions rather than by reusing the distinct-stress formulas unchanged. This is required for projector stability and exact branch closure near σ′₁ = σ′₂ or σ′₂ = σ′₃.

The mixed repeated shear-tension corners use branch-specific representative shear gradients in the repeated subspace; they are not treated as ad hoc stress clipping.

8.6 Equivalent plastic strain diagnostic

The public route reports an accumulated equivalent plastic strain for visualization and interpretation:

This is a diagnostic localization measure, not a direct serviceability quantity.

9.1 Linear elastic route

The linear elastic route is the baseline regression and comparison path. It uses the same mesh, the same geostatic preparation, and the same load model, but keeps the constitutive law purely elastic.

9.2 Stage 1 reduced-stiffness pseudo-plasticity

Stage 1 exists to make incremental loading mechanically meaningful before full plasticity. It is a pseudo-plastic route:

• Compute an elastic trial stress.
• Evaluate Mohr-Coulomb shear exceedance.
• If exceeded, recompute the increment with reduced shear stiffness.
• Do not return the stress exactly to the yield surface.
• Do not accumulate true plastic strain.

The current shipped Stage 1 branch remains monotonic or sticky within a run: once an element exceeds Mohr-Coulomb shear during that loading sequence, it remains on the reduced branch for the rest of that run. That is a documented implementation fact, not an idealized theory claim.

9.3 Stage 2 exact Mohr-Coulomb elastoplasticity

Stage 2 is the current public default and the exact elastoplastic constitutive route:

• It stores material-point plastic strain.
• It uses a local active-set return in principal effective stress space.
• It supports non-associated flow through ψ.
• It supports shear face, edge, apex, tension-face, mixed shear-tension, and triple tension-point branches.
• It returns the current elastoplastic tangent, with the default global solve using the symmetrized form for robustness unless the optional unsymmetric path is enabled.

The remaining limits are algorithmic and modelling limits, not missing core Mohr-Coulomb branch physics:

• it is not a large-deformation or post-failure formulation,
• it remains a drained small-strain section analysis on T3 / T6 triangles.

9A.1 Engineering positioning and intended use

The Hardening Soil (HS) model (Schanz, Vermeer, and Bonnier 1999) is the alternative drained constitutive route shipped alongside the linear and Mohr-Coulomb routes. It is positioned for problems where the constant Young modulus of Mohr-Coulomb is the dominant error source for displacement predictions: drained settlement of footings on real granular and lightly overconsolidated soils, excavation response in stiffness-stratified profiles, and any case where the primary-loading and unloading-reloading branches need to be distinguished.

• Use the linear elastic route for regression baselines and for elastic-comparison plots.
• Use Stage 2 Mohr-Coulomb elastoplasticity for limit-state screening, c-φ safety, and slope-stability-style analyses where the failure envelope dominates the engineering output.
• Use Hardening Soil when serviceability displacements matter and the stiffness-versus-stress behaviour of the soil is the controlling feature of the analysis.

9A.2 Stress-dependent power-law stiffness

HS replaces the constant E_mc of Mohr-Coulomb with three reference stiffnesses, each calibrated at a reference pressure p_ref (typically 100 kPa), and each scaled by a power law in confinement:

Typical exponent values: m ≈ 0.5 for cohesionless sands (Hardin and Drnevich 1972 scaling), m ≈ 1.0 for soft normally-consolidated clays, and m ≈ 0.7–0.9 for stiff clays and silts. E₅₀ and E_ur scale with the minor principal effective stress σ′₃ (the triaxial confining pressure), while E_oed scales with the major effective stress σ′₁ consistent with the one-dimensional oedometric loading direction.

9A.3 Hyperbolic primary-loading curve

HS reproduces the Duncan and Chang (1970) hyperbolic stress-strain relation in drained triaxial primary loading. For a deviator q below the Mohr-Coulomb failure deviator q_f, the axial strain follows:

The hyperbolic curve replaces the linear pre-yield response of Mohr-Coulomb with a continuously stiffening-then-softening secant stiffness as q approaches the Mohr-Coulomb envelope. The transition from elastic to plastic occurs at q = q_f as in Mohr-Coulomb, but the path to it is governed by the accumulated plastic shear strain γ^p through the cone hardening law of §9A.4.

9A.4 Two yield surfaces: shear cone and volumetric cap

HS uses two yield surfaces in principal effective stress space:

• A shear yield surface (the cone), parameterised by the plastic shear strain γ^p. The cone is Mohr-Coulomb-aligned: it depends on the major and minor principal effective stresses, mobilises with γ^p from a soft initial state, and saturates as the stress path approaches the Mohr-Coulomb failure envelope. Non-associated flow on the cone is governed by Rowe's (1962) stress-dilatancy through a mobilised dilation angle ψ_m.
• A volumetric cap yield surface, parameterised by the preconsolidation pressure p_p. The cap is an ellipse in the meridional plane (p′, q̃) and grows monotonically with cap-mode plastic volumetric strain. Associated flow on the cap means the plastic potential equals the yield function.

The two surfaces engage together along a K_0,nc normal-consolidation path: the cap controls the volumetric (oedometric) stiffness while the cone controls the deviatoric mobilisation, and the calibration in §9A.6 pins both the cap shape parameter M and the cap-hardening modulus H_cap from the user inputs E_oed^ref, K_0,nc, φ′, and ν_ur. Tension cutoff is enforced in the same hierarchy as the exact Mohr-Coulomb route: the tension surface dominates the HS surfaces and the return follows the existing project-and-check pattern.

9A.5 K_0,nc consolidation and cap-cone interaction

Under one-dimensional normally consolidated loading (ε_h = 0, σ_v ramped), both yield surfaces are active. The cap activation forces the volumetric stiffness to match E_oed^ref at p_ref; the cone activation, with mobilised dilatancy from Rowe, sets the lateral stress ratio σ_h/σ_v to the user-supplied K_0,nc. When K_0,nc is supplied as zero, the Jaky (1944) default K_0,nc = 1 − sin φ′ is used.

The unload-reload Young modulus E_ur is treated as the elastic stiffness inside both yield surfaces. The unload-reload Poisson ratio ν_ur is the elastic Poisson value (typically 0.2) used for the isotropic elastic tangent on which the return mapping is built.

9A.6 Parameter cookbook for typical soils

The parameter ranges below are reference values from Brinkgreve (2007), the standard PLAXIS HS calibration reference. They are starting points for preliminary analysis. Site-specific calibration against oedometer, triaxial, or field stiffness data should always supersede generic values.

Common supporting defaults: ν_ur ≈ 0.2, R_f = 0.9, p_ref = 100 kPa. K_0,nc defaults to 1 − sin φ′ when not entered explicitly. OCR = 1 places the initial cap exactly on the in-situ vertical effective stress; OCR > 1 shifts the cap outward and the cap remains inactive until subsequent loading reaches the preconsolidation level.

9A.7 Implementation status and current limits

The HS plugin is production-ready for uniform-load benchmark cases within the convergence envelope characterised in the implementation rollout. The current implementation has documented limits relative to the published reference formulation:

• WASM-only. The HS update is implemented in C++ and dispatched inside the WASM solver. The JavaScript reference path rejects HS analyses at run-start. The GPU pipeline does not support HS.
• Elastic-tangent globalisation. The shipped tangent assembly uses the elastic (E_ur-based) tangent for the global Newton solve rather than the consistent algorithmic tangent. Local admissibility per Gauss point is certified by the return mapping; the elastic globalisation costs a slower Newton convergence rate compared with a fully consistent tangent, but is robust and avoids the bookkeeping cost of differentiating the stress-dependent stiffness.
• Constant cap-hardening modulus. H_cap is calibrated once per region at material setup from the K_0,nc consistency condition and held constant during the analysis. PLAXIS uses a p_p- dependent H_cap that softens the cap as p_p grows; this refinement is tracked as a follow-up improvement (see §2.6 of the feature specification for the closed-form Brinkgreve relation).
• Narrow convergence envelope at extreme parameters. Combinations of low confinement, low cohesion, and high friction angle can drive σ′₃ close to zero, where the power-law denominator hits the numerical floor (0.5 kPa) and the cone hardening law becomes ill-conditioned. The tension cutoff usually engages before the floor matters, but for very soft surface layers the solver may need finer load-step settings.
• c-φ safety with HS is supported: c′, tan φ′, tan ψ, and the tension allowable are reduced by Σ_Msf at every active HS Gauss point; reference stiffnesses E₅₀^ref, E_oed^ref, E_ur^ref and the exponent m are not scaled, matching PLAXIS convention. The HS region constants M_cap, H_cap, and sin φ_cv are re-calibrated at every safety trial when φ′ changes via Σ_Msf.

10.1 Material-point state structure

The constitutive architecture uses committed and trial material-point state. This separates local constitutive updates from global equilibrium acceptance and supports:

• load stepping,
• commit or rollback,
• plastic strain accumulation only on accepted states,
• diagnostic comparison between displayed and committed states.

For plastic geostatic initialization the state container distinguishes predictorState, referenceState, committedState, and trialState. Predictor audit quantities remain attached to the seeded initial state, while equilibrated-initial audit quantities are stamped separately after a converged self-weight correction.

10.2 Global nonlinear solve

The global solve is residual-based:

Load stepping, line search, and cutback are used to keep the solve robust once plasticity activates.

10.3 Convergence and cutback

The public route exposes solver settings for nonlinear iteration limits, tolerances, initial load step, minimum load step, growth factors, cutback factors, and the optional unsymmetric plastic path. These are solver controls, not soil parameters.

• Stage 2 can use cautious growth after first yield.
• Plastic line-search failure or exhausted nonlinear iterations trigger cutback.
• Elastic low-load cases are intentionally kept close to the linear baseline.

10.4 Partial shown state near failure

If the nonlinear solve cannot fully converge but has already developed a physically meaningful near-failure state, the app now keeps and shows that best available state instead of discarding the entire run.

10.5 Three-phase staging: gravity → service → c-phi safety

The shipped solver runs as a three-phase staged analysis. Each phase reuses the previous phase's committed material-point state and plastic history; only the driving load and the convergence target change.

• initial-gravity — geostatic predictor and plastic self-weight equilibration (§6).
• service-load — drained service-load increment over the equilibrated initial state (§6.5, §6.7).
• safety-cphi — c-phi (φ′ / c′) strength reduction safety analysis. Strength is reduced by a per-phase factor; the solver follows the load-factor branch until equilibrium is no longer attainable. Reports a safety load factor with the meaning safety-strength-reduction, plus the safety-phase plastic increment and settlement.

10.6 Arc-length continuation in the safety phase

Prescribed-strength continuation (the default strength-control mode) parameterises the c-φ safety solve by the imposed Σ_Msf and Newton-iterates for displacement at each prescribed reduction level. The prescribed control variable becomes a poor path parameter near a limit point: the load step shrinks below minLoadStep · 4, line search stalls, the plastic active set oscillates near a fixed Σ_Msf, and the safety curve flattens before a coherent mechanism develops. At that point the solver can switch to a Crisfield–Riks spherical arc-length continuation in which displacement and continuation parameter are unknowns of the same Newton solve.

The corrector is the standard two-solve form (Crisfield 1981): at each Newton iteration the same tangent stiffness K is solved against −R and −∂R/∂λ, then the linearised constraint gives a scalar continuation correction and the displacement correction follows.

For c-φ safety ∂R/∂λ is computed by finite difference on Σ_Msf against the same committed material state — central difference when both probes lie in the admissible domain Σ_Msf ≥ 1, forward-only at Σ_Msf = 1 to respect the strength-reduction floor. Probe assemblies skip tangent construction (wantTangent=false) and never commit material state. The predictor sign is selected by dot-product against the previous accepted path increment, with positive Δλ forced on the first step within a phase; under the arcLengthAllowPostPeakSafetyPath option the predictor may descend (Σ_Msf decreasing with growing displacement) from the second step onward for diagnostic purposes.

The arc-length radius adapts after each accepted step from the achieved Newton iteration count relative to arcLengthTargetIterations, with separate shrink factors on rejection (controlled cutback) and on failure (faster cutback). Rejected steps restore committed displacement, material-point state, warm-start iterates, and active-set diagnostics exactly. Reported factor of safety remains the highest stable lower-bound Σ_Msf regardless of post-peak descent; post-peak curve points are recorded for diagnostic display only and never replace the peak value.

Implemented on the WASM CPU path only at present. The GPU v2 pipeline retains prescribed strength control; arc-length on GPU is deferred.

10.7 Material plugin registry

The deformation solver is constitutive-model-agnostic. Every constitutive route (linear elastic, Stage 1 reduced-stiffness, Stage 2 exact Mohr-Coulomb) is registered through the same material-plugin contract, dispatched by capability flags rather than by name string. Adding a new material model is a one-file change: write a factory returning a plugin object with the documented capability flags, then call registerMaterialPlugin(name, factory) once at module load. The solver picks up support for staged analysis, c-phi safety reduction, line-search, and predictor projection through the advertised flags.

Three plugins ship at present:

• linear-elastic — the baseline regression and comparison route.
• mc-reduced-stiffness — Stage 1 pseudo-plastic reduced-shear path.
• mc-plastic — Stage 2 exact Mohr-Coulomb elastoplasticity with tension cut-off.

10.8 GPU pipelines

Two WebGPU pipelines are shipped alongside the CPU solver, both opt-in and selected through the analysis options.

• v1 (runDeformationOnGpu) — linear-elastic only. Mirrors the relevant subset of the CPU analyzer's contract: K₀-controlled initial stress recovery, gravity, and surface load. Assembles a CSR-form K on device and runs CG with a 2×2 nodal block-Jacobi preconditioner. Bandwidth- bound at scale; remains the simpler audit path.
• v2 (runFullDeformationAnalysisOnGpuV2) — full elastoplastic resident pipeline. Single static-data handoff into device memory, then a resident Newton outer loop performs all three staged phases on device: geostatic predictor, service-load increment, and the optional c-phi safety reduction. The matrix-free K·x kernel recomputes element-wise B^TDB·u_e on every iteration through a precomputed incidence list, eliminating the global K assembly. Per-GP exact Mohr-Coulomb return mapping (face / edge / apex / tension / mixed branches) runs on device. CG is the primary linear solver; BiCGStab activates when the consistent tangent becomes unsymmetric under non-associated active plasticity, and as fallback. Jacobi or 2×2 block-Jacobi preconditioning per phase. Line-search residual trials happen on device; only one readback at the end.

Probing and fall-back: callers test WebGPU availability through probeGpuPipeline (v1) and the v2 controller's analogous probe; failure or unsupported scope automatically routes to the CPU path. The matrix-free Kx kernels and plastic Newton loop are parity-tested against the v1 and CPU references through scripts/verify_gpu_v2_parity.mjs.

11.1 Contour and legend system

The deformation workspace now supports high-contrast contour fills, isolines, and an in-canvas legend. The legend is placed inside the canvas and can be collapsed vertically.

11.2 Available displacement, strain, and stress fields

• u_x,fin, u_y,fin, |u|_fin, and settlement
• ε_xx,fin, ε_yy,fin, γ_xy,fin
• σ′_xx,init, σ′_xx,fin, σ_xx,init, σ_xx,fin
• σ′_yy,init, σ′_yy,fin, σ_yy,init, σ_yy,fin
• τ_xy
• η_MC and ε̄^p_acc

In the current effective-stress interpretation, pore pressure shifts only the normal stress components. Therefore τ_xy does not split into separate total and effective shear views.

When plastic geostatic equilibration is enabled and converges, the displayed displacement fields are the service-load increment relative to the equilibrated initial state. The solver still stores the total and initial displacement fields separately so the service increment can be reconstructed exactly as u_total − u₀.

11.3 Measurement line and line probe

The shared Stage 6 measurement line acts as a line probe in deformation. The graph and clipboard export follow the selected quantity and return distance-value data along the line.

Because the line probe samples the solved field, missing parts of the line outside the domain are shown as gaps rather than invented values.

11.4 Interpretation of η_MC and ε̄^p_acc

• η_MC is a pointwise reserve or utilization ratio, not a global factor of safety.
• In exact tension-cutoff-active zones, η_MC is suppressed because tension admissibility, not shear utilization, governs the accepted state.
• ε̄^p_acc is a plastic-history localization diagnostic, not a direct serviceability value.
• A plastic band is not yet automatically equivalent to a validated global failure plane.

12.1 Initial effective stress is stored as full stress₆

A major correctness fix in the current solver lineage was to stop storing only the in-plane effective stress components for the initial state. The public route now carries the full six-component effective stress vector through material-point seeding.

12.2 K_0,nc versus ν was separated explicitly

The current public route now uses K_0,nc for initial in-situ confinement and ν for elastic stiffness. This correction was necessary because using ν to control the initial lateral stress state can produce unrealistic blanket yielding in weak soils under light loading.

12.3 Stage 1 tension cut-off is diagnostic-only

Another important correction was to remove tension cut-off from Stage 1 reduced-shear activation. Stage 1 now interprets:

• MC-active zone as a reduced-shear pseudo-plastic branch,
• tension zone as a warning or diagnostic condition.

12.4 Stage 1 is still monotonic, not the ideal reversible active set

The correction plan for Stage 1 identified the mathematically cleaner final target as a reversible active-set problem. That is important theory, but it is not yet the shipped public behavior. The current public Stage 1 branch remains monotonic within a run for robustness.

12.5 Stage 2 is exact for the classical Mohr-Coulomb branch set

The current public default Stage 2 computes true plastic strain and performs an exact Mohr-Coulomb active-set return in principal effective stress space. The shipped branch set now includes the shear faces, repeated-eigenvalue shear edges, the formal shear apex where admissible, the tension face, the mixed shear-tension branches, and the triple tension point.

12.6 Tension cut-off is constitutive in Stage 2 and diagnostic in Stage 1

The public route now makes an explicit constitutive distinction:

• Stage 1: tension cut-off is diagnostic-only and never enters the reduced-stiffness branch.
• Stage 2: tension cut-off is part of the exact active-set return and is reported through the active surface label TENSION.

This prevents the earlier ambiguity where tension-governed accepted states could be displayed as if they were unresolved or infinite shear-utilization states.

12.7 Plastic geostatic initialization is a correction problem around the predictor

The initial plastic self-weight phase is formulated around the seeded geostatic predictor, not as a second gravity replay. The constitutive increment in that phase is driven by predictor-relative strain increments only. Predictor audit quantities and equilibrated-initial audit quantities are stored separately in the material-point state to avoid conflating inadmissible predictor diagnostics with the converged initial state.

14.1 Important present limits

• Section-based plane strain only.
• Small-strain kinematics only.
• T3 / T6 triangle formulations on the shared mesh; T6 with B-bar handles the near-incompressible case.
• Drained effective-stress loading step only; no coupled excess pore-pressure generation.
• No constitutive softening, fracture energy regularization, or strain localization control.
• No contact, interface, or structural-element coupling in the deformation solve.
• GPU v1 covers the linear-elastic route. GPU v2 covers the full elastoplastic pipeline including geostatic, service-load and c-phi safety phases with per-GP exact Mohr-Coulomb return mapping.

14.2 Natural next steps

• Consistent algorithmic tangent and p_p-dependent cap-hardening modulus for the Hardening Soil plugin (currently elastic-tangent globalisation with a constant H_cap); see §9A.7.
• HSsmall (small-strain stiffness) extension to the Hardening Soil family for cyclic and small-amplitude problems where the very-small-strain stiffness governs the response (Benz 2007).
• Arc-length continuation in the safety phase is available on the WASM CPU path (see §10.6) and is the default fallback when prescribed strength control approaches a limit point. The GPU v2 pipeline retains prescribed strength control only; arc-length on GPU is the next natural extension.
• Unsymmetric global Newton (GMRES with row/column equilibration) is the default linear solver whenever active non-associated plasticity is detected. The arc-length corrector reuses the same tangent across its two right-hand sides, with a scaling cache that skips redundant equilibration setup on the second solve.
• Hydro-mechanical and consolidation extensions beyond the present drained load-step assumption.
• T6 element coverage in the GPU v2 pipeline (currently T3-resident for the matrix-free Kx kernel).