I had Claude generate a mortar contact implementation and then asked what was awkward about it. This is what it wrote. Putting it up in case it could be useful.

What Ferrite needs to provide for contact mechanics on top

Ferrite's role is to provide general FEM building blocks; algorithms like contact belong in satellite packages (a future FerriteContact.jl / FerriteMortar.jl). This document derives, from the requirements of a proper contact solver, the primitives such a package needs from Ferrite core — what exists today, what is missing, and what the smallest useful core additions would be.

1. What a contact solver has to do

Any serious contact package ends up with this pipeline, independent of the specific method:

1. Search: find candidate surface pairs (broad phase), then closest-point / ray projections between non-matching surfaces (narrow phase), in the current (deformed) configuration.
2. Pairing/integration domains: node-to-segment pairs (NTS), or mortar segments (2D) / clipped polygons (3D), or quadrature-point-wise projections (Gauss-point-to-surface, Nitsche).
3. Two-sided evaluation: evaluate trial/test traces of both bodies (and possibly a multiplier space) at custom, per-pair quadrature points on facets.
4. Enforcement: penalty, Lagrange multipliers (+ active set / semismooth Newton / PDAS), augmented Lagrangian (Uzawa), dual mortar with condensation, or Nitsche. Friction adds a tangential complementarity problem per point.
5. System assembly: couple dofs of two unrelated cells (and multiplier dofs) into one matrix; boundary conditions must compose with the coupling terms.
6. Nonlinear loop: re-search / re-pair / re-assemble per Newton iteration or load step; the sparsity footprint changes as the active set and pairing evolve.

Notably, the enforcement choice decides how demanding the package is on core:

• Penalty / Nitsche need no new unknowns at all — only items 1–3 and 5–6.
• Multiplier methods (what one wants for robust, penetration-free contact) additionally need unknowns that live on a surface, not in a cell field — fully integrated with constraints, sparsity, and assembly. This is by far the largest core gap.

2. Gap A — fields on surfaces (multiplier/trace fields)

Needed by: Lagrange-multiplier and augmented-Lagrangian contact, mortar tying, and generally any mixed method with interface unknowns. The contact pressure is a trace field: dofs attached to the nodes (continuous variants) or facets (dual/discontinuous variants) of a facetset, not to cells.

Today: DofHandler distributes dofs for cell fields only, and everything downstream is scoped to those dofs:

• ConstraintHandler is sized to ndofs(dh) (close! builds isconstrained[1:ndofs], src/Dofs/ConstraintHandler.jl:289). Multiplier unknowns that are not DofHandler dofs cannot be Dirichlet-fixed (which is exactly what an inactive multiplier in an active-set method is), cannot appear in an AffineConstraint, and are invisible to apply!/apply_zero!.
• Sparsity patterns, assemblers, renumber!, and export are likewise scoped to the DofHandler's dofs, so multiplier couplings and their boundary-condition interactions all have to be re-implemented in user space (including the correction of coupling blocks for prescribed dofs — and doing that by hand is only straightforward for pure Dirichlet conditions, see Gap F).

The seemingly natural workaround — a second Grid of Line cells over the surface with its own DofHandler — is blocked twice: two DofHandlers cannot share a global numbering or a sparsity pattern/assembly, and ExclusiveTopology rejects embedded cells outright (src/Grid/topology.jl:56); dof distribution for embedded cells is itself flagged as incomplete (src/Dofs/DofHandler.jl:462-464).

What core should provide: facet-attached fields — trace fields as first-class citizens, e.g. add!(sdh, :λ, ip, facetset) distributing dofs on the sub-entities (vertices/edges/faces) of a facetset. This gives multiplier dofs celldofs-style access from the adjacent cells (correct sparsity for free), makes Dirichlet on multipliers natural, and also serves other mixed interface methods: hybridization and skeleton/static-condensation approaches, interface enrichment, and weak imposition of interface conditions in general.

The interface should be interpolation-based (the user picks the trace space) rather than hardcoding a nodal P1 space: dual (biorthogonal) mortar bases make the multiplier block diagonal so multipliers can be condensed out and the active set handled purely primal — the industry-standard mortar-contact formulation — and that only works if the package can choose the dof layout of the trace field.

This is the one gap that is a genuine design effort rather than an incremental API addition. There is adjacent prior work in the global-fields direction (#1395, and issues #1009/#596); a facet-field design should be developed with that machinery in mind so dof bookkeeping, sparsity, and constraint handling are solved once.

3. Gap B — re-initializing FEValues at custom quadrature points

Needed by: every contact discretization. Mortar integrates on segments (points differ per segment on both sides); NTS/Gauss-point-to-surface and Nitsche evaluate at projection points that change every Newton iteration; 3D mortar integrates on triangulated clipped polygons with a variable number of points per facet.

Today: arbitrary-point rules are constructible (QuadratureRule{RefLine}(weights, points), src/Quadrature/quadrature.jl:53-62) and FacetQuadratureRule accepts per-facet custom rules, but re-pointing an existing FEValues is internals-only: the idiom used inside InterfaceValues.reinit! — mutate getpoints(fv.fqr, facet) in place, call precompute_values! on get_fun_values/get_geo_mapping, then reinit! (src/FEValues/InterfaceValues.jl:120-145). Array sizes are frozen at construction (src/FEValues/FacetValues.jl:54-59), getnquadpoints is derived from the stored weights, and Base.copy(::QuadratureRule) is an aliasing no-op (src/Quadrature/quadrature.jl:324) that the unwary will trip over.

What core should provide: a supported, documented

reinit!(fv::FacetValues, cell, x, facet_nr, qr::QuadratureRule)   # and CellValues analog

that re-evaluates at the given reference points (resizing or bounding the internal storage), plus a documented statement of the weight conventions (a custom facet rule's weights are in the facet parametrization, including the refshape-dependent scaling — e.g. the ½ applied for RefTriangle facets, src/Quadrature/quadrature.jl:229). This one API unlocks mortar, XFEM, cut-cell/CutFEM, embedded interfaces, and transfer operators as external packages. It is the highest-leverage small addition.

4. Gap C — facet parametrization and geometry as public API

Needed by: any surface-coupling code that must relate facet-local coordinates, cell reference coordinates, and physical coordinates on both sides of a pair.

Today: everything required exists but is internal and undocumented: facet_to_element_transformation / element_to_facet_transformation (src/FEValues/facet_integrals.jl:103-152), weighted_normal, facets(cell) node tuples (src/Grid/grid.jl:132-134), PathOrientationInfo / SurfaceOrientationInfo. Surface-coupling code additionally needs orientation invariants — for example, that facet_to_element_transformation(ξ = -1) lands on the first node of facets(cell)[f] — which hold today but are stated nowhere, so external code cannot rely on them across Ferrite versions.

What core should provide: promote to documented API (not necessarily exported): the two transformations with an explicit orientation contract, weighted_normal, and small conveniences — facet_nodes(grid, fi::FacetIndex) and facet_coordinates(grid, fi). Cheap, zero-risk, and it converts a version-brittle package into a stable one.

5. Gap D — surface point search (closest-point projection)

Needed by: the narrow phase; in large deformation this runs every Newton iteration, and in 3D a brute-force all-pairs loop stops being viable.

Today: PointEvalHandler has the right ingredients — KDTree/knn broad phase, node→cell adjacency (_get_node_cell_map), and a damped-Newton local solve find_local_coordinate that is already dimension-generic (it uses the pseudo-inverse calculate_Jinv for embedded Jacobians, src/FEValues/FunctionValues.jl:161-169) — but the outer constructor demands rdim == sdim (src/PointEvalHandler.jl:58), so none of it is usable for "closest point on a facetset". There is also no way to run the search against a deformed configuration other than materializing a deformed grid.

What core should provide: a facetset point query, e.g. closest_facet_point(grid_or_cache, facetset, x_query; x = coords) returning (FacetIndex, ξ_facet, distance), built on the existing KDTree + Newton pieces, and accepting a coordinate vector (deformed configuration) instead of only the stored grid coordinates. Even just unlocking find_local_coordinate plus a documented adjacency helper would let packages assemble the rest.

6. Gap E — two-sided evaluation for non-conforming pairs

Needed by: everything that pairs facets of two different bodies.

Today: InterfaceIterator/InterfaceValues assume a conforming shared facet: iteration is driven by ExclusiveTopology and hard-errors on multi-neighbor facets (src/iterators.jl:374), and the quadrature transfer between sides is a pure orientation flip (transform_interface_points!, 2D case src/FEValues/InterfaceValues.jl:571-582). A contact package must drive InterfaceCache-like machinery manually with its own pair list.

What core should provide: nothing structural — Gap B already covers the evaluation side — but a documented "bring your own pairs" path (InterfaceCache reinit from two explicit (cell, facet) pairs without a topology) would make such packages idiomatic rather than hand-rolled. If Gap B lands for FacetValues, an InterfaceValues-with-custom-points variant is a natural follow-up: two-sided values on a shared set of physical points with per-side reference coordinates, which is exactly the object mortar and Nitsche methods integrate with.

7. Gap F — constraint composition with coupling blocks

Needed by: real models, which combine contact with symmetry planes, MPCs, periodic BCs, rigid drivers.

Today: correcting a hand-assembled coupling matrix for prescribed dofs is easy only for pure Dirichlet conditions. Ferrite condenses AffineConstraints as C'KC (_condense!, src/Dofs/ConstraintHandler.jl:698); the analogous transformation for a rectangular coupling block (B → BC, constant part into the gap/rhs) has no public entry point, and nested constraints are rejected (src/Dofs/ConstraintHandler.jl:323) — relevant because "multiplier-coupled dof that is also periodic" is exactly a nested case.

What core should provide: a utility to apply a ConstraintHandler's condensation to a user-supplied rectangular block and rhs — i.e. expose the C action (apply_constraint_matrix / condense_rhs!-style API). With Gap A in place, much of this falls out automatically because multiplier dofs become ordinary dofs that the existing condensation already covers.

8. Gap G — evolving sparsity over the nonlinear loop

Needed by: active-set methods (multiplier rows switch between identity and coupling), and large-deformation contact (the pairing — hence the coupling nonzeros — changes between load steps).

Today: patterns are static at allocate_matrix time; the assembler errors on missing entries. A workable strategy already exists: allocate the union pattern over all broad-phase candidate pairs once, and assemble zeros where inactive. What is missing is only convenience: an add_interface_entries!-style helper that adds cross-coupling entries for an arbitrary list of facet pairs (and, with Gap A, multiplier–cell couplings) into a SparsityPattern.

What core should provide: add_entries!(sp, dh, pairs::Vector{NTuple{2, FacetIndex}}) (name TBD) plus documentation of the union-pattern idiom. Small.

9. Explicitly not core's job

To keep the boundary honest, these stay in the package: contact search heuristics and segmentation/clipping algorithms (2D segments, 3D polygon clipping, integration-cell triangulation), averaged/continuous normal fields, gap functions and their linearizations (incl. the full large-deformation directional derivatives of the mortar operators), friction laws and return mappings, active-set / semismooth-Newton / Uzawa strategies, dual-basis construction, contact-specific stabilization (Nitsche parameters), and solver orchestration. All of it is ordinary Julia code over the primitives above; none of it needs to live in Ferrite.

10. Suggested roadmap

Priority	Item	Size	Unlocks
1	Gap B: custom-quadrature reinit! for FacetValues (+ documented weight conventions)	small PR	mortar, NTS, Nitsche, XFEM, cut-cell
2	Gap C: public facet parametrization contract + facet geometry helpers	small PR	version-stable surface packages
3	Gap D: facetset closest-point query (unlock find_local_coordinate, KDTree reuse, coordinate-vector input)	medium PR	3D + large-deformation search
4	Gap F: apply ConstraintHandler condensation to user blocks	small-medium PR	contact + MPC/periodic composition
5	Gap G: sparsity helpers for pair couplings	small PR	ergonomics
6	Gap A: facet-attached (trace) fields	design + large PR	multiplier methods, dual mortar condensation, mixed interface methods
7	Gap E: custom-points InterfaceValues / bring-your-own-pairs cache	follow-up to 1	idiomatic two-sided API

A penalty- or Nitsche-based contact package is fully enabled by items 1–2 (plus 3 for 3D and large deformation). A multiplier-based one — the method of choice for robust contact — additionally requires item 6. Everything else is quality of life.