Skip to content

Implicit contact on prescribed DOFs and sphere-indentation screen

Status: experimental source candidate, not a Hertz validation or a general contact promotion. The source branch follows merged PR #116.

Numerical correction

The serial finite-J2 rigid-pair route previously rejected contact trace DOFs shared with strong Dirichlet constraints. It now lifts the native contact tangent by the same prescribed correction g_bc - u_trial used by DOLFINx for the UFL material tangent. Native COO entries on constrained rows and columns are omitted after the material matrix has had its strong rows/columns set by DOLFINx. This preserves the total-operator action on free rows and the identity rows on constrained DOFs. It does not change the frictionless penalty law, rigid projection or material algorithm.

One active, tilted triangular-tool/J2 test imposes a nonzero displacement on the contact face, completes, and checks every shared strong row remains an identity row. Another fixed-old-state directional test compares the assembled total J2-plus-contact Newton tangent against a central difference of the complete correction residual, with the same trial material state and zero perturbations on constrained DOFs. The shared strong-DOF lifting regression also passes with CG2 displacement on a Q1 hexahedral geometry, checking active contact, unchanged identity rows and the free-row contact correction. This is a second supported discretization check, not evidence for arbitrary geometry orders or incompatible spaces.

A separate sphere-indentation test compares a full block with the symmetry-reduced quarter block at identical physical resolution: 256 versus 64 hexahedra give the same full-body force to better than 1e-9 relative difference, and four times the quarter contact potential equals the full-block potential. These are structural equivalence checks, not independent physical validation. Existing contact/MPC regression remains required.

Independent reference screen

examples/finite_j2_hertz_contact.py solves a prescribed rigid-sphere indentation of a finite, bottom-clamped block. Its very high yield threshold keeps the finite-strain J2 material in its elastic regime. The report keeps the model separate from the classical small-strain, elastic-half-space, exact-nonpenetration Hertz limit F = 4 E sqrt(R) delta^(3/2) / (3 (1-nu^2)). The Abaqus analytical Hertz relation states that reference relation and its domain; it does not certify the finite block or penalty contact used here.

Commands use the current source checkout and the reviewed FEniCSx runtime:

PYTHONPATH=src OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 python examples/finite_j2_hertz_contact.py --quarter --cells 8
PYTHONPATH=src OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 python examples/finite_j2_hertz_contact.py --quarter --cells 8 --increments 16

At E=1000, nu=0.3, R=2, indentation 0.02, penalty 1e5, four increments and unit half-width/depth, the half-space reference force is 5.86080586. Full-block mesh levels with 256, 2048, 6912 and 16384 hexahedra produce forces 7.99194, 7.32959, 7.16618, and 7.09758. The finest screened ratio is still 1.2110, so this is not a passed Hertz accuracy gate. Finite boundaries, contact quadrature, mesh resolution and finite penalty remain confounded. A deeper quarter block (half-width and depth 2, 4096 hexahedra) gives 6.65969 but changes both boundary distance and resolution, so it cannot independently identify the discrepancy.

At the same extent-3, 13,824-cell quarter-block resolution, separable grading 2.5 reduces the top-center reference edge spacing from 0.125 to 0.02945 without changing the cell count. The active quadrature count increases from 8 to 102, the maximum sampled active radius is 0.20194 versus the Hertz radius 0.2, and the full-body force changes from 6.45723 to 6.20730 (ratio to the half-space limit 1.1018 to 1.0591). This is evidence that under-resolving the contact patch materially biases the earlier comparison. It is not a pass: the mesh map changes the discrete problem, the result still differs by 5.9%, and penalty/domain limits remain unseparated. The peak value 69.737 is nominal force per reference area at a sampled quadrature point, not a smoothed nodal pressure. The Hertz peak for these inputs is about 69.96, so the sampled center peak is close while the integrated force is still high; pressure shape and reference-area activity therefore matter. The active reference area 0.13521 exceeds the Hertz circle area 0.12566; neither a matching peak nor a near-matching sampled radius establishes pressure-profile agreement. At fixed extent 3, grading 2.5 and penalty 1e5, a 4,096-cell quarter block gives full force 6.25502, versus 6.20730 for 13,824 cells; active quadrature points rise from 43 to 102. At the 4,096-cell level, raising the penalty to 2e5 changes full force to 6.34303 and maximum penetration from 0.0007015 to 0.0003581. Thus the penalty limit does not explain away the positive force mismatch, and the spatial sequence is not yet flat.

At 512 quarter-block cells, penalties 5e4, 1e5, 2e5 give full-body forces 7.10580, 7.32959, 7.38030 and maximum penetrations 0.0015424, 0.0010021, 0.00069635. At fixed mesh/penalty, 4, 8 and 16 increments leave the final force unchanged within the nonlinear tolerance. Their accepted trapezoidal work residual declines 1.8809e-4, 5.8245e-5, 2.0718e-5; this is a path-discretization diagnostic, not a general first law for plastic contact.

Linear-solver screen, not a new global default

The example now accepts an explicit --linear-policy so the same model can compare the existing direct LU route with PETSc GMRES/GAMG or GMRES/ILU. On this Apple Silicon FEniCSx environment, the 4,096-cell, extent-3 graded quarter model takes 32.84 s with direct LU and 11.65 s with GMRES/GAMG; the full-body force differs by about 1.3e-11 and both complete 17 Newton linear solves. Accepted linear-solve phases are 22.43 s versus 1.52 s; material update and line search are not materially accelerated. The 13,824- cell GMRES/GAMG case completes in 48.81 s with full force 6.20730, matching the earlier direct-route value at reported precision. All its 21 linear solves have positive convergence reasons. The 32,768-cell GMRES/GAMG case completes in 113.34 s with full force 6.18308, active reference area 0.13663 and peak nominal pressure 69.65; all 20 linear solves converge. Its force ratio to Hertz is still 1.0550, so the independent accuracy gate is not passed by mesh refinement alone. A 32,768-cell direct-LU screen was stopped after 20 minutes while PETSc was in numeric LU factorization; it produced no scientific result and is not a measured direct-versus-iterative speed ratio. These timings are case/hardware evidence, not a portable choice for every contact model.

One far-boundary screen increases the quarter-block half-width/depth from 3 to 4, retaining 13,824 cells and adjusting the separable grading from 2.5 to 3.0 so the top-center spacing remains of similar order (0.02945 versus 0.02791). With GMRES/GAMG it completes in 47.60 s, gives force 6.12050 (Hertz ratio 1.04431), active area 0.13266, and 102 active quadrature points. The decline from 6.20730 is consistent with finite-domain influence, but spacing and element distortion are not exactly identical; this is not an isolated domain-size convergence or a certified Hertz result. At extent 5, 13,824 cells and grading 3.35, the top-center spacing is 0.02723; the force is 6.07176 (ratio 1.03599), active area 0.13113, and 100 quadrature points are active. This strengthens the diagnosis that the earlier small finite block elevated the force, while still leaving a 3.6% gap and confounded mesh-map/penalty/finite-strain effects.

An additional fixed extent-8 and grading-4 spatial sequence with penalty 1e6 and four increments gives the following source-branch diagnostics:

Quarter cells Top-center spacing Active points Full force Hertz ratio Wall time
4,096 0.04239 32 6.28701 1.07272 15.98 s
13,824 0.02707 86 6.16946 1.05266 56.07 s
32,768 0.01987 156 6.11118 1.04272 143.23 s

The trend is monotone under this mesh map, and the last sampled active radius is 0.19840 versus the Hertz 0.2. However, the force still differs by 4.27%; the active reference area is 0.12044 versus the Hertz circle's 0.12566. This sequence does not separate penalty, far-boundary and finite-strain errors. It is therefore diagnostic evidence, not an independent pass or a solver speedup claim. These runs used the source branch with one process and PETSc GMRES/GAMG, not an installed release.

An opt-in Q1 hexahedral trace integration screen now permits degree 5 (3×3 points) or 7 (4×4 points) while keeping the established 2×2 rule as the default. It changes only contact integration, not the bulk J2 rule or the model geometry, and its quadrature identity enters restart checks. For the same extent-8/grade-4/penalty-1e6 4,096-cell quarter model, 3×3 gives full force 6.29810 (Hertz ratio 1.07461) versus default 6.28701 (1.07272); wall time is 17.89 s versus the previous 15.98 s source run. At 13,824 cells, 3×3 gives 6.18311 (ratio 1.05499) versus default 6.16946 (1.05266); wall time is 67.47 s versus 56.07 s. These are not paired performance trials, but they show no accuracy or speed benefit for this benchmark. More contact points alone do not resolve the remaining Hertz discrepancy. Keep the higher rule as an explicit quadrature-sensitivity diagnostic; do not promote it to a default or use it to claim validation. For both permitted higher rules, separate one-cell sphere and cylinder tests match an independently assembled UFL contact potential, residual and Hessian on the same CG1 or CG2 hexahedral trace. A two-rank trace test checks partition-independent point IDs and total reference weight for the 2×2, 3×3 and 4×4 CG1 rules and the default 3×3 and selected 4×4 CG2 rules. These are discrete Operator checks, not a Hertz continuum-limit proof.

A second Hertz-limit screen holds the classical contact radius sqrt(R delta)=0.2 fixed while increasing R and decreasing delta. It uses the same extent-8/grade-4/penalty-1e6 quarter block and default 2×2 contact rule, with four increments and serial GMRES/GAMG:

Sphere radius Indentation Quarter cells Full force Hertz force Ratio Wall time
2 0.02 13,824 6.16946 5.86081 1.05266 56.07 s
20 0.002 13,824 0.606669 0.586081 1.03513 53.87 s
20 0.002 32,768 0.600922 0.586081 1.02532 142.64 s
20 0.002 64,000 0.600599 0.586081 1.02477 342.08 s
200 0.0002 13,824 0.0605564 0.0586081 1.03324 53.20 s

At 64,000 quarter cells, the R=20 force changes by only about 0.055 percentage points from the 32,768-cell result, while still exceeding the half-space reference by 2.48%. Increasing the radius to 200 at the coarser 13,824-cell resolution likewise leaves a 3.32% excess. The small-strain limit therefore explains part, but not all, of the earlier discrepancy. The active area in the 64,000-cell R=20 run is 0.11689 versus the Hertz circle area 0.12566; finite penalty, finite-domain compliance, pressure integration and finite-strain bulk behavior remain to be separated. Continuing the same mesh sequence without changing the experiment is not justified. These are source-checkout diagnostics, not whole-solver validations or portable performance benchmarks.

Changing the bulk displacement order gives a more informative coarse-mesh screen at R=20, delta=0.002 and bulk quadrature degree 4. The rows without an extent qualifier use extent 8, grading 4; the extent-16 rows use grading 5:

Quarter cells Displacement Face rule Full force Hertz ratio Wall time
512 CG1 on Q1 geometry 2×2 0.652103 1.11265 7.62 s
512 CG1 on Q1 geometry 3×3 0.629171 1.07352 5.03 s
512 CG2 on Q1 geometry 3×3 0.603484 1.02969 20.06 s
512 CG2 on Q1 geometry 4×4 0.611064 1.04263 20.80 s
1,728 CG2 on Q1 geometry 3×3 0.600800 1.02511 69.62 s
1,728, extent 16 CG2 on Q1 geometry 3×3 0.597342 1.01921 107.05 s
4,096, extent 16 CG2 on Q1 geometry 3×3 0.594147 1.01376 308.31 s
8,000, extent 16 CG2 on Q1 geometry 3×3 0.591655 1.00951 673.28 s

The two nine-point 512-cell rows hold the face integration count fixed; quadratic displacement markedly reduces the coarse-grid force bias. Its 1,728-cell result approaches the 64,000-cell CG1 result (1.02477) at lower wall time in this environment, but both retain an approximately 2.5% reference gap. These are different discrete spaces and single-run wall times, not a portable solver speedup or evidence that the continuum Hertz limit has been reached. A 4×4 face rule does not monotonically improve the coarse CG2 result, so neither interpolation nor quadrature may be auto-selected from this screen. Bulk material, domain, penalty and contact-footprint effects still need independent limits.

The first extent-16 row increases block half-width/depth from 8 to 16 while changing grading from 4 to 5; top-center spacing stays similar (0.05905 versus 0.05610). The force ratio decreases from 1.02511 to 1.01921, consistent with finite-boundary stiffening. Because the separable mapping and remote cell shapes also change, this is a controlled-direction screen, not an isolated half-space extrapolation; the finite-domain gap remains unresolved.

The example now samples the analytic Hertz pressure on the actual reference trace and reports both its quadrature-integrated force and a weighted absolute pressure difference. This diagnostic deliberately does not identify nominal reference-area penalty pressure with Cauchy pressure. At the fixed-radius settings above, it reveals why total force alone is an unsafe promotion gate:

Quarter mesh / rule Sampled Hertz force / analytic force Nominal pressure L1 difference / analytic force
512 CG1, 2×2 1.0471 0.5198
512 CG2, 3×3 1.0194 0.8105
512 CG2, 4×4 0.9813 0.4490
1,728 CG2, extent 16, 3×3 0.9776 0.3759
1,728 CG2, extent 16, 4×4 1.0096 0.4301
4,096 CG2, extent 16, 3×3 1.0111 0.1267
8,000 CG2, extent 16, 3×3 1.0004 0.0866

The coarse CG2 3×3 total force is closest to Hertz in the table above, yet its sampled pressure-shape difference is the largest. A denser face rule improves that shape diagnostic while moving total force away on the coarse mesh. At 1,728 cells it instead improves the total force while worsening the pressure difference, so neither ordering is a robust degree-selection rule. This does not prove a defect in the law: finite penalty, finite support, a non-half-space material problem, distinct stress measures, and trace quadrature remain mixed. It does reject the tempting conclusion that a 1.9–3% force gap alone validates the contact discretization. The earlier rows are source-run diagnostics; the 8,000-cell row uses the isolated candidate wheel. Neither path supplies an external same-input validation or a portable performance claim.

The 4,096-cell row holds extent, grading, penalty, motion, and the 3×3 face rule fixed while refining the 1,728-cell mesh. Its near-contact reference spacing decreases from 0.05610 to 0.03982, the full-body force changes from 0.597342 to 0.594147 (Hertz ratio 1.01921 to 1.01376), and active reference area changes from 0.09840 to 0.12444 versus the Hertz area 0.12566. The pressure-difference reduction to 0.1267 supports a substantial spatial-resolution contribution, but a single pair of meshes is not convergence, and the finite penalty and bulk-domain limits are not yet separated. This finer source run took about 308 seconds on this Mac, versus about 107 seconds for the 1,728-cell run; these are not portable timings.

An additional 8,000-cell installed-wheel run at the same extent, grading, penalty, motion, bulk and face rules reduces the force ratio to 1.00951 and the pressure difference to 0.08664. The top-center spacing is 0.03083, and the sampled analytical Hertz force is 1.00041 times its continuous value. Yet the active reference area changes from the 4,096-cell value 0.12444 to 0.11888, moving away from the Hertz area 0.12566. Thus even this third spatial point does not justify a pressure-footprint or continuum-accuracy pass. It took 673.28 s on this Mac; accepted linear solves used 523.46 s of that time. This is an identified bottleneck, not a portable speed promise.

At the fixed 1,728-cell, extent-16, grade-five Q2 mesh with a 3×3 face rule, a three-point penalty screen gives:

Penalty Full force / Hertz Pressure L1 / Hertz force Maximum penetration
5e5 1.01444 0.3375 1.63e-5
1e6 1.01921 0.3759 not recorded in this screen
2e6 1.02153 0.4193 4.60e-6

Increasing penalty decreases penetration but does not drive the force or pressure diagnostic toward Hertz on this mesh. Choosing the softer penalty because its force ratio happens to be closer would tune against a mismatched half-space reference, not establish accuracy. A consistent zero-penetration limit still requires controlled mesh and domain refinement, and no threshold has been promoted from this screen.

Keeping the 1,728-cell, extent-16, grade-five, 1e6-penalty, CG2/3×3 problem fixed while changing the prescribed-motion path from four to eight equal increments changes final full-body force only from 0.597341748365 to 0.597341748239, and the nominal pressure L1 diagnostic only from 0.375891887091 to 0.375891886907. The eight-increment accepted trapezoidal work residual is 5.31e-7; it is a path-discretization ledger, not a general first-law balance. This final-state check suggests that increment count is not the main source of the present Hertz mismatch. It does not establish every intermediate friction/plastic path, which is a different validation problem.

The same four-increment source run now compares three distinct energy quantities: the analytical Hertz half-space displacement work (2/5) F_H delta = 0.00046886447; full-body accepted trapezoidal tool work 0.00048329004 (ratio 1.03077); and full-body accepted stored plus contact-penalty energy 0.00047528648 (ratio 1.01370). The quarter-model accepted work residual is 2.00089e-6, or four times that for the full body. The work ratio cannot be interpreted as a continuum error independent of increment size, and the recoverable-energy ratio still mixes finite domain, finite penalty and finite-strain differences with the Hertz limit. This comparison makes the boundary-work side visible without treating a trapezoidal ledger as a universal first-law proof. On the 8,000-cell installed-wheel screen, the recoverable-energy ratio decreases to 1.00766, while the four-increment tool-work ratio remains 1.02563 and the quarter-model trapezoidal residual is 2.10628e-6. The different movements of these energy channels reinforce the need to separate spatial error from path-work quadrature before setting a threshold.

Remaining gate

Run independently controlled mesh and domain-size refinement with a resolved contact patch, then penalty and increment sensitivity. Compare force, pressure footprint, penetration, reactions, tool work and plastic state against either the valid Hertz limit or a frozen external FE result. Do not turn the current force ratio into a release tolerance. Active noncoplanar triangle seams/free edges, distributed implicit assembly and broader friction remain separate work; the bounded serial multi-tool and fixed-plane friction candidates below do not close those general gates.

The benchmark now also reports reference-area activity, maximum sampled contact radius, maximum nominal penalty pressure, and the actual near-contact mesh spacing. An optional separable coordinate grading resolves the patch without uniformly refining the remote block; it changes only the example mesh, never the solver's geometry or quadrature. The independent force limit must be judged together with these spatial diagnostics, rather than as a single scalar. COMSOL's contact-fatigue example likewise identifies the contact zone and nearby subsurface as requiring much finer resolution than the remote structure.

Local-parity acceptance map

The deliberately bounded target is 3D solid versus prescribed rigid tools, not general contact. This record separates code gates from evidence gates:

Gate Current evidence Outstanding acceptance
Comparable problem Finite-block sphere setup is fixed and reproducible in source Freeze an independently run Abaqus deck with the same block, law, mesh, boundary path and penalty; no numerical Abaqus reference has been supplied
Existing contact composition Merged PR #116 has reviewed smooth/facet rigid geometry, affine MPC, moving-tool work and rollback; serial focused regressions pass Installed-package and post-merge release-gate run for the new changes
Facet features Open/nonconvex triangle tools still reject active noncoplanar ambiguity. A separately certified closed-convex tool now has a continuous interior support-plane gap, one deterministic limiting derivative at a crease, and an elastic/plastic finite-J2 feature-switch path with restart General nonconvex edge/vertex laws, external sharp-feature comparison, scalable search and mesh/penalty convergence; no normal averaging
Total implicit Operator Contact/strong-BC overlap now uses native contact lifting, identity rows, a complete fixed-old-state total-tangent difference check, failed-increment rollback and same-DOF checkpoint recovery. On the separate affine-periodic path, a certified convex tool's flat face matches the analytical-plane solution and contact potential, and the complete J2/contact residual satisfies TᵀR≈0 after MPC reduction Add an independently constructed residual/force check and expand beyond the current bounded serial discretizations
Multi-tool implicit Ordinary serial strong-BC J2 composes distinct rigid pairs: summed residual/tangent, independent moving-pair work, and atomic failed-increment/checkpoint recovery. Opposing moving frictionless planes check force balance. Two fixed frictional planes recover the equivalent single-law response while preserving separate State/dual/checkpoint identities; spatially opposed fixed frictional planes also check both force resultants and interrupted restart. Two separately moving convex frictional tools now keep pair-specific work, raw force reconstruction, friction State and restart identity; an injected second-pair commit failure restores both accepted states. The separate affine-periodic J2 route now admits named frictionless pairs: equivalent split-tool and spatially opposed tools check total MPC reduction, independent moving-tool work, failed projection rollback and interrupted restart Mixed curved/open-facet friction beyond the reviewed face branch, affine-periodic friction, broader plastic branch switching, external structural comparison and finer mesh/penalty checks
Implicit friction Fixed, translating and rotating analytical planes are admitted on serial strong-BC finite-J2. A certified closed-convex triangle solid shares the selected planar-branch tangent: a rotating single-face path matches the equivalent analytical plane, and one plastic noncoplanar feature-switch path preserves tangential State and restart. A reviewed open triangular tool now admits active face-interior planar branches: a coplanar two-triangle tool matches the analytical plane's force/tangent under translation and its solution, work and interrupted restart under rotation; active free edges still fail atomically. Stable-ID State, rigid pose, separate Coulomb dissipation and numerical separation release, failed-finalization rollback, serial checkpoint and result histories remain common Mesh/penalty/increment and external sharp-friction comparisons; active noncoplanar edges/vertices, curved/MPC/MPI friction remain unsupported; no general friction claim
Material/time composition Elastic-regime J2 sphere and plastic J2 facet cases complete; the convex sharp-feature path also evolves plastic State and restarts before a feature switch. A two-region finite-J2 material map now shares a moving convex frictional Step; uninterrupted and checkpoint-resumed displacement, quadrature State and dissipation agree. Explicit body contact has its own stable-step contract A genuinely different second history material on the same implicit seam and explicit contact stability under evolving state; two parameter sets of J2 are not two constitutive laws
Independent accuracy and speed Full/quarter symmetry, mesh, penalty, interpolation-order, face-quadrature and increment screens exist; the fixed-contact-radius screen reaches a 0.95% force gap and a 0.0866 nominal-pressure L1 difference on an installed-wheel 8,000-cell Q2 quarter model at extent 16. The active area is nonmonotone against the Hertz area; this remains directional, not a converged limit. The 673 s run spends 523 s in accepted linear solves. The sharp plastic screen below is also mesh/penalty sensitive Controlled far-boundary, graded mesh, pressure-footprint and penalty limits; external FE comparison and repeatable whole-step performance trials
Product route Ordinary model.step() and isolated installed-wheel closed-convex and rotating-plane-friction examples have run; final unsmoothed reference-trace pressure/penetration scalars are available; a serial accepted-point NPZ/JSON archive keeps pair identity, signed gap, reference weights, pressure, friction traction and accepted tangential state without smoothing. An optional vertex-only VTP view puts the same samples at current slave positions; an independent VTK reader checks coordinates and pressure without interpolating. The public capability matrix is explicit Broad installed-platform and candidate-bound release gate; the VTP point cloud is not a continuous contact-surface pressure map

The output candidate at commit e054df76 was built as an isolated native Apple Silicon wheel, installed under a temporary path and imported from that path while the process ran outside the source checkout. The ordinary moving closed-convex-friction example completed four increments, reported tool work 2.2357364364 and friction dissipation 1.1808584038, and wrote accepted_contact_points.npz plus its schema JSON. The archive opens with numpy.load(..., allow_pickle=False) and contains six raw quadrature points for this tiny fixture. This is an installed-path/format smoke, not a Windows/Linux test or external contact validation.

A second isolated native Apple Silicon wheel built from the later source candidate at 25f3c7d1 was imported outside the checkout and ran the ordinary closed-convex moving-friction model.step() example. Its four accepted increments reproduce tool work 2.2357364364 and Coulomb loss 1.1808584038; the new optional accepted_contact_points.vtp has six vertex cells and an independent VTK reader reproduces the unsmoothed NPZ pressure values. This tests installed packaging and visualization syntax, not physical convergence, non-Mac portability or a release promotion. One bounded two-rank MPI regression group additionally confirms that the closed-convex friction route fails at capability/Step construction rather than silently entering an unsupported distributed implicit solve.

From a fresh archive of source commit c65294cc, a native Apple Silicon wheel and matching source distribution were built in an isolated temporary packaging directory. The installed wheel SHA-256 is 3e1f3d1528e5cce19a5ca204297820dd069f19ede50311a56a395216186b5db1. The standard release_gate.py --dist ... --smoke passed its installed-wheel, artifact-identity, flagship, entrypoint and project-template checks without an editable import. An additional installed-path model.step() sphere contact smoke with the explicit 3×3 Q1 face rule completed two increments; its eight-cell force is deliberately not an accuracy result. This is a Mac candidate packaging check, not a Windows/Linux gate, a two-rank release smoke, an Abaqus comparison, or a promotion of general implicit contact.

From a second fresh archive of source commit 60215b2d, a native Apple Silicon wheel and matching source distribution were built and passed the standard installed-wheel release_gate.py --smoke artifact, flagship, entrypoint and project-template checks. The wheel SHA-256 is f99526ef3937b18267357e9f4a13c772da8a79b6b5181eaf120d673e3cd8a80b. Outside the source checkout, the installed package completed an ordinary CG2 finite-J2 sphere-contact Step with the selected degree-seven/16-point face rule. The resulting force on its deliberately coarse two-cell fixture is not an accuracy result. This gate proves installed packaging and bounded execution on this Mac, not Linux/Windows, external comparison, general friction, or an actual release promotion.

The small CG2 hexahedral tool-contact fixture exposed a different failure: with the former inherited bulk degree-two default it accepted the first of two increments, then PETSc reported a failed preconditioner on the second. Holding geometry, material, contact face rule and load path fixed while raising only the bulk degree to four completed both increments. The ordinary Step now selects degree four for CG2 hexahedra, and explicit lower bulk degrees are rejected at capability preflight. This removes one known underintegrated path; it is not a general proof of Q2 contact stability.

At source commit 9965dcce, a fresh Git archive built a native macOS arm64 wheel and sdist. The wheel SHA-256 is be28b2830c2ed17bd3068357325fe5366d553517458c096b74da2817b18bb745. The first standard platform gate ran from the archive directory and correctly returned failed because that directory had no clean Git identity, even though the installed-wheel smoke passed. Rerunning the same wheel from the clean source checkout recorded that exact commit, source_dirty=false, and passed release_gate.py --dist ... --smoke --require-platform macos. Separately importing the wheel from a temporary install target outside the checkout completed the CG2 sphere-contact Step with default bulk degree four and face degree four. A bounded two-rank regression gave 7 passes per rank for the explicit MPI contact refusal. This is a macOS candidate gate, not a distributed implicit implementation or cross-platform release approval.

At source commit 2287ac05, a fresh Git archive produced a native Apple Silicon wheel with SHA-256 5c3bcc1752b56792a10d16e73ae6d7aa9f39c37cd6b32c900a5c194d4cda6100. The standard installed-wheel macOS gate reported status=passed, source_dirty=false, and the exact source commit. Importing this wheel from an isolated target outside the checkout passed both open/closed vertex-fan regressions and completed the ordinary CG2 sphere-contact Step with default bulk degree four and selected 4×4 face rule. The subsequent 8,000-cell fixed-physics Hertz screen above also ran from this installed wheel. These are local Mac execution and scientific-direction screens, not external same-input validation, Windows/Linux qualification, or general-contact promotion.

Bounded rigid-mode GAMG improvement

The installed-wheel 8,000-cell Q2 Hertz screen above spent 523.46 s of its 673.28 s wall time in accepted linear solves. For the existing opt-in GMRES/GAMG policy on the serial native rigid-contact finite-J2 displacement route, the matrix now carries a PETSc rigid-body near nullspace formed from the reference displacement coordinates. This supplies GAMG with six elastic low-energy modes; it does not alter the contact law, material update, assembled residual, strong constraints, checkpoint identity or the converged solution. Incompatible matrix layouts leave the old preconditioner untouched. The actual attachment is reported as linear_preconditioner_near_nullspace in the Step and example summaries. The construction follows the DOLFINx elasticity GAMG example and PETSc near-nullspace contract.

On the same machine and fixed 1,728-cell quarter-block inputs (E=1000, nu=0.3, R=20, delta=0.002, extent 16, grading 5, penalty 1e6, CG2, bulk degree 4, face degree 5, four increments), the earlier isolated-wheel baseline gave full force 0.5973417483646, 5,352 linear iterations, 73.33 s accepted linear solves and 108.26 s wall time. The integrated source route gave 0.5973417483273, 4,094 iterations, 62.00 s linear solves and 94.37 s wall time. The active reference area was identical and the nominal pressure L1 ratio changed below 1e-10. A separate 8,000-cell installed-wheel monkeypatch prototype reduced single-run wall time from 673.28 s to 571.99 s with the same force to about 4e-11.

The integrated code at bfdbe26e was then built from a clean Git archive as a native Apple Silicon wheel (SHA-256 b3d7a04a8a62a45065f1a1f30f1fdb75dc9248618f415eccd6e540ecead2940a). The standard installed-wheel macOS gate passed with source_dirty=false and that exact source commit. An ordinary CG2 rigid-sphere model.step() run from the isolated install target also completed. On the same fixed 8,000-cell inputs, that installed integrated wheel gave full force 0.5916554842398 versus the pre-change 0.5916554841966, 4,977 versus 7,169 linear iterations, 410.86 s versus 523.46 s accepted linear solves and 553.43 s versus 673.28 s single-run wall time. The active area, maximum penetration, pressure L1 metric and accepted energy ratios agree to reported precision. The pressure L1 metric remains 0.08664, however, so this is not a Hertz convergence certificate. A separate plastic J2 test crosses a noncoplanar convex-tool feature: direct LU and GMRES/GAMG agree on displacement, committed plastic State and tool work, and the iterative path survives checkpoint/restart. These are one-machine timings, not a portable speedup promise. There is no new global solver default, no automatic geometry/penalty change and no claim that an elastic near-nullspace guarantees convergence in a strongly plastic or more complex contact state.

The Abaqus/Standard general-contact description uses separate surface-, edge- and vertex-contact formulations. This supports treating sharp features as a separate numerical law, not treating the current selected-facet normal as a differentiable smooth surface. No row in this map is promoted merely because a single visual plot or Newton solve looks plausible.

Next implementation decisions

The implicit friction route reuses the accepted, stable-ID TangentialContactState and return map used by explicit contact. At fixed old State and a current prescribed analytical-plane pose, the tangent must include both the pressure cap and the tangent projection of the transported master material point. Writing P=I-n⊗n and old/current rigid rotations as R_old, R_new, that relative tangential-increment derivative is A=P R_old R_new^T P. Thus a stationary or purely translating plane has A=P; replacing A with P for a rotating tool would give the wrong Newton matrix even though the return-map value looked plausible. The current-pose projection is checked against the supplied pose to reject a stale trial. The slip branch also differentiates mu p and slip direction. At opening and the exact stick/slip cap the chosen semismooth derivative is not a classical two-sided derivative.

implicit_friction.py maps this point derivative through the reference trace to the normal-plus-friction Newton Operator. Solver-neutral and real DOLFINx stick/slip direction checks pass for fixed, translating and rotating planar tools, including a nonzero accepted old rotation. A serial strong-BC finite-J2 rotation path has nonzero generalized moment, accepted tool work, positive Coulomb dissipation and interrupted restart agreement. Multiple independently named planar pairs share the same Procedure transaction; a second-pair commit fault restores both accepted states. Curved, faceted, affine-MPC and MPI friction remain rejected. The accepted energy record reports Coulomb dissipation separately from recoverable tangential penalty energy and numerical separation release; its trapezoidal residual is a diagnostic, not a universal first-law proof. For the dimensionless prescribed 0.15-radian plane rotation and 0.4 tangent translation, a 2/4/8/16-increment screen gives tool work 2.3731/2.2357/2.2688/2.2749 and friction loss 1.7735/1.1809/1.0495/0.9835. The signed mechanical-energy residual changes from -0.5847 to -0.0628 over that scan. The finite-increment values are not yet converged and no physical-unit accuracy claim follows from this screen. The sliding tangent includes a pressure-cap cross term and is generally non-symmetric. Public admission rejects symmetry-only linear solver selections such as CG or Cholesky; the default direct LU route is accepted. At the exact stick/slip cap or opening the selected semismooth branch is not a classical two-sided derivative. This narrow global solve and return-map agreement do not certify a broad industrial friction capability. The FEBio frictional-contact formulation uses an explicit stick/slip return map and verifies patch, Hertz and large-sliding cases; its numerical results are not AgentFEM evidence.

The source example finite_j2_fixed_plane_friction.py deliberately crosses the slip cap under combined normal and shear loading. For 4, 16, 64 and 128 accepted load increments on its one-cell-per-axis fixture, cumulative Coulomb dissipation is 0.18695, 0.12040, 0.10915, 0.10762; the trapezoidal energy residual is -0.07612, -0.01356, -0.00305, -0.00149. The full contact force changes by less than 0.05% between 16 and 128 increments. These dimensionless numbers document strong coarse-step dissipation sensitivity and apparent path refinement, not mesh or external benchmark convergence. A weaker synchronous shear load of 15 showed vanishing slip upon refinement; it was discarded as a misleading sliding fixture. No exact first-law or general friction claim follows from this screen.

An isolated native 0.4.3 wheel built from source commit 60bb91e1 was installed into a temporary directory and imported from that directory, outside the source checkout. The ordinary model.step() examples for the moving, rotating Coulomb plane (four increments) and closed convex tool (two cells, sixteen increments) both completed. The plane reported tool work 2.23573644 and Coulomb dissipation 1.18085840; the convex example reported right-end displacement -0.0009990007775003966 versus its narrow analytical bar value -0.0009990007774995924. This proves the local installed-package path for these examples, not Windows/Linux portability, curved friction accuracy or a release-grade scientific validation.

Purely translating planar Coulomb tool

The same serial strong-BC J2 Operator now accepts a prescribed translation of an analytical frictional plane on the shared normalized Step coordinate; rotation still fails closed. relative_contact_displacement_increment tracks the current projected master material point back to the accepted rigid pose, so tangential tool motion is not erased by subtracting gap vectors. With old State and current pose fixed, the planar stick/slip tangent remains the same derivative with respect to slave displacement. The accepted pose factor is stored in point-keyed friction kinematics; corrupted pose-factor checkpoints are rejected before replacing State. Tool work uses the total normal-plus- friction resultant and moment, while Coulomb dissipation and recoverable tangential penalty remain separate channels.

Focused DOLFINx tests check translated normal and tangential forces, fixed-old-State directional tangents, rotating-tool rejection, an ordinary moving-tool J2 Step, accepted tool work and interrupted restart. The public two-increment fixture deliberately crosses the Coulomb cap; a weaker path reached the cap without sliding and was not used to claim dissipative validation. This bounded translation path is not a curved-tool, rotating- frame, general finite-sliding or affine-MPC friction solver.

The same 16-increment fixture completed from an isolated locally built macOS wheel with no source checkout on PYTHONPATH; agentfem.__file__ resolved inside the isolated wheel target. This is an installed-package smoke, not a Windows/Linux or release qualification. An independently calculated, displacement-controlled FE patch checks the complete contact integration and public reaction field on an area-one face: with normal penalty 500, penetration 0.05, tangential penalty 100, and mu=0.4, the normal resultant is -25. At tangential displacement 0.02 the force is -2 and sliding loss zero; at 0.20 the force is capped at -10, cumulative Coulomb loss is 1.0, and recoverable total contact potential is 1.125. Both 3D J2/strong-BC solves match these constructed values within 1e-10 absolute tolerance. This is a bounded patch identity, not an independent large-sliding or forming validation. This patch also exposed and fixed a result-contract defect: the public nodal RF field previously assembled only the UFL material/external form while a native rigid-pair contact Operator was added separately to Newton. It now uses the same full-residual helper as prescribed-work integration. The patch checks that strong reactions plus the independently reported contact resultant vanish without external traction; native contact is no longer missing from RF. Legacy UFL contact remains assembled once.

In-memory increment rollback now retains the existing read-only, point-keyed friction records instead of converting every field to nested Python lists. Only a written checkpoint serializes JSON-safe values. A targeted conversion microbenchmark with 50,000 points measured about 8.0 ms per former list-producing snapshot on this machine; the new runtime snapshot stores two object references. This removes that conversion from each accepted attempt, but is not a measured whole-solve acceleration and must not be marketed as one. Single- and multi-pair failure rollback still pass after the change.

For multiple rigid tools, the ordinary serial strong-BC J2 Step now adds independently named frictionless pair Operators. Each pair retains its own projection identity and prescribed pose, while one owning Procedure accepts or rolls back the increment. SimulationResult carries separate duals, nodal reaction distributions, quadrature diagnostics and accepted cumulative tool work. The full-space contact tangent sum is compared with a central difference in an elastic-regime J2 fixture. Two independently moving tools also pass interrupted/restarted path equivalence; malformed per-tool work is rejected on checkpoint load. A separate opposed planar-die test checks that both tools remain active, force directions oppose, support-inclusive force balance closes and result metadata retains both identities. This is not yet evidence for arbitrary shaped dies, mixed sharp/curved branch switches, moving/curved friction, affine MPC multi-pair reduction or MPI.

examples/finite_j2_opposing_rigid_tools.py is a compact ordinary model.step() example. A locally built candidate wheel was installed with --no-deps --no-index --target into an isolated directory and that example completed from outside the source checkout; the imported package resolved to the isolated wheel target. The result carried four accepted increments, two separate tool resultants and accepted works, and an available constraint work balance. This is a local Apple Silicon installed-wheel smoke, not a Windows/Linux platform or final release acceptance.

A targeted topology audit found an active-vertex loophole in the bounded planar-facet admission: a selected triangle can have two coplanar neighbors at its vertex while the full vertex fan is still open elsewhere. With same-point/same-normal facet ties and a middle facet having the lowest ID, the old two-edge check admitted that free boundary vertex. A failing regression reproduced it before the fix. The Operator now checks every incident facet normal and every edge containing an active vertex; a free or noncoplanar member rejects the trial, while a closed coplanar fan remains admissible. This strengthens fail-closed semantics, not a new edge/vertex contact law. Forty-two focused solver-neutral/DOLFINx facet and friction tests pass after the correction.

For a truly sharp noncoplanar tool edge, a bounded feature-contact law needs face/edge/vertex classification plus an admissible-side rule and a semismooth/generalized derivative. A declared physical fillet instead uses the smooth-surface path and must keep its radius in the model identity. The FEBio paper reports that a slight radius was needed to solve one sharp-corner example; this is motivation for an explicit geometry choice, not permission to round an imported mesh behind the user's back.

For an ideal noncoplanar sharp edge, one selected-facet signed plane gap does not have a continuous normal derivative across the edge at finite penetration. Requiring a continuously rotating force there without a declared fillet or a feature-specific contact potential would silently change the geometry or the law. The proper sharp-feature gate is: identify face/edge/vertex ownership, define one potential with the correct admissible side, check potential continuity and one-sided or generalized residual derivatives on each branch, then verify stable branch switching and work. For a physically rounded die, continuity should instead be tested against the declared smooth radius. The current fail-closed error is preferable to choosing an arbitrary adjacent facet normal, but it is not the finished feature-contact algorithm.