From 9b1829d6ea35e2a73609c3edfc67d4c5f17310ad Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 16 Apr 2026 15:20:42 +0100 Subject: [PATCH 1/7] Fix Netgen permutation --- firedrake/mesh.py | 14 ++++++++------ firedrake/netgen.py | 16 +++++++++++++--- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 394f0cb592..65a47fdabf 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3020,10 +3020,13 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent nodes = fiat_element.dual_basis() ref_pts = [] - for node in nodes: - # Assert singleton point for each node. - pt, = node.get_point_dict().keys() - ref_pts.append(pt) + entity_ids = fiat_element.entity_dofs() + for dim in sorted(entity_ids): + for entity in sorted(entity_ids[dim]): + for i in entity_ids[dim][entity]: + # Assert singleton point for each node. + pt, = nodes[i].get_point_dict().keys() + ref_pts.append(pt) reference_points = np.array(ref_pts) # Construct numpy arrays for physical domain data @@ -3033,8 +3036,8 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): curved_points = np.zeros( (ng_dimension, reference_points.shape[0], self.geometric_dimension) ) + self.netgen_mesh.Curve(1) self.netgen_mesh.CalcElementMapping(reference_points, physical_points) - # NOTE: This will segfault for MeshHierarchy on a netgen CSG geometry self.netgen_mesh.Curve(order) self.netgen_mesh.CalcElementMapping(reference_points, curved_points) curved = ng_element.NumPy()["curved"] @@ -3059,7 +3062,6 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): permutation = find_permutation( own_physical_points, new_coordinates.dat.data_ro_with_halos[broken_indices].real, - tol=permutation_tol, ) self.comm.Barrier() # Apply the permutation to each cell in turn diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 57b442b83a..178dbb6a6d 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -80,8 +80,7 @@ def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, @PETSc.Log.EventDecorator() -def find_permutation(points_a: np.ndarray, points_b: np.ndarray, - tol: float = 1e-5): +def find_permutation(points_a: np.ndarray, points_b: np.ndarray): """ Find all permutations between a list of two sets of points. Given two numpy arrays of shape (ncells, npoints, dim) containing @@ -95,7 +94,18 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray, if points_a.shape != points_b.shape: raise ValueError("`points_a` and `points_b` must have the same shape.") - p = [np.where(cdist(a, b).T < tol)[1] for a, b in zip(points_a, points_b)] + dim = points_a.shape[-1] + vids = list(range(dim+1)) + + bs = points_a[:, vids[0], :] + As = points_a[:, vids[1:], :] + As -= bs[:, None, :] + Ainvs = np.linalg.inv(As) + + ref_points_a = np.matmul(points_a - bs[:, None, :], Ainvs) + ref_points_b = np.matmul(points_b - bs[:, None, :], Ainvs) + + p = [np.argmin(cdist(a, b), axis=0) for a, b in zip(ref_points_a, ref_points_b)] if len(p) == 0: return p From 79f23e0778fc3cf0432dbe149bd2c5f25aed546e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 Jul 2026 11:03:49 +0100 Subject: [PATCH 2/7] Deprecate permutation_tol, add test --- firedrake/mesh.py | 12 +++++++++--- firedrake/mg/netgen.py | 2 +- tests/firedrake/regression/test_netgen.py | 24 +++++++++++++++++++++-- 3 files changed, 32 insertions(+), 6 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 65a47fdabf..c2733f48b5 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2979,14 +2979,14 @@ def refine_marked_elements(self, mark, netgen_flags=None): netgen_flags=netgen_flags) @PETSc.Log.EventDecorator() - def curve_field(self, order, permutation_tol=1e-8, cg_field=None): + def curve_field(self, order, permutation_tol=None, cg_field=None): '''Return a function containing the curved coordinates of the mesh. This method requires that the mesh has been constructed from a netgen mesh. :arg order: the order of the curved mesh. - :arg permutation_tol: tolerance used to construct the permutation of the reference element. + :arg permutation_tol: ignored. :arg cg_field: return a CG function field representing the mesh, as opposed to a DG field. Defaults to the continuity of the coordinates of the original mesh. @@ -2998,6 +2998,12 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): if not hasattr(self, "netgen_mesh"): raise ValueError("Cannot curve a mesh that has not been generated by netgen.") + if permutation_tol is not None: + warnings.warn( + "permutation_tol is no longer required to obtain the curved coordinates. " + "This kwarg will be removed in a future release.", + FutureWarning, + ) if cg_field is None: cg_field = not self.coordinates.function_space().finat_element.is_dg() @@ -3427,7 +3433,7 @@ def Mesh(meshfile, **kwargs): # Curve the mesh, if requested degree = netgen_flags.get("degree", 1) if degree != 1: - permutation_tol = netgen_flags.get("permutation_tol", 1e-8) + permutation_tol = netgen_flags.get("permutation_tol", None) cg = netgen_flags.get("cg", None) coordinates = mesh.curve_field( order=degree, diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 83f98bbbef..f6088528ed 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -242,7 +242,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): order = flags.get("degree", 1) if isinstance(order, int): order = [order]*(levs+1) - permutation_tol = flags.get("permutation_tol", 1e-8) + permutation_tol = flags.get("permutation_tol", None) refType = flags.get("refinement_type", "uniform") optMoves = flags.get("optimisation_moves", False) snap = flags.get("snap_to", "geometry") diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 8524b52748..64abf21c31 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -34,14 +34,34 @@ def test_netgen_csg_mesh_high_order(): assert mesh3.coordinates.function_space().ufl_element().degree() == order -def square_geometry(h): +def square_geometry(h, L=np.pi): from netgen.geom2d import SplineGeometry geo = SplineGeometry() - geo.AddRectangle((0, 0), (np.pi, np.pi), bc="rect") + geo.AddRectangle((0, 0), (L, L), bc="rect") ngmesh = geo.GenerateMesh(maxh=h) return ngmesh +def circle_geometry(h, R=1.0): + from netgen.geom2d import SplineGeometry + geo = SplineGeometry() + geo.AddCircle((0, 0), R, bc="circ") + ngmesh = geo.GenerateMesh(maxh=h) + return ngmesh + + +@pytest.mark.parametrize("scale", (1E-5, 1E5)) +def test_high_order(scale): + # Test scale independence of high-order geometry + expected = np.pi * scale * scale + ngmesh = circle_geometry(h=scale/4, R=scale) + + degree = 3 + msh = Mesh(ngmesh, netgen_flags={"degree": degree}) + assert msh.coordinates.function_space().ufl_element().degree() == degree + assert np.isclose(assemble(1*dx(domain=msh)), expected) + + def poisson(h, degree=2): import netgen comm = COMM_WORLD From 44e9db0d1e6d21948af4372724f565325c3208ac Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 Jul 2026 11:07:08 +0100 Subject: [PATCH 3/7] comment --- firedrake/netgen.py | 1 + 1 file changed, 1 insertion(+) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 178dbb6a6d..60073e85d6 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -94,6 +94,7 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray): if points_a.shape != points_b.shape: raise ValueError("`points_a` and `points_b` must have the same shape.") + # Match reference points instead of physical points to ensure scale invariance dim = points_a.shape[-1] vids = list(range(dim+1)) From 02b85c78e1ee16a73604b03895c01eddf0d6e346 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 Jul 2026 13:15:49 +0100 Subject: [PATCH 4/7] Apply suggestions from code review Co-authored-by: Pablo Brubeck --- firedrake/netgen.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 60073e85d6..2c2a66c603 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -102,7 +102,7 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray): As = points_a[:, vids[1:], :] As -= bs[:, None, :] Ainvs = np.linalg.inv(As) - + # x_phys = A * x_ref + b <==> x_ref = inv(A) * (x_phys - b) ref_points_a = np.matmul(points_a - bs[:, None, :], Ainvs) ref_points_b = np.matmul(points_b - bs[:, None, :], Ainvs) From 7a2cc77e70e442c505cccf385c1da148a714c1c6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 Jul 2026 13:18:55 +0100 Subject: [PATCH 5/7] Apply suggestion from @pbrubeck --- firedrake/netgen.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 2c2a66c603..9d19a861c7 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -97,7 +97,7 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray): # Match reference points instead of physical points to ensure scale invariance dim = points_a.shape[-1] vids = list(range(dim+1)) - + # Infer the affine mapping (A, b) from the image of the vertices bs = points_a[:, vids[0], :] As = points_a[:, vids[1:], :] As -= bs[:, None, :] From 93642576a847d78a5cdc4595a9c81c4d7aede785 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 Jul 2026 13:44:11 +0100 Subject: [PATCH 6/7] Apply suggestion from @pbrubeck --- firedrake/netgen.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 9d19a861c7..e0ac8ebd52 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -98,13 +98,13 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray): dim = points_a.shape[-1] vids = list(range(dim+1)) # Infer the affine mapping (A, b) from the image of the vertices - bs = points_a[:, vids[0], :] + bs = points_a[:, vids[:1], :] As = points_a[:, vids[1:], :] - As -= bs[:, None, :] + As -= bs Ainvs = np.linalg.inv(As) # x_phys = A * x_ref + b <==> x_ref = inv(A) * (x_phys - b) - ref_points_a = np.matmul(points_a - bs[:, None, :], Ainvs) - ref_points_b = np.matmul(points_b - bs[:, None, :], Ainvs) + ref_points_a = np.matmul(points_a - bs, Ainvs) + ref_points_b = np.matmul(points_b - bs, Ainvs) p = [np.argmin(cdist(a, b), axis=0) for a, b in zip(ref_points_a, ref_points_b)] From f11a016e76faddc0a31c716118556e9a8a63072e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 Jul 2026 14:36:35 +0100 Subject: [PATCH 7/7] Apply suggestions from code review Co-authored-by: Pablo Brubeck --- firedrake/netgen.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index e0ac8ebd52..e44af74140 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -97,12 +97,13 @@ def find_permutation(points_a: np.ndarray, points_b: np.ndarray): # Match reference points instead of physical points to ensure scale invariance dim = points_a.shape[-1] vids = list(range(dim+1)) - # Infer the affine mapping (A, b) from the image of the vertices + # Infer the affine mapping (A, b) from the image of the vertices (first dim+1 dofs) bs = points_a[:, vids[:1], :] As = points_a[:, vids[1:], :] As -= bs Ainvs = np.linalg.inv(As) # x_phys = A * x_ref + b <==> x_ref = inv(A) * (x_phys - b) + # Multiply inv(A) from the right, since the data is row-major ref_points_a = np.matmul(points_a - bs, Ainvs) ref_points_b = np.matmul(points_b - bs, Ainvs)