create_custom_element fails LU factorization for specific anisotropic tensor-product entity layouts, independent of any rank deficiency
- Dominant language
- C++
- Stars
- 145
- Forks
- 47
- Avg merge
- 3h 5m
- Merged PRs (30d)
- 8
Description
## Summary
Building custom hexahedron elements as 3-fold anisotropic tensor products of
1D elements (a C^1 cubic Hermite factor, a C^0 quadratic Lagrange factor, and
a discontinuous linear Lagrange factor — degrees {3, 2, 1} respectively),
`basix.create_custom_element` raises `RuntimeError: LU decomposition failed: N`
for a specific subset of axis orderings of the same three factors, while other
orderings of the *identical* three factors succeed.
Example: for factors (S1, S0, S-1) placed along (x, y, z):
- `(x,y,z) = (S1, S-1, S-1)` → succeeds
- `(x,y,z) = (S-1, S1, S-1)` → succeeds
- `(x,y,z) = (S-1, S-1, S1)` → `LU decomposition failed: 5`
All three orderings define mathematically equivalent (up to relabeling)
anisotropic polynomial spaces of the same dimension.
## What's been ruled out
I built the interpolation matrix (dof functionals applied to the raw
tensor-product basis) by hand in numpy for both passing and failing
orderings, and confirmed:
1. **No rank deficiency.** The hand-built interpolation matrix has full rank
for every failing case, with singular values and condition number
matching the passing cases to machine precision.
2. **No entity collisions.** Entity assignment via `basix.cell.topology`
(matching each combined dof group to a unique face/edge/vertex by its
fixed coordinates) is 1:1 for every case — verified by exhaustively
checking for duplicate entity indices across all dof groups.
3. **`wcoeffs` itself is not the problem.** I pre-inverted the interpolation
matrix by hand (i.e., computed the nodal basis explicitly, expressed in
basix's own orthonormal Legendre polyset coordinates, sanity-checked to
`1e-8` against the identity relation `D @ D_inv ≈ I`) and passed that
already-solved `wcoeffs` into `create_custom_element`. It **still**
produces the identical `LU decomposition failed: N` error code, with the
same `N`, for the same entity layouts. This rules out the wcoeffs →
interpolation-matrix step specifically — whatever internal solve is
failing, it's happening downstream of `wcoeffs`, most likely in
dof-transformation / entity-orientation setup. Possibly related to #2586.
## Environment
- basix 0.11.0
- Cell: hexahedron
- Element family: custom (`basix.create_custom_element`), built from
`ElementFamily.Hermite` degree 3, `ElementFamily.P` (Lagrange, continuous)
degree 2, and `ElementFamily.P` (Lagrange, `discontinuous=True`,
`LagrangeVariant.legendre`) degree 1, combined via quadrature-projected
tensor product.
## Reproducer
Minimal scripts attached: `stage1_1d_elements.py` (the three 1D element
factories), `stage2_tensor_product_builder.py` (the tensor-product custom
element builder), `stage2b_full_coverage_test.py` (sweeps all triples and
reports pass/fail).
Run `python3 stage2b_full_coverage_test.py`.
**Failing triples** (axis order, `1`=Hermite deg 3, `0`=Lagrange deg 2,
`-1`=discontinuous Lagrange deg 1):
`(-1,-1,1)`, `(-1,0,0)`, `(-1,0,1)`, `(-1,1,0)`, `(0,-1,0)`, `(0,-1,1)`
**Passing triples**, including permutations of the same three factors as
several of the failing cases above:
`(0,0,0)`, `(0,0,-1)`, `(0,1,-1)`, `(1,-1,-1)`, `(1,-1,0)`, `(1,0,-1)`,
`(-1,1,-1)`
Contributor guide
Research direction
Run stage2b_full_coverage_test.py using the attached stage1_1d_elements.py and stage2_tensor_product_builder.py scripts to reproduce the ordering-dependent failures. Then inspect basix.create_custom_element and the downstream dof-transformation or entity-orientation setup, including the possible connection to #2586. Done means the failing triples construct successfully while the existing passing cases remain unchanged, with regression coverage for the sweep.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- cpp, python
- Domain
- backend, testing-qa
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Active
- Clarity
- Mostly clear
- Newbie friendliness
- 45/100