FEniCS / FEniCS/basix

create_custom_element fails LU factorization for specific anisotropic tensor-product entity layouts, independent of any rank deficiency

Open
#1,067 2 comments 0 reactions 0 assignees View on GitHub
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

Open the contributing 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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.