Tabulate high-order expansion derivatives by recurrence by pbrubeck · Pull Request #256 · firedrakeproject/fiat

pbrubeck · 2026-07-10T21:58:24Z

Fixes firedrakeproject/firedrake#4152

Summary

Tabulate arbitrary-order Dubiner expansion derivatives by applying Leibniz to the recurrence factors.
Default simplex expansion tabulation to the recurrence path instead of dense differentiation matrices.
Add order-4 triangle/tetrahedron coverage for orthonormal and bubble variants.

Tests

python -m py_compile FIAT/expansions.py test/FIAT/unit/test_polynomial.py
python -m pytest test/FIAT/unit/test_polynomial.py test/FIAT/unit/test_hierarchical.py::test_hierarchical_sparsity test/FIAT/regression/test_regression.py::test_expansions_jet -q

AI tool used: Codex.

pbrubeck · 2026-07-10T22:27:40Z

Updated verification after the HCT degree-14 follow-up.

What changed:

Stabilized high-order Ck vertex supersmoothness constraints by imposing each constraint block on the current nullspace, instead of burying high-order vertex constraints in one large SVD.
Added high-order HCT regression coverage for degrees 13 and 14.

Verification:

python -m py_compile FIAT/expansions.py FIAT/macro.py test/FIAT/unit/test_hct.py test/FIAT/unit/test_polynomial.py
python -m pytest test/FIAT/unit/test_hct.py test/FIAT/unit/test_macro.py -q
python -m pytest test/FIAT/unit/test_polynomial.py::test_high_order_expansion_derivatives test/FIAT/regression/test_regression.py::test_expansions_jet -q
Firedrake exact issue command for HCT degree 13: passed, V.dim() = 803
Firedrake HCT degree 14: passed, V.dim() = 923
Firedrake HCT degree sweep 13..17: passed with dimensions 803, 923, 1051, 1187, 1331

AI tool used: Codex.

pbrubeck · 2026-07-10T22:37:26Z

    return pts


+@pytest.mark.parametrize("degree, space_dimension", [(13, 127), (14, 144)])


These are rather expensive tests, maybe we don't want them

We could use a lower-degree version of the tests? But are they causing any problems in run-time on the test suite?

pbrubeck · 2026-07-13T19:35:23Z

@@ -95,7 +142,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
    results = [numpy.zeros((num_members,) + (dim,)*k + phi0.shape[1:], dtype=phi0.dtype)


Storage blows up for large k, we need to exploit symmetry

Didn't we already have that problem? Or is this intermediate storage as a result of doing Leibniz instead of hand-coding?

We only ran recurrences for k<=2

pbrubeck · 2026-07-13T19:36:56Z

-                ddphi[inext] += fprev * ddphi[iprev]
+                cur = [result[icur] for result in results]
+                prev = [result[iprev] for result in results]
+                for rank in range(1, order+1):


We should only take the k-th derivative if the polynomial degree >= k

We know the answer for degree < k without computing.

rckirby

After we resolve the comments, I think this makes sense to do. The derivative matrices are kind of clunky and tend to go unstable at high order.

rckirby · 2026-07-13T21:02:55Z

@@ -95,7 +142,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
    results = [numpy.zeros((num_members,) + (dim,)*k + phi0.shape[1:], dtype=phi0.dtype)


Didn't we already have that problem? Or is this intermediate storage as a result of doing Leibniz instead of hand-coding?

rckirby · 2026-07-13T21:04:52Z

+                            for axis in range(rank))
+        product = product.transpose(permutation + tuple(range(rank, product.ndim)))
+        product = product.reshape((dim,) * rank + operand.shape[operand_rank:])
+


What if not rank? Do we handle None being returned correctly?

rckirby · 2026-07-13T21:06:02Z

-                ddphi[inext] += fprev * ddphi[iprev]
+                cur = [result[icur] for result in results]
+                prev = [result[iprev] for result in results]
+                for rank in range(1, order+1):


We know the answer for degree < k without computing.

rckirby · 2026-07-13T21:08:01Z

            return sum(ai != bi for ai, bi in zip(alpha, beta))

-        # Only use dmats if tabulate failed
+        # Use dmats only for derivatives above the configured recurrence order.


Do we have cases (outside of tests) where recurrence_order isn't infinite?

rckirby · 2026-07-13T21:08:38Z

    return pts


+@pytest.mark.parametrize("degree, space_dimension", [(13, 127), (14, 144)])


We could use a lower-degree version of the tests? But are they causing any problems in run-time on the test suite?

pbrubeck added 2 commits July 10, 2026 22:56

Tabulate high-order expansion derivatives by recurrence

cf37617

Stabilize high-order HCT supersmooth constraints

db31c4b

pbrubeck commented Jul 10, 2026

View reviewed changes

pbrubeck mentioned this pull request Jul 11, 2026

DO NOT MERGE Test FIAT 256 firedrakeproject/firedrake#5237

Open

pbrubeck commented Jul 13, 2026

View reviewed changes

rckirby reviewed Jul 13, 2026

View reviewed changes

		return pts


		@pytest.mark.parametrize("degree, space_dimension", [(13, 127), (14, 144)])

		@@ -95,7 +142,7 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
		results = [numpy.zeros((num_members,) + (dim,)*k + phi0.shape[1:], dtype=phi0.dtype)

Uh oh!

Conversation

pbrubeck commented Jul 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Tests

Uh oh!

pbrubeck commented Jul 10, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pbrubeck Jul 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rckirby left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

pbrubeck commented Jul 10, 2026 •

edited

Loading

pbrubeck Jul 13, 2026 •

edited

Loading