Gamma_P1 projected normals cause tensor dimension mismatch in solvers

[lmoresi]

Summary

PR #119 switches the default Nitsche surface normals from Gamma_N to Gamma_P1 (projected P1 normals). This causes ValueError: cannot contract or diagonalize between axes of different dimension in sympy tensor operations during solver assembly.

Failing tests (CI on PR #119)

• test_0504_projections.py — scalar projection, gradient recovery
• test_0720_coordinate_units_gradients.py — all gradient tests
• test_1050_VEstokesCart.py — all VE Stokes box mesh tests

All fail with the same error from sympy.tensor.array.arrayop.py:105.

Root cause

The Gamma_P1 components likely have a different sympy array shape/rank than Gamma_N, so when they're substituted into constitutive model tensor expressions (which expect the Gamma_N shape), the contraction dimensions don't match.

Context

PR #119 bundles several independent features:

1. Submesh infrastructure (DMPlexFilter wrapper, region DS)
2. Region labels for internal boundaries
3. mesh.Gamma_P1 projected normals
4. Default normals switch from Gamma_N to Gamma_P1

Items 1-2 are independent and working. Items 3-4 need the tensor shape issue resolved before merging.

Suggestion

Either fix the Gamma_P1 shape to match Gamma_N's tensor structure, or split PR #119 so the working submesh/region infrastructure can merge independently.

Underworld development team with AI support from Claude Code

[lmoresi]

Fixed in PR #119 — the issue was _check_expression_meshes accessing cm.flux (triggering tensor contraction) rather than checking parameter expressions directly.