[jmartfrut]

∇∇ was compiled by the kernels in cell and boundary integrals only (0.2.1); a skeleton term with Hessians — jump(∇∇(u)) ⊙ jump(∇∇(v)), mean(Δ(u)), … — sent that skeleton to the Gridap fallback ("second derivatives (∇∇) on skeleton integrals are not supported yet"). They are now compiled by the facet kernels. The trigger is the third-medium regularisation of the gear contact problem (G1): β ∫ ∇∇u ⊙ ∇∇v dΩc was already in the kernels, and the only term left to Gridap was c₂ β h_F ∫ ⟦∇∇u⟧ : ⟦∇∇v⟧ on the skeleton of the third medium.

What changes

• SkeletonGroup no longer rejects leaves of deriv = 2. Everything else was already in place from the two earlier pieces of work: the probes and the slot layout give Hessian leaves a side (plus/minus) like any other leaf, and the facet tables are ordinary FieldTable/GeoTable at the points of each side, which carry the reference Hessians (ddNt, ddM).
• Push-forward per side. The Hessian of each side uses the Jt⁻¹ and the curvature Φ of that side's cell map at its own mapped points (side_geometry: point_curvature of geo for the plus cell and of geom for the minus cell, evaluated only in classes with Hessians — hess_use(::SkeletonGroup) —, zero on affine cells). It enters the same three places as in the volume kernels: the values of the unknown (field_values_sides with one Φ per (side, field)), the physical weights of test/trial functions (point_weights! with the Hessian columns, only for the fields with Hessian slots) and the pull-back of the residual coefficients (contract_residual!, with the curvature term added to the gradient coefficients). So ⟦∇∇u⟧ is exact on non-affine quadrilaterals and hexahedra too, where Gridap's ∇∇ is not defined.
• Matrix-free Jacobian. The two GroupData of a skeleton class store the curvature of their side (Phi, only with Hessians on non-affine cells); the dense apply of the facets reads it like the volume groups. Weight buffers of the skeleton classes are sized with wcols (Hessian columns).
• No change to the sparsity pattern, colouring, facet positions or update (the signature already carried deriv and side).

KernelFunctional with skeleton integrals (second commit)

The energy of the same problem (for the Armijo line search of InjectivityLS) has the potentials of both jump penalties, ½ ∫ c_J ⟦∇u⟧:⟦∇u⟧ + ½ ∫ c_J2 ⟦∇∇u⟧:⟦∇∇u⟧ on the skeleton. KernelFunctional did not evaluate skeleton integrals: compiling the skeleton failed (the facet colouring needs the test rows of an operator) and those terms were summed by the Gridap fallback. Now:

• compile_groups(...; facet_colors = false) builds the skeleton data with every facet in one colour (skeleton_data(...; colored = false)): a reduction has no scatter and no races.
• functional_cells for a SkeletonGroup: both cells gathered, side_geometry and field_values_sides (values, gradients and Hessians of both sides) at the facet points, the integrand times the facet measure (:sum) or the extremum over the points (:min/:max).
• KernelFunctional keeps its skeleton data (skel), reused by update; aux_data of a skeleton class works for functionals too.

Tests (test/solvers/test_functionals.jl): the energy of the third-medium regularisation (jumps of ∇u and ∇∇u on the skeleton of a subdomain, facet size as data) plus a mean term on the whole skeleton, against Gridap's lazy sum at 1e-12 on triangles, Cartesian and non-affine quadrilaterals (hessian_nonaffine); update with another β (groups rebound, tables and skeleton data reused); one thread vs default; the gradient of the energy equals the residual of the operator (finite differences); :min/:max of a facet quantity against Gridap; no allocations in the facet loop.

Still going to the fallback: ∇∇ of Raviart–Thomas fields; third derivatives. Not supported (unchanged): BlockJacobi with skeleton integrals.

Tests (test/test_hessian.jl; the "fallback of ∇∇ on a skeleton" test is replaced)

• C⁰ interior penalty (scalar): γ/h_F ⟦∂ₙu⟧⟦∂ₙv⟧ − {Δu}⟦∂ₙv⟧ − {Δv}⟦∂ₙu⟧ plus a nonlinear Hessian-jump term 0.2 h_F (1 + {Δu}²) ⟦∇∇u⟧:⟦∇∇v⟧ (dual-number path) and the volume ∇∇v ⊙ ∇∇u.
• G1 regularisation (vector): St Venant–Kirchhoff on the mesh, β ∫ ∇∇u ⊙ ∇∇v on a subdomain and σβ/h_F ⟦∇u⟧:⟦∇v⟧ + c₂βh_F ⟦∇∇u⟧:⟦∇∇v⟧ on the skeleton of that subdomain (composite triangulation, facet size as data).
• Against Gridap at 1e-13 (residual and assembled Jacobian; matrix-free product at 1e-11) on triangles P2/P3, Cartesian quadrilaterals Q2 and tetrahedra P2; linear maps on/off, 1 thread vs default, zero allocations, no fallback.
• Non-affine quadrilaterals and hexahedra Q2 against the same forms written with hessian_nonaffine (HyperFEM's push-forward, test/common.jl) at 1e-12; the curvature of both sides stored in the matrix-free data.
• A quadratic field (constant Hessian) has zero Hessian jump and the exact mean Hessian on every facet (Q2 Cartesian, Q2 non-affine, P2).
• Multi-field with Hessians of both fields on the skeleton, Hessian only in the Jacobian, update with another coefficient (same compiled kernels and tables) and the matrix-free Jacobian rebound.

Measurements

TMC_gears migration/check_forms_G1.jl (our measurement: Mac, Julia 1.13, 18 threads; mesh SpurGear_z25_nt1_hc0.02, P2, 33 330 dofs, HESS=nabla), with this branch:

• no fallback warning: the whole weak form (volume, ⟦∇u⟧ and ⟦∇∇u⟧ on the skeleton) and the energy are evaluated by the kernels;
• check against Gridap: residual 1.4e-15, Jacobian 1.4e-16 relative;
• residual + Jacobian 0.024 s (Gridap 0.141 s); energy 0.0035 s (Gridap 0.017 s), difference 1.1e-15 relative.

Test suite on the Mac before the second commit: 72 471 passed, 4 failed (the matrix-free curvature expectation for Cartesian quadrilaterals in the new test, fixed: their Q1 reference Hessians do not vanish); test_hessian.jl after the fix: 3 891 passed; test_functionals.jl with the second commit: 51 passed.