Skip to content

Support physical_normals on manifolds #4994

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
pbrubeck wants to merge 3 commits into
base: main
Choose a base branch
from
+22 −14
Draft

Support physical_normals on manifolds #4994

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f090a2a
Support physical_normals on manifolds
pbrubeck Mar 26, 2026
508828b
Remove NaN from gem.Conditional, this triggers a loopy warning
pbrubeck Mar 27, 2026
671b615
DROP BEFORE MERGE
pbrubeck Mar 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .github/workflows/core.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -210,6 +210,8 @@ jobs:
--extra-index-url https://download.pytorch.org/whl/cpu \
"$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]"

: # UNDO ME
pip install -v --no-deps --ignore-installed git+https://github.com/firedrakeproject/fiat.git@pbrubeck/zany-manifold
firedrake-clean
pip list

Expand Down
29 changes: 19 additions & 10 deletions tsfc/fem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,19 +268,28 @@ def physical_tangents(self):
return rts @ jac.T

def physical_normals(self):
domain = extract_unique_domain(self.mt.terminal)
gdim = domain.geometric_dimension
tdim = domain.topological_dimension
cell = self.interface.fiat_cell
sd = cell.get_spatial_dimension()
num_faces = len(cell.get_topology()[sd-1])
if isinstance(cell, UFCSimplex) and sd == 2:
num_faces = len(cell.get_topology()[tdim-1])
if isinstance(cell, UFCSimplex) and tdim == 1:
pts = self.physical_tangents()
return gem.ListTensor([[pts[0, j] for j in range(gdim)] for i in range(num_faces)])
elif isinstance(cell, UFCSimplex) and tdim == 2 and gdim == 2:
pts = self.physical_tangents()
return gem.ListTensor([[pts[i, 1], -1*pts[i, 0]] for i in range(num_faces)])
elif isinstance(cell, UFCSimplex) and sd == 3:
t = ufl.classes.CellEdgeVectors(extract_unique_domain(self.mt.terminal))
edges = cell.get_connectivity()[(sd-1, 1)]
return gem.ListTensor([[pts[i, 1], -pts[i, 0]] for i in range(num_faces)])
elif isinstance(cell, UFCSimplex) and gdim == 3:
t = ufl.classes.CellEdgeVectors(domain)
normalize = lambda x: x / ufl.sqrt(ufl.dot(x, x))
expr = ufl.as_tensor([-2.0*normalize(ufl.cross(t[edges[i][0], :], t[edges[i][1], :]))
for i in range(num_faces)])
return self.translate_point_expression(expr)
if tdim == 2:
R = ufl.cross(t[0, :], t[1, :])
exprs = [normalize(ufl.cross(t[i, :], R)) for i in range(num_faces)]
else:
edges = cell.get_connectivity()[(tdim-1, 1)]
exprs = [-2.0*normalize(ufl.cross(t[edges[i][0], :], t[edges[i][1], :]))
for i in range(num_faces)]
return self.translate_point_expression(ufl.as_tensor(exprs))
else:
raise NotImplementedError("Can't do physical normals on that cell yet")

Expand Down
5 changes: 1 addition & 4 deletions tsfc/kernel_interface/common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,10 +73,7 @@ def cell_orientation(self, domain, restriction):
f = {None: 0, '+': 0, '-': 1}[restriction]
co_int = self._cell_orientations[domain][f]
return gem.Conditional(gem.Comparison("==", co_int, gem.Literal(1)),
gem.Literal(-1),
gem.Conditional(gem.Comparison("==", co_int, gem.Zero()),
gem.Literal(1),
gem.Literal(numpy.nan)))
Copy link
Copy Markdown
Contributor Author

@pbrubeck pbrubeck Mar 27, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

This code generates a gem expression that takes self._cell_orientations on a given side of a facet and returns -1 if co == 1, +1 if co == 0, and NaN otherwise.

I am removing the NaN branch because it generates a loopy warning at compile time:

/home/pbrubeck/firedrake-venv/lib/python3.10/site-packages/pytools/persistent_dict.py:273: NanKeyWarning: Encountered a NaN while hashing. Since NaNs compare unequal to themselves, the resulting key can not be retrieved from a PersistentDict and will lead to a collision error on retrieval.
  warn("Encountered a NaN while hashing. Since NaNs compare unequal "
/home/pbrubeck/firedrake-venv/lib/python3.10/site-packages/loopy/target/c/codegen/expression.py:472: UserWarning: Encountered 'bare' floating point NaN value. Since NaN != NaN, this leads to problems with cache retrieval. Consider using `pymbolic.primitives.NaN` instead of `math.nan`. The generated code will be equivalent with the added benefit of sound pickling/unpickling of kernel objects.
  warn("Encountered 'bare' floating point NaN value. Since NaN != NaN,"

Copy link
Copy Markdown
Contributor 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.

Maybe the right thing is to

Consider using pymbolic.primitives.NaN instead of math.nan

@connorjward what do you think?

Copy link
Copy Markdown
Contributor

@connorjward connorjward Mar 27, 2026

Choose a reason for hiding this comment

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

I would suggest trying it out and seeing if it works. Seems reasonable to try.

Copy link
Copy Markdown
Contributor 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.

pymbolic.primitives.NaN is a type, and we cannot just create a numpy.ndarray with dtype=object. Ideally we would intercept NaN in the gem -> loopy translation.

firedrake/tsfc/loopy.py

Lines 239 to 240 in 43a5050

if isinstance(temp, gem.Constant):
data.append(lp.TemporaryVariable(name, shape=temp.shape, dtype=dtype, initializer=temp.array, address_space=lp.AddressSpace.LOCAL, read_only=True))

Copy link
Copy Markdown
Contributor

@connorjward connorjward Mar 27, 2026

Choose a reason for hiding this comment

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

Add a gem.NaN type?

gem.Literal(-1), gem.Literal(1))

def cell_size(self, domain, restriction):
if not hasattr(self, "_cell_sizes"):
Expand Down
Loading