🐛[BUG]: Float16 cell areas, cell normals, and Butterfly subdivision overflow, though the results fit in float16 #2248

Version

main (2.3.0a0, ce0eaf2)

On which installation method(s) does this occur?

Source

Describe the issue

While making mesh transforms, Loop subdivision, and sampling work on float16 meshes in #2230, I found more places where float16 intermediates overflow, though the results fit in float16. #2230 does not fix them.

• cell_areas and cell_normals build cross products and determinants in float16, whose maximum is 65504. A right triangle with legs of 300 has an area of 45000, which fits in float16. Its area is inf in 2D, and in 3D both its area and its normal are nan.
• Butterfly subdivision adds positions in float16, so an open mesh near x = 33000 gets infinite edge points. Linear subdivision is not affected.

With #2230, transform updates cached areas and normals in float32. So scale(300.0) on a float16 triangle with legs of 1 gives a finite area if it was cached before the call, and inf (2D) or nan (3D) if not.

Expected: finite areas, normals, and subdivided points whenever they fit in float16. One fix is to compute in at least float32 and cast the result back once, as #2230 does for Loop subdivision.

Minimum reproducible example

import torch
from physicsnemo.mesh import Mesh

points = torch.tensor([[0.0, 0.0, 0.0], [300.0, 0.0, 0.0], [0.0, 300.0, 0.0]])
cells = torch.tensor([[0, 1, 2]])

surface = Mesh(points=points.half(), cells=cells)
print(surface.cell_areas, surface.cell_normals)

planar = Mesh(points=points[:, :2].half(), cells=cells)
print(planar.cell_areas)

grid = torch.tensor([[0.0, 0.0], [1.0, 0.0], [0.5, 1.0], [1.5, 1.0]])
open_mesh = Mesh(
    points=(1000 * grid + torch.tensor([33000.0, 0.0])).half(),
    cells=torch.tensor([[0, 1, 2], [1, 3, 2]]),
)
print(open_mesh.subdivide(1, filter="butterfly").points.isinf().any())

Relevant log output

tensor([nan], dtype=torch.float16) tensor([[nan, nan, nan]], dtype=torch.float16)
tensor([inf], dtype=torch.float16)
tensor(True)

Expected:

tensor([44992.], dtype=torch.float16) tensor([[0., 0., 1.]], dtype=torch.float16)
tensor([44992.], dtype=torch.float16)
tensor(False)

Environment details

PyTorch 2.13.0+cu130, CPU