Skip to content

JohnsonMercier on manifolds#239

Open
pbrubeck wants to merge 4 commits intomainfrom
pbrubeck/zany-manifold
Open

JohnsonMercier on manifolds#239
pbrubeck wants to merge 4 commits intomainfrom
pbrubeck/zany-manifold

Conversation

@pbrubeck
Copy link
Copy Markdown

@pbrubeck pbrubeck commented Mar 25, 2026

This PR enables FIAT to handle physical cells with higher codimension, and is by no means required to support manifolds in Firedrake (which we already do for some elements, including Johnson-Mercier). This is ONLY to allow testing the basis transformation in FIAT.

@pbrubeck pbrubeck force-pushed the pbrubeck/zany-manifold branch 5 times, most recently from e97f118 to 4874c71 Compare March 25, 2026 14:43
@pbrubeck pbrubeck force-pushed the pbrubeck/zany-manifold branch from 4874c71 to f7a8460 Compare March 25, 2026 14:53
@pbrubeck pbrubeck marked this pull request as ready for review March 25, 2026 21:17
@pbrubeck pbrubeck requested a review from rckirby March 25, 2026 21:22
Comment thread FIAT/expansions.py Outdated
Comment thread FIAT/expansions.py Outdated
Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>
finat.JohnsonMercier,
finat.Morley,
])
def test_piola_manifold(ref_to_phys, element):
Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
def test_piola_manifold(ref_to_phys, element):
def test_manifold(ref_to_phys, element):

Comment thread FIAT/reference_element.py
Comment on lines -1614 to -1615
dim_x = len(xs[0])
dim_y = len(ys[0])
Copy link
Copy Markdown
Author

@pbrubeck pbrubeck Mar 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How does the old code handle different topological dimensions? It seems to work fine going from low dimension to a higher one, but I think you need least-squares in order to go the other way around. This code always builds a square matrix mat, but it might be singular.

The new code builds a possibly rectangular matrix DX if xs defines a simplex with tdim != gdim. If tdim == gdim, DX is square it inverts it with numpy.linalg.solve, and uses numpy.linalg.lstsq otherwise.

Comment thread FIAT/reference_element.py
Comment on lines -1617 to -1618
if len(xs) != len(ys):
raise Exception("")
Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this exception be considered in the new code?

Comment thread FIAT/reference_element.py
if DX.shape[0] == DX.shape[1]:
AT = numpy.linalg.solve(DX, DY)
else:
AT, *_ = numpy.linalg.lstsq(DX, DY)
Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we check for a small residual in the least-squares case?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant