From 8c6f918e43ee3a97bccaa3cac4e533aa1d86c3f6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 27 Feb 2026 09:39:13 +0000 Subject: [PATCH 01/25] Netgen: curved meshes from Mesh constructor and refine_marked_elements --- firedrake/mesh.py | 33 +++++++++++++++++++++++++++------ firedrake/mg/netgen.py | 26 +++++++++++--------------- 2 files changed, 38 insertions(+), 21 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index a05d93c3ed..2c4177db2a 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2994,9 +2994,10 @@ def refine_marked_elements(self, mark, netgen_flags=None): utils.check_netgen_installed() if netgen_flags is None: - netgen_flags = {} + netgen_flags = self.netgen_flags DistParams = self._distribution_parameters - els = {2: self.netgen_mesh.Elements2D, 3: self.netgen_mesh.Elements3D} + netgen_mesh = self.netgen_mesh.Copy() + els = {2: netgen_mesh.Elements2D, 3: netgen_mesh.Elements3D} dim = self.geometric_dimension refine_faces = netgen_flags.get("refine_faces", False) if dim in [2, 3]: @@ -3019,12 +3020,14 @@ def refine_marked_elements(self, mark, netgen_flags=None): else: el.refine = False if not refine_faces and dim == 3: - self.netgen_mesh.Elements2D().NumPy()["refine"] = 0 - self.netgen_mesh.Refine(adaptive=True) + netgen_mesh.Elements2D().NumPy()["refine"] = 0 + netgen_mesh.Refine(adaptive=True) mark = mark-np.ones(mark.shape) - return fd.Mesh(self.netgen_mesh, distribution_parameters=DistParams, comm=self.comm) + return fd.Mesh(netgen_mesh, distribution_parameters=DistParams, + netgen_flags=netgen_flags, comm=self.comm) return fd.Mesh(netgen.libngpy._meshing.Mesh(dim), - distribution_parameters=DistParams, comm=self.comm) + distribution_parameters=DistParams, + netgen_flags=netgen_flags, comm=self.comm) else: raise NotImplementedError("No implementation for dimension other than 2 and 3.") @@ -3431,6 +3434,24 @@ def Mesh(meshfile, **kwargs): if from_netgen: mesh.netgen_mesh = netgen_firedrake_mesh.meshMap.ngMesh + mesh.netgen_flags = netgen_flags + + # Curve the mesh, if requested + degree = netgen_flags.get("degree", 1) + if degree != 1: + permutation_tol = netgen_flags.get("permutation_tol", 1e-8) + cg = netgen_flags.get("cg", False) + ho_field = mesh.curve_field( + order=degree, + permutation_tol=permutation_tol, + cg_field=cg + ) + temp = Mesh(ho_field, distribution_parameters=mesh._distribution_parameters, + comm=mesh.comm) + temp.netgen_mesh = mesh.netgen_mesh + temp.netgen_flags = netgen_flags + temp._tolerance = mesh.tolerance + mesh = temp mesh.submesh_parent = submesh_parent mesh._tolerance = tolerance diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 82acab25aa..ae3edb7a6b 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -246,8 +246,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): logger.info(f"\tOrder of the hierarchy: {order}") if isinstance(order, int): order = [order]*(levs+1) - permutation_tol = flags.get("permutation_tol", 1e-8) - location_tol = flags.get("location_tol", 1e-8) refType = flags.get("refinement_type", "uniform") logger.info(f"\tRefinement type: {refType}") optMoves = flags.get("optimisation_moves", False) @@ -266,15 +264,11 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): parameters.update(mesh._distribution_parameters) parameters["partition"] = False # Curve the mesh - if order[0] > 1: - ho_field = mesh.curve_field( - order=order[0], - permutation_tol=permutation_tol, - cg_field=cg - ) - temp = fd.Mesh(ho_field, distribution_parameters=parameters, comm=comm) - temp.netgen_mesh = mesh.netgen_mesh - temp._tolerance = mesh.tolerance + if order[0] != mesh.coordinates.function_space().ufl_element().degree(): + temp_flags = dict(flags) + temp_flags['degree'] = order[0] + temp = fd.Mesh(mesh.netgen_mesh, distribution_parameters=parameters, + netgen_flags=temp_flags, comm=comm) mesh = temp # Make a plex (cdm) without overlap. dm_cell_type, = mesh.dm_cell_types @@ -313,17 +307,19 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): raise ValueError("Only 2D and 3D meshes can be optimised.") mesh.netgen_mesh = ngmesh # Curve the mesh - if order[l+1] > 1: + if order[l+1] != mesh.coordinates.function_space().ufl_element().degree(): logger.info("\t\t\tCurving the mesh ...") tic = time.time() if snap == "geometry": + temp_flags = dict(flags) + temp_flags['degree'] = order[l+1] mesh = fd.Mesh( - mesh.curve_field(order=order[l+1], - location_tol=location_tol, - permutation_tol=permutation_tol), + mesh.netgen_mesh, distribution_parameters=parameters, + netgen_flags=temp_flags, comm=comm) elif snap == "coarse": + ho_field = meshes[0].coordinates mesh = snapToCoarse(ho_field, mesh, order[l+1], snap_smoothing, cg) toc = time.time() logger.info(f"\t\t\tMeshed curved. Time taken: {toc-tic}") From db3396fd5282c56ee3b2f8312b11457b6b425c5c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 27 Feb 2026 10:07:02 +0000 Subject: [PATCH 02/25] tests --- firedrake/mesh.py | 10 ++++++++-- tests/firedrake/multigrid/test_netgen_gmg.py | 2 ++ tests/firedrake/regression/test_netgen.py | 21 ++++++++++++++++++++ 3 files changed, 31 insertions(+), 2 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 2c4177db2a..a80f69e09c 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3446,11 +3446,17 @@ def Mesh(meshfile, **kwargs): permutation_tol=permutation_tol, cg_field=cg ) - temp = Mesh(ho_field, distribution_parameters=mesh._distribution_parameters, + # Do not redistribute the mesh + distribution_parameters_noop = {"partition": False, + "overlap_type": (DistributedMeshOverlapType.NONE, 0)} + reorder_noop = None + temp = Mesh(ho_field, + reorder=reorder_noop, + distribution_parameters=distribution_parameters_noop, comm=mesh.comm) temp.netgen_mesh = mesh.netgen_mesh temp.netgen_flags = netgen_flags - temp._tolerance = mesh.tolerance + temp._distribution_parameters = mesh._distribution_parameters mesh = temp mesh.submesh_parent = submesh_parent diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 2b4dcc6602..1093dec38e 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -22,6 +22,8 @@ def test_netgen_mg_circle(): mesh = Mesh(ngmesh) nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3}) mesh = nh[-1] + for m in nh: + assert m.coordinates.function_space().ufl_element().degree() == 3 V = FunctionSpace(mesh, "CG", 3) diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index dca547a9f0..12a29e2202 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -3,6 +3,27 @@ import pytest +@pytest.mark.skipnetgen +def test_create_netgen_mesh_high_order(): + from netgen.geom2d import Circle, CSG2d + geo = CSG2d() + + circle = Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle") + geo.Add(circle) + + ngmesh = geo.GenerateMesh(maxh=0.75) + + # Test that setting the degree in netgen_flags produces a high-order mesh + mesh = Mesh(ngmesh, netgen_flags={"degree": 3}) + assert mesh.coordinates.function_space().ufl_element().degree() == 3 + + # Test that refining a high-order mesh gives a high-order mesh + DG0 = FunctionSpace(mesh, "DG", 0) + markers = Function(DG0).assign(1) + mesh2 = mesh.refine_marked_elements(markers) + assert mesh2.coordinates.function_space().ufl_element().degree() == 3 + + def square_geometry(h): from netgen.geom2d import SplineGeometry geo = SplineGeometry() From 86a1af2226f72013cbcd58fa898a4693241a5c0f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 11 Mar 2026 12:26:48 +0000 Subject: [PATCH 03/25] update --- firedrake/mesh.py | 43 +++++++++++++++++++++---------------------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 90f561570f..3de0d6989f 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -543,6 +543,7 @@ def __init__(self, topology_dm, name, reorder, sfXB, perm_is, distribution_name, self.sfXB = sfXB r"The PETSc SF that pushes the global point number slab [0, NX) to input (naive) plex." self.submesh_parent = submesh_parent + self.sfBC_orig = None # User comm self.user_comm = comm dmcommon.label_facets(self.topology_dm) @@ -1145,6 +1146,7 @@ def _distribute(self): sfBC = plex.distribute(overlap=0) plex.setName(original_name) self.sfBC = sfBC + self.sfBC_orig = sfBC # plex carries a new dm after distribute, which # does not inherit partitioner from the old dm. # It probably makes sense as chaco does not work @@ -2928,7 +2930,6 @@ def refine_marked_elements(self, mark, netgen_flags=None): if netgen_flags is None: netgen_flags = self.netgen_flags - DistParams = self._distribution_parameters netgen_mesh = self.netgen_mesh.Copy() els = {2: netgen_mesh.Elements2D, 3: netgen_mesh.Elements3D} dim = self.geometric_dimension @@ -2937,30 +2938,28 @@ def refine_marked_elements(self, mark, netgen_flags=None): with mark.dat.vec as marked: marked0 = marked getIdx = self._cell_numbering.getOffset - if self.sfBC is not None: - sfBCInv = self.sfBC.createInverse() + sf = self.sfBC_orig + if sf is not None: + sfBCInv = sf.createInverse() getIdx = lambda x: x _, marked0 = self.topology_dm.distributeField(sfBCInv, self._cell_numbering, marked) - if self.comm.Get_rank() == 0: + if self.comm.rank == 0: mark = marked0.getArray() - max_refs = np.max(mark) - for _ in range(int(max_refs)): + max_refs = int(np.max(mark)) + for _ in range(max_refs): for i, el in enumerate(els[dim]()): - if mark[getIdx(i)] > 0: - el.refine = True - else: - el.refine = False + el.refine = mark[getIdx(i)] > 0 if not refine_faces and dim == 3: netgen_mesh.Elements2D().NumPy()["refine"] = 0 netgen_mesh.Refine(adaptive=True) - mark = mark-np.ones(mark.shape) - return fd.Mesh(netgen_mesh, distribution_parameters=DistParams, - netgen_flags=netgen_flags, comm=self.comm) - return fd.Mesh(netgen.libngpy._meshing.Mesh(dim), - distribution_parameters=DistParams, - netgen_flags=netgen_flags, comm=self.comm) + mark -= 1 + + return fd.Mesh(netgen_mesh, + reorder=self._did_reordering, + distribution_parameters=self._distribution_parameters, + comm=self.comm) else: raise NotImplementedError("No implementation for dimension other than 2 and 3.") @@ -3374,22 +3373,22 @@ def Mesh(meshfile, **kwargs): if degree != 1: permutation_tol = netgen_flags.get("permutation_tol", 1e-8) cg = netgen_flags.get("cg", False) - ho_field = mesh.curve_field( + coordinates = mesh.curve_field( order=degree, permutation_tol=permutation_tol, - cg_field=cg + cg_field=cg, ) # Do not redistribute the mesh - distribution_parameters_noop = {"partition": False, - "overlap_type": (DistributedMeshOverlapType.NONE, 0)} reorder_noop = None - temp = Mesh(ho_field, + temp = Mesh(coordinates, reorder=reorder_noop, - distribution_parameters=distribution_parameters_noop, + perm_is=mesh._dm_renumbering, + distribution_parameters=DISTRIBUTION_PARAMETERS_NOOP, comm=mesh.comm) temp.netgen_mesh = mesh.netgen_mesh temp.netgen_flags = netgen_flags temp._distribution_parameters = mesh._distribution_parameters + temp._did_reordering = mesh._did_reordering mesh = temp mesh.submesh_parent = submesh_parent From 3488bc4b6cf71a299c6cab8a8be50e678edf19f2 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 11 Mar 2026 12:34:00 +0000 Subject: [PATCH 04/25] clean up --- firedrake/mesh.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 3de0d6989f..11b8983916 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2959,7 +2959,8 @@ def refine_marked_elements(self, mark, netgen_flags=None): return fd.Mesh(netgen_mesh, reorder=self._did_reordering, distribution_parameters=self._distribution_parameters, - comm=self.comm) + comm=self.comm, + netgen_flags=netgen_flags) else: raise NotImplementedError("No implementation for dimension other than 2 and 3.") @@ -3386,7 +3387,7 @@ def Mesh(meshfile, **kwargs): distribution_parameters=DISTRIBUTION_PARAMETERS_NOOP, comm=mesh.comm) temp.netgen_mesh = mesh.netgen_mesh - temp.netgen_flags = netgen_flags + temp.netgen_flags = mesh.netgen_flags temp._distribution_parameters = mesh._distribution_parameters temp._did_reordering = mesh._did_reordering mesh = temp From bb668a5ed495c1cb3a0d18de4ef27871bb705bdc Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 11 Mar 2026 12:53:28 +0000 Subject: [PATCH 05/25] sanitise refine_marked_elements --- firedrake/mesh.py | 63 ++++++++++++++++++++++------------------------- 1 file changed, 30 insertions(+), 33 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 11b8983916..a88ae25745 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2924,46 +2924,43 @@ def refine_marked_elements(self, mark, netgen_flags=None): - refine_faces, which is a boolean specifying if you want to refine faces. """ - import firedrake as fd - utils.check_netgen_installed() if netgen_flags is None: netgen_flags = self.netgen_flags - netgen_mesh = self.netgen_mesh.Copy() - els = {2: netgen_mesh.Elements2D, 3: netgen_mesh.Elements3D} dim = self.geometric_dimension - refine_faces = netgen_flags.get("refine_faces", False) - if dim in [2, 3]: - with mark.dat.vec as marked: - marked0 = marked - getIdx = self._cell_numbering.getOffset - sf = self.sfBC_orig - if sf is not None: - sfBCInv = sf.createInverse() - getIdx = lambda x: x - _, marked0 = self.topology_dm.distributeField(sfBCInv, - self._cell_numbering, - marked) - if self.comm.rank == 0: - mark = marked0.getArray() - max_refs = int(np.max(mark)) - for _ in range(max_refs): - for i, el in enumerate(els[dim]()): - el.refine = mark[getIdx(i)] > 0 - if not refine_faces and dim == 3: - netgen_mesh.Elements2D().NumPy()["refine"] = 0 - netgen_mesh.Refine(adaptive=True) - mark -= 1 - - return fd.Mesh(netgen_mesh, - reorder=self._did_reordering, - distribution_parameters=self._distribution_parameters, - comm=self.comm, - netgen_flags=netgen_flags) - else: + if dim not in {2, 3}: raise NotImplementedError("No implementation for dimension other than 2 and 3.") + with mark.dat.vec as mvec: + if self.sfBC_orig is None: + perm = list(map(self._cell_numbering.getOffset, range(mvec.getSize()))) + mark_np = mvec.getArray()[perm] + else: + sfBCInv = self.sfBC_orig.createInverse() + _, marked0 = self.topology_dm.distributeField(sfBCInv, + self._cell_numbering, + mvec) + mark_np = marked0.getArray() + + netgen_mesh = self.netgen_mesh.Copy() + refine_faces = netgen_flags.get("refine_faces", False) + if self.comm.rank == 0: + max_refs = int(mark_np.max()) + for _ in range(max_refs): + netgen_cells = netgen_mesh.Elements2D() if dim == 2 else netgen_mesh.Elements3D() + netgen_cells.NumPy()["refine"][:mark_np.size] = mark_np > 0 + if not refine_faces and dim == 3: + netgen_mesh.Elements2D().NumPy()["refine"] = 0 + netgen_mesh.Refine(adaptive=True) + mark_np -= 1 + + return Mesh(netgen_mesh, + reorder=self._did_reordering, + distribution_parameters=self._distribution_parameters, + comm=self.comm, + netgen_flags=netgen_flags) + @PETSc.Log.EventDecorator() def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=False): '''Return a function containing the curved coordinates of the mesh. From 74ea35d30e0a66306c840f6866995df355227cc2 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 11 Mar 2026 16:50:15 +0000 Subject: [PATCH 06/25] attempt to fix netgen MeshHierarchy --- firedrake/mesh.py | 27 +++--- firedrake/mg/netgen.py | 86 +++++++++++++------- firedrake/netgen.py | 4 +- tests/firedrake/multigrid/test_netgen_gmg.py | 38 +++++---- 4 files changed, 94 insertions(+), 61 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index a88ae25745..22ff2a9da7 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2948,7 +2948,7 @@ def refine_marked_elements(self, mark, netgen_flags=None): if self.comm.rank == 0: max_refs = int(mark_np.max()) for _ in range(max_refs): - netgen_cells = netgen_mesh.Elements2D() if dim == 2 else netgen_mesh.Elements3D() + netgen_cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() netgen_cells.NumPy()["refine"][:mark_np.size] = mark_np > 0 if not refine_faces and dim == 3: netgen_mesh.Elements2D().NumPy()["refine"] = 0 @@ -2981,12 +2981,13 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F from firedrake.netgen import find_permutation + netgen_mesh = self.netgen_mesh # Check if the mesh is a surface mesh or two dimensional mesh - if len(self.netgen_mesh.Elements3D()) == 0: - ng_element = self.netgen_mesh.Elements2D + if len(netgen_mesh.Elements3D()) == 0: + ng_element = netgen_mesh.Elements2D() else: - ng_element = self.netgen_mesh.Elements3D - ng_dimension = len(ng_element()) + ng_element = netgen_mesh.Elements3D() + ng_dimension = len(ng_element) geom_dim = self.geometric_dimension # Construct the mesh as a Firedrake function @@ -3020,11 +3021,11 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F curved_space_points = np.zeros( (ng_dimension, reference_space_points.shape[0], geom_dim) ) - self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) - # NOTE: This will segfault! - self.netgen_mesh.Curve(order) - self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) - curved = ng_element().NumPy()["curved"] + netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) + # NOTE: This will do nothing for OCC or segfault for CSG! + netgen_mesh.Curve(order) + netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) + curved = ng_element.NumPy()["curved"] # Broadcast a boolean array identifying curved cells curved = self.comm.bcast(curved, root=0) physical_space_points = physical_space_points[curved] @@ -3041,13 +3042,13 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F ) # Broadcast curved cell point data - self.comm.Bcast(physical_space_points, root=0) - self.comm.Bcast(curved_space_points, root=0) + physical_space_points = self.comm.bcast(physical_space_points, root=0) + curved_space_points = self.comm.bcast(curved_space_points, root=0) cell_node_map = new_coordinates.cell_node_map() # Select only the points in curved cells barycentres = np.average(physical_space_points, axis=1) - ng_index = [*map(lambda x: self.locate_cell(x, tolerance=location_tol), barycentres)] + ng_index = list(map(lambda x: self.locate_cell(x, tolerance=location_tol), barycentres)) # Select only the indices of points owned by this rank owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index ae3edb7a6b..17f7733697 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -10,6 +10,7 @@ from petsc4py import PETSc import firedrake as fd +from firedrake.mesh import DISTRIBUTION_PARAMETERS_NOOP from firedrake.cython import mgimpl as impl, dmcommon from firedrake import dmhooks from firedrake.logging import logger @@ -254,22 +255,21 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): logger.info(f"\tSnap to {snap} using {snap_smoothing} smoothing (if snapping to coarse)") cg = flags.get("cg", False) nested = flags.get("nested", snap in ["coarse"]) + permutation_tol = flags.get("permutation_tol", 1e-8) + location_tol = flags.get("location_tol", 1e-8) # Firedrake quantities meshes = [] lgmaps = [] - parameters = {} - if distribution_parameters is not None: - parameters.update(distribution_parameters) + if distribution_parameters is None: + distribution_parameters = mesh._distribution_parameters else: - parameters.update(mesh._distribution_parameters) - parameters["partition"] = False + distribution_parameters.update(mesh._distribution_parameters) # Curve the mesh - if order[0] != mesh.coordinates.function_space().ufl_element().degree(): - temp_flags = dict(flags) - temp_flags['degree'] = order[0] - temp = fd.Mesh(mesh.netgen_mesh, distribution_parameters=parameters, - netgen_flags=temp_flags, comm=comm) - mesh = temp + if mesh.coordinates.function_space().ufl_element().degree() != order[0]: + mesh = curved_mesh(mesh, + order=order[0], + cg_field=cg, + permutation_tol=permutation_tol) # Make a plex (cdm) without overlap. dm_cell_type, = mesh.dm_cell_types tdim = mesh.topology_dm.getDimension() @@ -281,22 +281,15 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): o = impl.create_lgmap(mesh.topology_dm) lgmaps.append((no, o)) mesh.topology_dm.setRefineLevel(0) - meshes.append(mesh) ngmesh = mesh.netgen_mesh + meshes.append(mesh) + for l in range(levs): # Straighten the mesh ngmesh.Curve(1) rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm) cdm = rdm - # Snap the mesh to the Netgen mesh - if snap == "geometry": - snapToNetgenDMPlex(ngmesh, rdm, comm) - # We construct a Firedrake mesh from the DMPlex mesh - no = impl.create_lgmap(rdm) - mesh = fd.Mesh(rdm, dim=meshes[-1].geometric_dimension, reorder=False, - distribution_parameters=parameters, comm=comm) - o = impl.create_lgmap(mesh.topology_dm) - lgmaps.append((no, o)) + if optMoves: # Optimises the mesh, for example smoothing if ngmesh.dim == 2: @@ -305,29 +298,60 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): ngmesh.OptimizeVolumeMesh(MeshingParameters(optimize3d=optMoves)) else: raise ValueError("Only 2D and 3D meshes can be optimised.") + + # Snap the mesh to the Netgen mesh + if snap == "geometry": + snapToNetgenDMPlex(ngmesh, rdm, comm) + + # We construct a Firedrake mesh from the DMPlex mesh + no = impl.create_lgmap(rdm) + mesh = fd.Mesh(rdm, dim=meshes[0].geometric_dimension, + reorder=False, + distribution_parameters=distribution_parameters, + comm=comm) mesh.netgen_mesh = ngmesh + mesh.netgen_flags = flags + + o = impl.create_lgmap(mesh.topology_dm) + lgmaps.append((no, o)) + # Curve the mesh if order[l+1] != mesh.coordinates.function_space().ufl_element().degree(): logger.info("\t\t\tCurving the mesh ...") tic = time.time() if snap == "geometry": - temp_flags = dict(flags) - temp_flags['degree'] = order[l+1] - mesh = fd.Mesh( - mesh.netgen_mesh, - distribution_parameters=parameters, - netgen_flags=temp_flags, - comm=comm) + mesh = curved_mesh(mesh, + order=order[l+1], + location_tol=location_tol, + permutation_tol=permutation_tol, + cg_field=cg, + ) elif snap == "coarse": - ho_field = meshes[0].coordinates - mesh = snapToCoarse(ho_field, mesh, order[l+1], snap_smoothing, cg) + mesh = snapToCoarse(meshes[l].coordinates, mesh, order[l+1], snap_smoothing, cg) + mesh.netgen_mesh = ngmesh + mesh.netgen_flags = flags + toc = time.time() logger.info(f"\t\t\tMeshed curved. Time taken: {toc-tic}") logger.info(f"\t\tLevel {l+1}: with {ngmesh.Coordinates().shape[0]}\ vertices, with order {order[l+1]}, snapping to {snap}\ and optimisation moves {optMoves}.") - mesh.topology_dm.setRefineLevel(1 + l) + mesh.topology_dm.setRefineLevel(l+1) meshes.append(mesh) + # Populate the coarse to fine map coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes, lgmaps) return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, 1, nested=nested) + + +def curved_mesh(mesh, **kwargs): + coordinates = mesh.curve_field(**kwargs) + temp = fd.Mesh(coordinates, reorder=False, + perm_is=mesh._dm_renumbering, + distribution_parameters=DISTRIBUTION_PARAMETERS_NOOP, + comm=mesh.comm) + temp._distribution_parameters = mesh._distribution_parameters + temp._did_reordering = mesh._did_reordering + temp.netgen_mesh = mesh.netgen_mesh + temp._tolerance = mesh.tolerance + return temp diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 1033b2f627..83e6476a20 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -91,10 +91,10 @@ class FiredrakeMesh: :param netgen_flags: The dictionary of flags to be passed to ngsPETSc. :arg comm: the MPI communicator. ''' - def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD): + def __init__(self, mesh, netgen_flags=None, user_comm=fd.COMM_WORLD): self.comm = user_comm # Parsing netgen flags - if not isinstance(netgen_flags, dict): + if netgen_flags is None: netgen_flags = {} split2tets = netgen_flags.get("split_to_tets", False) split = netgen_flags.get("split", False) diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 1093dec38e..a0e4208fba 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -3,12 +3,20 @@ from firedrake import * -def create_netgen_mesh_circle(): - from netgen.geom2d import Circle, CSG2d - geo = CSG2d() - - circle = Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle") - geo.Add(circle) +@pytest.fixture(params=["occ"]) +def circle_geometry(request): + backend = request.param + if backend == "occ": + from netgen.occ import Circle, OCCGeometry + circle = Circle((0, 0), 1.0).Face() + geo = OCCGeometry(circle, dim=2) + elif backend == "csg": + from netgen.geom2d import Circle, CSG2d + circle = Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle") + geo = CSG2d() + geo.Add(circle) + else: + raise ValueError(f"Unexpected backend {backend}") ngmesh = geo.GenerateMesh(maxh=0.75) return ngmesh @@ -17,9 +25,9 @@ def create_netgen_mesh_circle(): @pytest.mark.skip(reason="See https://github.com/firedrakeproject/firedrake/issues/4784") @pytest.mark.skipcomplex @pytest.mark.skipnetgen -def test_netgen_mg_circle(): - ngmesh = create_netgen_mesh_circle() - mesh = Mesh(ngmesh) +def test_netgen_mg_circle(circle_geometry): + mesh = Mesh(circle_geometry) + ngmesh = mesh.netgen_mesh nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3}) mesh = nh[-1] for m in nh: @@ -49,9 +57,9 @@ def test_netgen_mg_circle(): @pytest.mark.skip(reason="See https://github.com/firedrakeproject/firedrake/issues/4784") @pytest.mark.skipcomplex @pytest.mark.skipnetgen -def test_netgen_mg_circle_non_uniform_degree(): - ngmesh = create_netgen_mesh_circle() - mesh = Mesh(ngmesh) +def test_netgen_mg_circle_non_uniform_degree(circle_geometry): + mesh = Mesh(circle_geometry) + ngmesh = mesh.netgen_mesh nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": [1, 2, 3]}) mesh = nh[-1] @@ -80,9 +88,9 @@ def test_netgen_mg_circle_non_uniform_degree(): @pytest.mark.skipcomplex @pytest.mark.skipnetgen @pytest.mark.parallel -def test_netgen_mg_circle_parallel(): - ngmesh = create_netgen_mesh_circle() - mesh = Mesh(ngmesh) +def test_netgen_mg_circle_parallel(circle_geometry): + mesh = Mesh(circle_geometry) + ngmehs = mesh.netgen_mesh nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3}) mesh = nh[-1] From c592b66ba7af0388572bcb95a0bf052ffcce3aac Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 12 Mar 2026 21:17:05 +0000 Subject: [PATCH 07/25] fixes --- firedrake/mg/kernels.py | 2 +- firedrake/mg/netgen.py | 63 ++++---- firedrake/pointquery_utils.py | 7 +- tests/firedrake/multigrid/test_netgen_gmg.py | 143 +++++++------------ tests/firedrake/regression/test_netgen.py | 26 ++++ 5 files changed, 117 insertions(+), 124 deletions(-) diff --git a/firedrake/mg/kernels.py b/firedrake/mg/kernels.py index a0159428c3..f4ce35328f 100644 --- a/firedrake/mg/kernels.py +++ b/firedrake/mg/kernels.py @@ -53,7 +53,7 @@ def to_reference_coordinates(ufl_coordinate_element, parameters=None): "topological_dimension": cell.topological_dimension, "to_reference_coords_newton_step": to_reference_coords_newton_step_body(ufl_coordinate_element, parameters, x0_dtype=ScalarType, dX_dtype="double"), "init_X": init_X(element.cell, parameters), - "max_iteration_count": 1 if is_affine(ufl_coordinate_element) else 16, + "max_iteration_count": 1 if is_affine(ufl_coordinate_element) else 20, "convergence_epsilon": 1e-12, "dX_norm_square": dX_norm_square(cell.topological_dimension), "X_isub_dX": X_isub_dX(cell.topological_dimension), diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 17f7733697..3e60bb1f92 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -260,16 +260,16 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): # Firedrake quantities meshes = [] lgmaps = [] - if distribution_parameters is None: - distribution_parameters = mesh._distribution_parameters - else: - distribution_parameters.update(mesh._distribution_parameters) # Curve the mesh if mesh.coordinates.function_space().ufl_element().degree() != order[0]: - mesh = curved_mesh(mesh, - order=order[0], - cg_field=cg, - permutation_tol=permutation_tol) + coordinates = mesh.curve_field( + order=order[0], + location_tol=location_tol, + permutation_tol=permutation_tol, + cg_field=cg, + ) + mesh = reconstruct_mesh(mesh, coordinates) + # Make a plex (cdm) without overlap. dm_cell_type, = mesh.dm_cell_types tdim = mesh.topology_dm.getDimension() @@ -280,6 +280,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): no = impl.create_lgmap(cdm) o = impl.create_lgmap(mesh.topology_dm) lgmaps.append((no, o)) + mesh.topology_dm.setRefineLevel(0) ngmesh = mesh.netgen_mesh meshes.append(mesh) @@ -304,14 +305,17 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): snapToNetgenDMPlex(ngmesh, rdm, comm) # We construct a Firedrake mesh from the DMPlex mesh - no = impl.create_lgmap(rdm) - mesh = fd.Mesh(rdm, dim=meshes[0].geometric_dimension, - reorder=False, - distribution_parameters=distribution_parameters, - comm=comm) + parameters = {} + if distribution_parameters is not None: + parameters.update(distribution_parameters) + else: + parameters.update(mesh._distribution_parameters) + parameters["partition"] = False + mesh = fd.Mesh(rdm, dim=mesh.geometric_dimension, distribution_parameters=parameters) mesh.netgen_mesh = ngmesh mesh.netgen_flags = flags + no = impl.create_lgmap(rdm) o = impl.create_lgmap(mesh.topology_dm) lgmaps.append((no, o)) @@ -320,16 +324,15 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): logger.info("\t\t\tCurving the mesh ...") tic = time.time() if snap == "geometry": - mesh = curved_mesh(mesh, + coordinates = mesh.curve_field( order=order[l+1], location_tol=location_tol, permutation_tol=permutation_tol, cg_field=cg, ) + mesh = reconstruct_mesh(mesh, coordinates) elif snap == "coarse": - mesh = snapToCoarse(meshes[l].coordinates, mesh, order[l+1], snap_smoothing, cg) - mesh.netgen_mesh = ngmesh - mesh.netgen_flags = flags + mesh = snapToCoarse(meshes[0].coordinates, mesh, order[l+1], snap_smoothing, cg) toc = time.time() logger.info(f"\t\t\tMeshed curved. Time taken: {toc-tic}") @@ -344,14 +347,18 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, 1, nested=nested) -def curved_mesh(mesh, **kwargs): - coordinates = mesh.curve_field(**kwargs) - temp = fd.Mesh(coordinates, reorder=False, - perm_is=mesh._dm_renumbering, - distribution_parameters=DISTRIBUTION_PARAMETERS_NOOP, - comm=mesh.comm) - temp._distribution_parameters = mesh._distribution_parameters - temp._did_reordering = mesh._did_reordering - temp.netgen_mesh = mesh.netgen_mesh - temp._tolerance = mesh.tolerance - return temp +def reconstruct_mesh(mesh, *args, **kwargs): + """Reconstruct a mesh.""" + kwargs.setdefault("dim", mesh.geometric_dimension) + kwargs.setdefault("reorder", False) + kwargs.setdefault("distribution_parameters", DISTRIBUTION_PARAMETERS_NOOP) + kwargs.setdefault("comm", mesh.comm) + kwargs.setdefault("tolerance", mesh.tolerance) + kwargs.setdefault("perm_is", mesh._dm_renumbering) + + tmesh = fd.Mesh(*args, **kwargs) + tmesh._distribution_parameters = mesh._distribution_parameters + tmesh._did_reordering = mesh._did_reordering + tmesh.netgen_mesh = mesh.netgen_mesh + tmesh.netgen_flags = mesh.netgen_flags + return tmesh diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 141eb3a939..4102d0f020 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -128,7 +128,10 @@ def init_X(fiat_cell, parameters): @PETSc.Log.EventDecorator() -def to_reference_coords_newton_step(ufl_coordinate_element, parameters, x0_dtype="double", dX_dtype=ScalarType): +def to_reference_coords_newton_step(ufl_coordinate_element, parameters, + x0_dtype="double", + dX_dtype=ScalarType, + kernel_name="to_reference_coords_newton_step"): # Set up UFL form cell = ufl_coordinate_element.cell domain = ufl.Mesh(ufl_coordinate_element) @@ -191,7 +194,7 @@ def predicate(index): impero_c = impero_utils.compile_gem(assignments, ()) kernel, _ = tsfc.loopy.generate( impero_c, loopy_args, ScalarType, - kernel_name="to_reference_coords_newton_step") + kernel_name=kernel_name) return lp.generate_code_v2(kernel).device_code() diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index a0e4208fba..97b66cf9b7 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -3,113 +3,70 @@ from firedrake import * -@pytest.fixture(params=["occ"]) -def circle_geometry(request): - backend = request.param - if backend == "occ": +@pytest.fixture(params=[2, 3]) +def ngmesh(request): + dim = request.param + if dim == 2: from netgen.occ import Circle, OCCGeometry circle = Circle((0, 0), 1.0).Face() + circle.edges.name = "surface" geo = OCCGeometry(circle, dim=2) - elif backend == "csg": - from netgen.geom2d import Circle, CSG2d - circle = Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle") - geo = CSG2d() - geo.Add(circle) + elif dim == 3: + from netgen.occ import Sphere, OCCGeometry + sphere = Sphere((0, 0, 0), 1.0) + sphere.faces.name = "surface" + geo = OCCGeometry(sphere, dim=3) else: - raise ValueError(f"Unexpected backend {backend}") - + raise ValueError(f"Unexpected dimension {dim}") ngmesh = geo.GenerateMesh(maxh=0.75) return ngmesh -@pytest.mark.skip(reason="See https://github.com/firedrakeproject/firedrake/issues/4784") @pytest.mark.skipcomplex @pytest.mark.skipnetgen -def test_netgen_mg_circle(circle_geometry): - mesh = Mesh(circle_geometry) - ngmesh = mesh.netgen_mesh - nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3}) +# @pytest.mark.parallel([1, 3]) +@pytest.mark.parametrize("netgen_degree", [1, 3, (1, 2, 3)]) +def test_netgen_mg(ngmesh, netgen_degree): + dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)} + mesh = Mesh(ngmesh, distribution_parameters=dparams) + nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": netgen_degree}) mesh = nh[-1] - for m in nh: - assert m.coordinates.function_space().ufl_element().degree() == 3 - - V = FunctionSpace(mesh, "CG", 3) - - u = TrialFunction(V) - v = TestFunction(V) - - a = inner(grad(u), grad(v))*dx - labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]] - bcs = DirichletBC(V, zero(), labels) - x, y = SpatialCoordinate(mesh) - - f = 4+0*x - L = f*v*dx - exact = (1-x**2-y**2) + try: + len(netgen_degree) + except TypeError: + netgen_degree = (netgen_degree,)*len(nh) - u = Function(V) - solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg", - "pc_type": "mg"}) - expect = Function(V).interpolate(exact) - assert (norm(assemble(u - expect)) <= 1e-6) - - -@pytest.mark.skip(reason="See https://github.com/firedrakeproject/firedrake/issues/4784") -@pytest.mark.skipcomplex -@pytest.mark.skipnetgen -def test_netgen_mg_circle_non_uniform_degree(circle_geometry): - mesh = Mesh(circle_geometry) - ngmesh = mesh.netgen_mesh - nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": [1, 2, 3]}) - mesh = nh[-1] + for m, deg in zip(nh, netgen_degree): + assert m.coordinates.function_space().ufl_element().degree() == deg V = FunctionSpace(mesh, "CG", 3) - u = TrialFunction(V) v = TestFunction(V) - a = inner(grad(u), grad(v))*dx - labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]] - bcs = DirichletBC(V, zero(), labels) - x, y = SpatialCoordinate(mesh) - - f = 4+0*x - L = f*v*dx - exact = (1-x**2-y**2) - - u = Function(V) - solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg", - "pc_type": "mg"}) - expect = Function(V).interpolate(exact) - assert (norm(assemble(u - expect)) <= 1e-6) - - -@pytest.mark.skip(reason="See https://github.com/firedrakeproject/firedrake/issues/4784") -@pytest.mark.skipcomplex -@pytest.mark.skipnetgen -@pytest.mark.parallel -def test_netgen_mg_circle_parallel(circle_geometry): - mesh = Mesh(circle_geometry) - ngmehs = mesh.netgen_mesh - nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": 3}) - mesh = nh[-1] - - V = FunctionSpace(mesh, "CG", 3) - - u = TrialFunction(V) - v = TestFunction(V) - - a = inner(grad(u), grad(v))*dx - labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["circle"]] - bcs = DirichletBC(V, zero(), labels) - x, y = SpatialCoordinate(mesh) - - f = 4+0*x - L = f*v*dx - exact = (1-x**2-y**2) - - u = Function(V) - solve(a == L, u, bcs=bcs, solver_parameters={"ksp_type": "cg", - "pc_type": "mg"}) - expect = Function(V).interpolate(exact) - assert norm(assemble(u - expect)) <= 1e-6 + a = inner(grad(u), grad(v)) * dx + labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["surface"]] + + x = SpatialCoordinate(mesh) + uexact = dot(x, x) + bcs = DirichletBC(V, uexact, labels) + L = a(uexact, v) + uh = Function(V) + + rtol = 1E-8 + uerr = uexact - uh + err0 = assemble(a(uerr, uerr)) + solve(a == L, uh, bcs=bcs, solver_parameters={ + "ksp_type": "cg", + "ksp_norm_type": "natural", + "ksp_max_it": 12, + "ksp_rtol": rtol, + "ksp_monitor": None, + "pc_type": "mg", + "mg_levels_pc_type": "python", + "mg_levels_pc_python_type": "firedrake.ASMStarPC", + "mg_levels_pc_star_backend": "tinyasm", + "mg_coarse_pc_type": "lu", + "mg_coarse_pc_factor_mat_solver_type": "mumps", + }) + errf = assemble(a(uerr, uerr)) + assert errf < err0 * (2*rtol) diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 12a29e2202..91d5c8dacd 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -404,3 +404,29 @@ def adapt(mesh, eta): dofs.append(uh.function_space().dim()) mesh = adapt(mesh, eta) assert error_estimators[-1] < 0.06 + + +def test_netgen_manifold(): + from netgen.meshing import MeshingStep + from netgen.occ import Pnt, SplineApproximation, Face, Wire, Axis, OCCGeometry, X, Z, Circle + # Ellipsoid with center (0,R,0) with Y-radius a, Z-radius b + R = 3.0 + a = 1.5 + b = 1.6 + + def Curve(t): + return Pnt(0, R+a*np.cos(t), b*np.sin(t)) + + n = 100 + pnts = [Curve(2*np.pi*t/n) for t in range(n+1)] + + spline = SplineApproximation(pnts) + f = Face(Wire(spline)) + + torus = f.Revolve(Axis((0,0,0), Z), 360) + geo = OCCGeometry(torus, dim=3) + ngmesh = geo.GenerateMesh(maxh=0.5, perfstepsend=MeshingStep.MESHSURFACE) + mesh = Mesh(ngmesh) + + assert mesh.topological_dimension == 2 + assert mesh.geometric_dimension == 3 From 574be15ba2e3042a1da27d46f6bb1137b7e6960c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 12 Mar 2026 22:25:30 +0000 Subject: [PATCH 08/25] fix --- firedrake/mesh.py | 41 ++++++++++---------- firedrake/mg/netgen.py | 13 ++++--- tests/firedrake/multigrid/test_netgen_gmg.py | 1 - tests/firedrake/regression/test_netgen.py | 4 +- 4 files changed, 31 insertions(+), 28 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 22ff2a9da7..848f254859 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3022,10 +3022,11 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F (ng_dimension, reference_space_points.shape[0], geom_dim) ) netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) - # NOTE: This will do nothing for OCC or segfault for CSG! + # NOTE: This will segfault for CSG! netgen_mesh.Curve(order) netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) curved = ng_element.NumPy()["curved"] + # Broadcast a boolean array identifying curved cells curved = self.comm.bcast(curved, root=0) physical_space_points = physical_space_points[curved] @@ -3050,36 +3051,36 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F barycentres = np.average(physical_space_points, axis=1) ng_index = list(map(lambda x: self.locate_cell(x, tolerance=location_tol), barycentres)) - # Select only the indices of points owned by this rank + # Select only the indices of cells owned by this rank owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] + # Get the PyOP2 indices corresponding to the netgen indices + ng_index = [idx for idx, o in zip(ng_index, owned) if o] + pyop2_index = cell_node_map.values[ng_index].flatten() + # Select only the points owned by this rank physical_space_points = physical_space_points[owned] curved_space_points = curved_space_points[owned] barycentres = barycentres[owned] - ng_index = [idx for idx, o in zip(ng_index, owned) if o] - # Get the PyOP2 indices corresponding to the netgen indices - pyop2_index = [] - for ngidx in ng_index: - pyop2_index.extend(cell_node_map.values[ngidx]) - - # Find the correct coordinate permutation for each cell - permutation = find_permutation( - physical_space_points, - new_coordinates.dat.data[pyop2_index].real.reshape( - physical_space_points.shape - ), - tol=permutation_tol - ) + if any(owned): + # Find the correct coordinate permutation for each cell + permutation = find_permutation( + physical_space_points, + new_coordinates.dat.data[pyop2_index].real.reshape( + physical_space_points.shape + ), + tol=permutation_tol + ) - # Apply the permutation to each cell in turn - for ii, p in enumerate(curved_space_points): - curved_space_points[ii] = p[permutation[ii]] + # Apply the permutation to each cell in turn + for ii, p in enumerate(curved_space_points): + curved_space_points[ii] = p[permutation[ii]] + else: + print("barf") # Assign the curved coordinates to the dat new_coordinates.dat.data[pyop2_index] = curved_space_points.reshape(-1, geom_dim) - return new_coordinates diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 3e60bb1f92..9356938de1 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -263,10 +263,10 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): # Curve the mesh if mesh.coordinates.function_space().ufl_element().degree() != order[0]: coordinates = mesh.curve_field( - order=order[0], - location_tol=location_tol, - permutation_tol=permutation_tol, - cg_field=cg, + order=order[0], + location_tol=location_tol, + permutation_tol=permutation_tol, + cg_field=cg, ) mesh = reconstruct_mesh(mesh, coordinates) @@ -311,7 +311,10 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): else: parameters.update(mesh._distribution_parameters) parameters["partition"] = False - mesh = fd.Mesh(rdm, dim=mesh.geometric_dimension, distribution_parameters=parameters) + mesh = fd.Mesh(rdm, dim=mesh.geometric_dimension, + reorder=False, + distribution_parameters=parameters, + tolerance=mesh.tolerance) mesh.netgen_mesh = ngmesh mesh.netgen_flags = flags diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 97b66cf9b7..3e788318b0 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -1,5 +1,4 @@ import pytest - from firedrake import * diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 91d5c8dacd..2d5820d4b6 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -408,7 +408,7 @@ def adapt(mesh, eta): def test_netgen_manifold(): from netgen.meshing import MeshingStep - from netgen.occ import Pnt, SplineApproximation, Face, Wire, Axis, OCCGeometry, X, Z, Circle + from netgen.occ import Pnt, SplineApproximation, Face, Wire, Axis, OCCGeometry, Z # Ellipsoid with center (0,R,0) with Y-radius a, Z-radius b R = 3.0 a = 1.5 @@ -423,7 +423,7 @@ def Curve(t): spline = SplineApproximation(pnts) f = Face(Wire(spline)) - torus = f.Revolve(Axis((0,0,0), Z), 360) + torus = f.Revolve(Axis((0, 0, 0), Z), 360) geo = OCCGeometry(torus, dim=3) ngmesh = geo.GenerateMesh(maxh=0.5, perfstepsend=MeshingStep.MESHSURFACE) mesh = Mesh(ngmesh) From f56305ccec46d0f84515a58b5d13114cf6663e7e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Mar 2026 17:58:58 +0000 Subject: [PATCH 09/25] WIP --- firedrake/mesh.py | 136 ++++++++++------------ firedrake/mg/netgen.py | 3 - firedrake/netgen.py | 66 +++++++++++ tests/firedrake/regression/test_netgen.py | 109 +---------------- 4 files changed, 131 insertions(+), 183 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 848f254859..a58308f8ac 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2926,6 +2926,8 @@ def refine_marked_elements(self, mark, netgen_flags=None): """ utils.check_netgen_installed() + if not hasattr(self, "netgen_mesh"): + raise ValueError("Adaptive refinement requires a netgen mesh.") if netgen_flags is None: netgen_flags = self.netgen_flags dim = self.geometric_dimension @@ -2947,13 +2949,22 @@ def refine_marked_elements(self, mark, netgen_flags=None): refine_faces = netgen_flags.get("refine_faces", False) if self.comm.rank == 0: max_refs = int(mark_np.max()) - for _ in range(max_refs): - netgen_cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() - netgen_cells.NumPy()["refine"][:mark_np.size] = mark_np > 0 + for r in range(max_refs): + cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() + cells.NumPy()["refine"] = (mark_np > 0) if not refine_faces and dim == 3: - netgen_mesh.Elements2D().NumPy()["refine"] = 0 + netgen_mesh.Elements2D().NumPy()["refine"] = False netgen_mesh.Refine(adaptive=True) mark_np -= 1 + if r < max_refs - 1: + parents = netgen_mesh.parentelements if dim == 3 else netgen_mesh.parentsurfaceelements + parents = parents.NumPy().astype(int).flatten() + num_coarse_cells = mark_np.size + num_fine_cells = parents.shape[0] + indices = np.arange(num_fine_cells, dtype=int) + fine_cells = indices > num_coarse_cells + indices[fine_cells] = parents[indices[fine_cells]] + mark_np = mark_np[indices] return Mesh(netgen_mesh, reorder=self._did_reordering, @@ -2962,7 +2973,7 @@ 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, location_tol=1e-1, cg_field=False): + def curve_field(self, order, permutation_tol=1e-8, cg_field=False): '''Return a function containing the curved coordinates of the mesh. This method requires that the mesh has been constructed from a @@ -2970,34 +2981,30 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F :arg order: the order of the curved mesh. :arg permutation_tol: tolerance used to construct the permutation of the reference element. - :arg location_tol: tolerance used to locate the cell a point belongs to. :arg cg_field: return a CG function field representing the mesh, rather than the default DG field. ''' - import firedrake as fd + from firedrake.netgen import find_permutation, netgen_to_plex_numbering, netgen_distribute + from firedrake.functionspace import VectorFunctionSpace, FunctionSpace + from firedrake.function import Function utils.check_netgen_installed() - - from firedrake.netgen import find_permutation - - netgen_mesh = self.netgen_mesh # Check if the mesh is a surface mesh or two dimensional mesh - if len(netgen_mesh.Elements3D()) == 0: - ng_element = netgen_mesh.Elements2D() + if self.topological_dimension == 2: + ng_element = self.netgen_mesh.Elements2D() else: - ng_element = netgen_mesh.Elements3D() + ng_element = self.netgen_mesh.Elements3D() ng_dimension = len(ng_element) - geom_dim = self.geometric_dimension # Construct the mesh as a Firedrake function if cg_field: - firedrake_space = fd.VectorFunctionSpace(self, "CG", order) + coords_space = VectorFunctionSpace(self, "CG", order) else: low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] ufl_element = low_order_element.reconstruct(degree=order) - firedrake_space = fd.VectorFunctionSpace(self, fd.BrokenElement(ufl_element)) - new_coordinates = fd.assemble(fd.interpolate(self.coordinates, firedrake_space)) + coords_space = VectorFunctionSpace(self, finat.ufl.BrokenElement(ufl_element)) + new_coordinates = Function(coords_space).interpolate(self.coordinates) # Compute reference points using fiat fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent @@ -3012,75 +3019,50 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F ref.append(pt) reference_space_points = np.array(ref) - # Curve the mesh on rank 0 only - if self.comm.rank == 0: - # Construct numpy arrays for physical domain data - physical_space_points = np.zeros( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.zeros( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) - # NOTE: This will segfault for CSG! - netgen_mesh.Curve(order) - netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) - curved = ng_element.NumPy()["curved"] - - # Broadcast a boolean array identifying curved cells - curved = self.comm.bcast(curved, root=0) - physical_space_points = physical_space_points[curved] - curved_space_points = curved_space_points[curved] - else: - curved = self.comm.bcast(None, root=0) - # Construct numpy arrays as buffers to receive physical domain data - ncurved = np.sum(curved) - physical_space_points = np.zeros( - (ncurved, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.zeros( - (ncurved, reference_space_points.shape[0], geom_dim) - ) + # Construct numpy arrays for physical domain data + physical_space_points = np.zeros( + (ng_dimension, reference_space_points.shape[0], self.geometric_dimension) + ) + curved_space_points = np.zeros( + (ng_dimension, reference_space_points.shape[0], self.geometric_dimension) + ) + self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) + # NOTE: This will segfault for CSG! + self.netgen_mesh.Curve(order) + self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) + curved = ng_element.NumPy()["curved"] + + # Get numbering + DG0 = FunctionSpace(self, "DG", 0) + rstart, rend = DG0.dof_dset.layout_vec.getOwnershipRange() + num_cells = rend - rstart + _, iperm = netgen_to_plex_numbering(self) + iperm -= rstart + + # Distribute curved cell data + own_curved = netgen_distribute(self, curved) + own_curved = np.flatnonzero(own_curved[:num_cells]) - # Broadcast curved cell point data - physical_space_points = self.comm.bcast(physical_space_points, root=0) - curved_space_points = self.comm.bcast(curved_space_points, root=0) cell_node_map = new_coordinates.cell_node_map() + pyop2_index = cell_node_map.values[iperm[own_curved]] - # Select only the points in curved cells - barycentres = np.average(physical_space_points, axis=1) - ng_index = list(map(lambda x: self.locate_cell(x, tolerance=location_tol), barycentres)) - - # Select only the indices of cells owned by this rank - owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] - - # Get the PyOP2 indices corresponding to the netgen indices - ng_index = [idx for idx, o in zip(ng_index, owned) if o] - pyop2_index = cell_node_map.values[ng_index].flatten() + # Distribute coordinate data + own_curved_points = netgen_distribute(self, curved_space_points)[own_curved] + own_physical_points = netgen_distribute(self, physical_space_points)[own_curved] - # Select only the points owned by this rank - physical_space_points = physical_space_points[owned] - curved_space_points = curved_space_points[owned] - barycentres = barycentres[owned] - - if any(owned): + if any(own_curved): # Find the correct coordinate permutation for each cell permutation = find_permutation( - physical_space_points, - new_coordinates.dat.data[pyop2_index].real.reshape( - physical_space_points.shape - ), - tol=permutation_tol + own_physical_points, + new_coordinates.dat.data[pyop2_index].real, + tol=permutation_tol, ) - # Apply the permutation to each cell in turn - for ii, p in enumerate(curved_space_points): - curved_space_points[ii] = p[permutation[ii]] - else: - print("barf") + for ii, p in enumerate(own_curved_points): + own_curved_points[ii] = p[permutation[ii]] # Assign the curved coordinates to the dat - new_coordinates.dat.data[pyop2_index] = curved_space_points.reshape(-1, geom_dim) + new_coordinates.dat.data[pyop2_index] = own_curved_points return new_coordinates diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 9356938de1..516e94a956 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -256,7 +256,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): cg = flags.get("cg", False) nested = flags.get("nested", snap in ["coarse"]) permutation_tol = flags.get("permutation_tol", 1e-8) - location_tol = flags.get("location_tol", 1e-8) # Firedrake quantities meshes = [] lgmaps = [] @@ -264,7 +263,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): if mesh.coordinates.function_space().ufl_element().degree() != order[0]: coordinates = mesh.curve_field( order=order[0], - location_tol=location_tol, permutation_tol=permutation_tol, cg_field=cg, ) @@ -329,7 +327,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): if snap == "geometry": coordinates = mesh.curve_field( order=order[l+1], - location_tol=location_tol, permutation_tol=permutation_tol, cg_field=cg, ) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 83e6476a20..f908edbe07 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -29,6 +29,72 @@ class comp: Mesh = type(None) +def netgen_distribute(mesh, netgen_data): + from firedrake import FunctionSpace + # Create Netgen to Plex reordering + plex = mesh.topology_dm + sf = mesh.sfBC_orig + perm, iperm = netgen_to_plex_numbering(mesh) + if sf is not None: + netgen_data = np.asarray(netgen_data) + dtype = netgen_data.dtype + dtype = mesh.comm.bcast(dtype, root=0) + + netgen_data = netgen_data.transpose() + shp = netgen_data.shape[:-1] + shp = mesh.comm.bcast(shp, root=0) + if mesh.comm.rank != 0: + netgen_data = np.empty((*shp, 0), dtype=dtype) + + M = FunctionSpace(mesh, "DG", 0) + marked = M.dof_dset.layout_vec.copy() + marked.set(0) + + sfBCInv = sf.createInverse() + section, marked0 = plex.distributeField(sfBCInv, mesh._cell_numbering, marked) + plex_data = None + for i in np.ndindex(shp): + marked0[:netgen_data.shape[-1]] = netgen_data[i] + _, marked = plex.distributeField(sf, section, marked0) + arr = marked.getArray() + if plex_data is None: + plex_data = np.empty(shp + arr.shape, dtype=dtype) + plex_data[i] = arr.astype(dtype) + + plex_data = plex_data.transpose() + else: + plex_data = netgen_data + return plex_data + + +def netgen_to_plex_numbering(mesh): + from firedrake import FunctionSpace + + sf = mesh.sfBC_orig + plex = mesh.topology_dm + cellNum = plex.getCellNumbering().indices + cellNum[cellNum < 0] = -cellNum[cellNum < 0]-1 + fstart, fend = plex.getHeightStratum(0) + cids = list(map(mesh._cell_numbering.getOffset, range(fstart, fend))) + + # Create Netgen to Plex reordering + M = FunctionSpace(mesh, "DG", 0) + marked = M.dof_dset.layout_vec.copy() + marked.set(0) + + cstart, cend = marked.getOwnershipRange() + iperm = cellNum[cids[:cend-cstart]] + marked.setValues(iperm, np.arange(cstart, cend)) + marked.assemble() + marked0 = marked + if sf is not None: + sfBCInv = sf.createInverse() + _, marked0 = plex.distributeField(sfBCInv, mesh._cell_numbering, marked) + + perm = marked0.getArray()[:M.dim()].astype(PETSc.IntType) + return perm, iperm + + @PETSc.Log.EventDecorator() def find_permutation(points_a: npt.NDArray[np.inexact], points_b: npt.NDArray[np.inexact], tol: float = 1e-5): diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 2d5820d4b6..21a054b643 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -211,44 +211,20 @@ def test_firedrake_integral_ball_netgen(): @pytest.mark.skipnetgen +@pytest.mark.parallel([1, 3]) def test_firedrake_integral_sphere_high_order_netgen(): from netgen.csg import CSGeometry, Pnt, Sphere import netgen comm = COMM_WORLD - if comm.Get_rank() == 0: - geo = CSGeometry() - geo.Add(Sphere(Pnt(0, 0, 0), 1).bc("sphere")) - ngmesh = geo.GenerateMesh(maxh=0.1) - else: - ngmesh = netgen.libngpy._meshing.Mesh(3) - - msh = Mesh(ngmesh) - homsh = Mesh(msh.curve_field(4)) - V = FunctionSpace(homsh, "CG", 4) - x, y, z = SpatialCoordinate(homsh) - f = assemble(interpolate(1+0*x, V)) - assert abs(assemble(f * dx) - (4/3)*np.pi) < 1.e-4 - - -@pytest.mark.skipnetgen -@pytest.mark.parallel -def test_firedrake_integral_sphere_high_order_netgen_parallel(): - from netgen.csg import CSGeometry, Pnt, Sphere - import netgen - - comm = COMM_WORLD - if comm.Get_rank() == 0: + if comm.rank == 0: geo = CSGeometry() geo.Add(Sphere(Pnt(0, 0, 0), 1).bc("sphere")) ngmesh = geo.GenerateMesh(maxh=0.7) else: ngmesh = netgen.libngpy._meshing.Mesh(3) - msh = Mesh(ngmesh) - # The default value for location_tol is much too large (see https://github.com/NGSolve/ngsPETSc/issues/76) - # TODO: Once the default value is adjusted this can be removed - homsh = Mesh(msh.curve_field(2, location_tol=1e-8)) + homsh = Mesh(ngmesh, netgen_flags={"degree": 2}) V = FunctionSpace(homsh, "CG", 2) x, y, z = SpatialCoordinate(homsh) f = assemble(interpolate(1+0*x, V)) @@ -257,6 +233,7 @@ def test_firedrake_integral_sphere_high_order_netgen_parallel(): @pytest.mark.skipcomplex @pytest.mark.skipnetgen +@pytest.mark.parallel([1, 3]) def test_firedrake_Adaptivity_netgen(): from netgen.occ import WorkPlane, OCCGeometry, Axes from netgen.occ import X, Z @@ -326,82 +303,8 @@ def adapt(mesh, eta): (eta, error_est) = estimate_error(mesh, uh) error_estimators.append(error_est) dofs.append(uh.function_space().dim()) - mesh = adapt(mesh, eta) - assert error_estimators[-1] < 0.05 - - -@pytest.mark.skipcomplex -@pytest.mark.skipnetgen -@pytest.mark.parallel -def test_firedrake_Adaptivity_netgen_parallel(): - from netgen.occ import WorkPlane, OCCGeometry, Axes - from netgen.occ import X, Z - - def solve_poisson(mesh): - V = FunctionSpace(mesh, "CG", 1) - uh = Function(V, name="Solution") - v = TestFunction(V) - bc = DirichletBC(V, 0, "on_boundary") - f = Constant(1) - F = inner(grad(uh), grad(v))*dx - inner(f, v)*dx - solve(F == 0, uh, bc) - return uh - - def estimate_error(mesh, uh): - W = FunctionSpace(mesh, "DG", 0) - eta_sq = Function(W) - w = TestFunction(W) - f = Constant(1) - h = CellDiameter(mesh) - n = FacetNormal(mesh) - v = CellVolume(mesh) - - # Compute error indicator cellwise - G = inner(eta_sq / v, w)*dx - G = G - inner(h**2 * (f + div(grad(uh)))**2, w) * dx - G = G - inner(h('+')/2 * jump(grad(uh), n)**2, w('+')) * dS - - # Each cell is an independent 1x1 solve, so Jacobi is exact - sp = {"mat_type": "matfree", - "ksp_type": "richardson", - "pc_type": "jacobi"} - solve(G == 0, eta_sq, solver_parameters=sp) - eta = Function(W) - eta.interpolate(sqrt(eta_sq)) # the above computed eta^2 - - with eta.dat.vec_ro as eta_: - error_est = sqrt(eta_.dot(eta_)) - return (eta, error_est) - - def adapt(mesh, eta): - W = FunctionSpace(mesh, "DG", 0) - markers = Function(W) - with eta.dat.vec_ro as eta_: - eta_max = eta_.max()[1] - - theta = 0.5 - should_refine = conditional(gt(eta, theta*eta_max), 1, 0) - markers.interpolate(should_refine) - - refined_mesh = mesh.refine_marked_elements(markers) - return refined_mesh - - rect1 = WorkPlane(Axes((0, 0, 0), n=Z, h=X)).Rectangle(1, 2).Face() - rect2 = WorkPlane(Axes((0, 1, 0), n=Z, h=X)).Rectangle(2, 1).Face() - L = rect1 + rect2 - - geo = OCCGeometry(L, dim=2) - ngmsh = geo.GenerateMesh(maxh=0.1) - mesh = Mesh(ngmsh) - - max_iterations = 10 - error_estimators = [] - dofs = [] - for i in range(max_iterations): - uh = solve_poisson(mesh) - (eta, error_est) = estimate_error(mesh, uh) - error_estimators.append(error_est) - dofs.append(uh.function_space().dim()) + if error_est < 0.05: + break mesh = adapt(mesh, eta) assert error_estimators[-1] < 0.06 From e877b90f2941cf0a11cd850e0615c894fa78f54b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 13 Mar 2026 17:58:58 +0000 Subject: [PATCH 10/25] Fixes --- firedrake/mesh.py | 140 ++++++++++------------ firedrake/mg/netgen.py | 3 - firedrake/netgen.py | 66 ++++++++++ tests/firedrake/regression/test_netgen.py | 109 +---------------- 4 files changed, 134 insertions(+), 184 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 848f254859..8ef5c6cc9c 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2926,6 +2926,8 @@ def refine_marked_elements(self, mark, netgen_flags=None): """ utils.check_netgen_installed() + if not hasattr(self, "netgen_mesh"): + raise ValueError("Adaptive refinement requires a netgen mesh.") if netgen_flags is None: netgen_flags = self.netgen_flags dim = self.geometric_dimension @@ -2947,14 +2949,25 @@ def refine_marked_elements(self, mark, netgen_flags=None): refine_faces = netgen_flags.get("refine_faces", False) if self.comm.rank == 0: max_refs = int(mark_np.max()) - for _ in range(max_refs): - netgen_cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() - netgen_cells.NumPy()["refine"][:mark_np.size] = mark_np > 0 + for r in range(max_refs): + cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() + cells.NumPy()["refine"] = (mark_np[:len(cells)] > 0) if not refine_faces and dim == 3: - netgen_mesh.Elements2D().NumPy()["refine"] = 0 + netgen_mesh.Elements2D().NumPy()["refine"] = False netgen_mesh.Refine(adaptive=True) mark_np -= 1 - + if r < max_refs - 1: + parents = netgen_mesh.parentelements if dim == 3 else netgen_mesh.parentsurfaceelements + parents = parents.NumPy().astype(int).flatten() + cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() + num_coarse_cells = len(cells) + num_fine_cells = parents.shape[0] + indices = np.arange(num_fine_cells, dtype=int) + fine_cells = indices > num_coarse_cells + indices[fine_cells] = parents[indices[fine_cells]] + mark_np = mark_np[indices] + + self.comm.Barrier() return Mesh(netgen_mesh, reorder=self._did_reordering, distribution_parameters=self._distribution_parameters, @@ -2962,7 +2975,7 @@ 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, location_tol=1e-1, cg_field=False): + def curve_field(self, order, permutation_tol=1e-8, cg_field=False): '''Return a function containing the curved coordinates of the mesh. This method requires that the mesh has been constructed from a @@ -2970,34 +2983,30 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F :arg order: the order of the curved mesh. :arg permutation_tol: tolerance used to construct the permutation of the reference element. - :arg location_tol: tolerance used to locate the cell a point belongs to. :arg cg_field: return a CG function field representing the mesh, rather than the default DG field. ''' - import firedrake as fd + from firedrake.netgen import find_permutation, netgen_to_plex_numbering, netgen_distribute + from firedrake.functionspace import VectorFunctionSpace, FunctionSpace + from firedrake.function import Function utils.check_netgen_installed() - - from firedrake.netgen import find_permutation - - netgen_mesh = self.netgen_mesh # Check if the mesh is a surface mesh or two dimensional mesh - if len(netgen_mesh.Elements3D()) == 0: - ng_element = netgen_mesh.Elements2D() + if self.topological_dimension == 2: + ng_element = self.netgen_mesh.Elements2D() else: - ng_element = netgen_mesh.Elements3D() + ng_element = self.netgen_mesh.Elements3D() ng_dimension = len(ng_element) - geom_dim = self.geometric_dimension # Construct the mesh as a Firedrake function if cg_field: - firedrake_space = fd.VectorFunctionSpace(self, "CG", order) + coords_space = VectorFunctionSpace(self, "CG", order) else: low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] ufl_element = low_order_element.reconstruct(degree=order) - firedrake_space = fd.VectorFunctionSpace(self, fd.BrokenElement(ufl_element)) - new_coordinates = fd.assemble(fd.interpolate(self.coordinates, firedrake_space)) + coords_space = VectorFunctionSpace(self, finat.ufl.BrokenElement(ufl_element)) + new_coordinates = Function(coords_space).interpolate(self.coordinates) # Compute reference points using fiat fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent @@ -3012,75 +3021,50 @@ def curve_field(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=F ref.append(pt) reference_space_points = np.array(ref) - # Curve the mesh on rank 0 only - if self.comm.rank == 0: - # Construct numpy arrays for physical domain data - physical_space_points = np.zeros( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.zeros( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) - # NOTE: This will segfault for CSG! - netgen_mesh.Curve(order) - netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) - curved = ng_element.NumPy()["curved"] - - # Broadcast a boolean array identifying curved cells - curved = self.comm.bcast(curved, root=0) - physical_space_points = physical_space_points[curved] - curved_space_points = curved_space_points[curved] - else: - curved = self.comm.bcast(None, root=0) - # Construct numpy arrays as buffers to receive physical domain data - ncurved = np.sum(curved) - physical_space_points = np.zeros( - (ncurved, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.zeros( - (ncurved, reference_space_points.shape[0], geom_dim) - ) + # Construct numpy arrays for physical domain data + physical_space_points = np.zeros( + (ng_dimension, reference_space_points.shape[0], self.geometric_dimension) + ) + curved_space_points = np.zeros( + (ng_dimension, reference_space_points.shape[0], self.geometric_dimension) + ) + self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) + # NOTE: This will segfault for CSG! + self.netgen_mesh.Curve(order) + self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) + curved = ng_element.NumPy()["curved"] + + # Get numbering + DG0 = FunctionSpace(self, "DG", 0) + rstart, rend = DG0.dof_dset.layout_vec.getOwnershipRange() + num_cells = rend - rstart + _, iperm = netgen_to_plex_numbering(self) + iperm -= rstart + + # Distribute curved cell data + own_curved = netgen_distribute(self, curved) + own_curved = np.flatnonzero(own_curved[:num_cells]) - # Broadcast curved cell point data - physical_space_points = self.comm.bcast(physical_space_points, root=0) - curved_space_points = self.comm.bcast(curved_space_points, root=0) cell_node_map = new_coordinates.cell_node_map() + pyop2_index = cell_node_map.values[iperm[own_curved]] - # Select only the points in curved cells - barycentres = np.average(physical_space_points, axis=1) - ng_index = list(map(lambda x: self.locate_cell(x, tolerance=location_tol), barycentres)) - - # Select only the indices of cells owned by this rank - owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] + # Distribute coordinate data + own_curved_points = netgen_distribute(self, curved_space_points)[own_curved] + own_physical_points = netgen_distribute(self, physical_space_points)[own_curved] - # Get the PyOP2 indices corresponding to the netgen indices - ng_index = [idx for idx, o in zip(ng_index, owned) if o] - pyop2_index = cell_node_map.values[ng_index].flatten() - - # Select only the points owned by this rank - physical_space_points = physical_space_points[owned] - curved_space_points = curved_space_points[owned] - barycentres = barycentres[owned] - - if any(owned): + if any(own_curved): # Find the correct coordinate permutation for each cell permutation = find_permutation( - physical_space_points, - new_coordinates.dat.data[pyop2_index].real.reshape( - physical_space_points.shape - ), - tol=permutation_tol + own_physical_points, + new_coordinates.dat.data[pyop2_index].real, + tol=permutation_tol, ) - # Apply the permutation to each cell in turn - for ii, p in enumerate(curved_space_points): - curved_space_points[ii] = p[permutation[ii]] - else: - print("barf") + for ii, p in enumerate(own_curved_points): + own_curved_points[ii] = p[permutation[ii]] # Assign the curved coordinates to the dat - new_coordinates.dat.data[pyop2_index] = curved_space_points.reshape(-1, geom_dim) + new_coordinates.dat.data[pyop2_index] = own_curved_points return new_coordinates diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 9356938de1..516e94a956 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -256,7 +256,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): cg = flags.get("cg", False) nested = flags.get("nested", snap in ["coarse"]) permutation_tol = flags.get("permutation_tol", 1e-8) - location_tol = flags.get("location_tol", 1e-8) # Firedrake quantities meshes = [] lgmaps = [] @@ -264,7 +263,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): if mesh.coordinates.function_space().ufl_element().degree() != order[0]: coordinates = mesh.curve_field( order=order[0], - location_tol=location_tol, permutation_tol=permutation_tol, cg_field=cg, ) @@ -329,7 +327,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): if snap == "geometry": coordinates = mesh.curve_field( order=order[l+1], - location_tol=location_tol, permutation_tol=permutation_tol, cg_field=cg, ) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 83e6476a20..f908edbe07 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -29,6 +29,72 @@ class comp: Mesh = type(None) +def netgen_distribute(mesh, netgen_data): + from firedrake import FunctionSpace + # Create Netgen to Plex reordering + plex = mesh.topology_dm + sf = mesh.sfBC_orig + perm, iperm = netgen_to_plex_numbering(mesh) + if sf is not None: + netgen_data = np.asarray(netgen_data) + dtype = netgen_data.dtype + dtype = mesh.comm.bcast(dtype, root=0) + + netgen_data = netgen_data.transpose() + shp = netgen_data.shape[:-1] + shp = mesh.comm.bcast(shp, root=0) + if mesh.comm.rank != 0: + netgen_data = np.empty((*shp, 0), dtype=dtype) + + M = FunctionSpace(mesh, "DG", 0) + marked = M.dof_dset.layout_vec.copy() + marked.set(0) + + sfBCInv = sf.createInverse() + section, marked0 = plex.distributeField(sfBCInv, mesh._cell_numbering, marked) + plex_data = None + for i in np.ndindex(shp): + marked0[:netgen_data.shape[-1]] = netgen_data[i] + _, marked = plex.distributeField(sf, section, marked0) + arr = marked.getArray() + if plex_data is None: + plex_data = np.empty(shp + arr.shape, dtype=dtype) + plex_data[i] = arr.astype(dtype) + + plex_data = plex_data.transpose() + else: + plex_data = netgen_data + return plex_data + + +def netgen_to_plex_numbering(mesh): + from firedrake import FunctionSpace + + sf = mesh.sfBC_orig + plex = mesh.topology_dm + cellNum = plex.getCellNumbering().indices + cellNum[cellNum < 0] = -cellNum[cellNum < 0]-1 + fstart, fend = plex.getHeightStratum(0) + cids = list(map(mesh._cell_numbering.getOffset, range(fstart, fend))) + + # Create Netgen to Plex reordering + M = FunctionSpace(mesh, "DG", 0) + marked = M.dof_dset.layout_vec.copy() + marked.set(0) + + cstart, cend = marked.getOwnershipRange() + iperm = cellNum[cids[:cend-cstart]] + marked.setValues(iperm, np.arange(cstart, cend)) + marked.assemble() + marked0 = marked + if sf is not None: + sfBCInv = sf.createInverse() + _, marked0 = plex.distributeField(sfBCInv, mesh._cell_numbering, marked) + + perm = marked0.getArray()[:M.dim()].astype(PETSc.IntType) + return perm, iperm + + @PETSc.Log.EventDecorator() def find_permutation(points_a: npt.NDArray[np.inexact], points_b: npt.NDArray[np.inexact], tol: float = 1e-5): diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 2d5820d4b6..21a054b643 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -211,44 +211,20 @@ def test_firedrake_integral_ball_netgen(): @pytest.mark.skipnetgen +@pytest.mark.parallel([1, 3]) def test_firedrake_integral_sphere_high_order_netgen(): from netgen.csg import CSGeometry, Pnt, Sphere import netgen comm = COMM_WORLD - if comm.Get_rank() == 0: - geo = CSGeometry() - geo.Add(Sphere(Pnt(0, 0, 0), 1).bc("sphere")) - ngmesh = geo.GenerateMesh(maxh=0.1) - else: - ngmesh = netgen.libngpy._meshing.Mesh(3) - - msh = Mesh(ngmesh) - homsh = Mesh(msh.curve_field(4)) - V = FunctionSpace(homsh, "CG", 4) - x, y, z = SpatialCoordinate(homsh) - f = assemble(interpolate(1+0*x, V)) - assert abs(assemble(f * dx) - (4/3)*np.pi) < 1.e-4 - - -@pytest.mark.skipnetgen -@pytest.mark.parallel -def test_firedrake_integral_sphere_high_order_netgen_parallel(): - from netgen.csg import CSGeometry, Pnt, Sphere - import netgen - - comm = COMM_WORLD - if comm.Get_rank() == 0: + if comm.rank == 0: geo = CSGeometry() geo.Add(Sphere(Pnt(0, 0, 0), 1).bc("sphere")) ngmesh = geo.GenerateMesh(maxh=0.7) else: ngmesh = netgen.libngpy._meshing.Mesh(3) - msh = Mesh(ngmesh) - # The default value for location_tol is much too large (see https://github.com/NGSolve/ngsPETSc/issues/76) - # TODO: Once the default value is adjusted this can be removed - homsh = Mesh(msh.curve_field(2, location_tol=1e-8)) + homsh = Mesh(ngmesh, netgen_flags={"degree": 2}) V = FunctionSpace(homsh, "CG", 2) x, y, z = SpatialCoordinate(homsh) f = assemble(interpolate(1+0*x, V)) @@ -257,6 +233,7 @@ def test_firedrake_integral_sphere_high_order_netgen_parallel(): @pytest.mark.skipcomplex @pytest.mark.skipnetgen +@pytest.mark.parallel([1, 3]) def test_firedrake_Adaptivity_netgen(): from netgen.occ import WorkPlane, OCCGeometry, Axes from netgen.occ import X, Z @@ -326,82 +303,8 @@ def adapt(mesh, eta): (eta, error_est) = estimate_error(mesh, uh) error_estimators.append(error_est) dofs.append(uh.function_space().dim()) - mesh = adapt(mesh, eta) - assert error_estimators[-1] < 0.05 - - -@pytest.mark.skipcomplex -@pytest.mark.skipnetgen -@pytest.mark.parallel -def test_firedrake_Adaptivity_netgen_parallel(): - from netgen.occ import WorkPlane, OCCGeometry, Axes - from netgen.occ import X, Z - - def solve_poisson(mesh): - V = FunctionSpace(mesh, "CG", 1) - uh = Function(V, name="Solution") - v = TestFunction(V) - bc = DirichletBC(V, 0, "on_boundary") - f = Constant(1) - F = inner(grad(uh), grad(v))*dx - inner(f, v)*dx - solve(F == 0, uh, bc) - return uh - - def estimate_error(mesh, uh): - W = FunctionSpace(mesh, "DG", 0) - eta_sq = Function(W) - w = TestFunction(W) - f = Constant(1) - h = CellDiameter(mesh) - n = FacetNormal(mesh) - v = CellVolume(mesh) - - # Compute error indicator cellwise - G = inner(eta_sq / v, w)*dx - G = G - inner(h**2 * (f + div(grad(uh)))**2, w) * dx - G = G - inner(h('+')/2 * jump(grad(uh), n)**2, w('+')) * dS - - # Each cell is an independent 1x1 solve, so Jacobi is exact - sp = {"mat_type": "matfree", - "ksp_type": "richardson", - "pc_type": "jacobi"} - solve(G == 0, eta_sq, solver_parameters=sp) - eta = Function(W) - eta.interpolate(sqrt(eta_sq)) # the above computed eta^2 - - with eta.dat.vec_ro as eta_: - error_est = sqrt(eta_.dot(eta_)) - return (eta, error_est) - - def adapt(mesh, eta): - W = FunctionSpace(mesh, "DG", 0) - markers = Function(W) - with eta.dat.vec_ro as eta_: - eta_max = eta_.max()[1] - - theta = 0.5 - should_refine = conditional(gt(eta, theta*eta_max), 1, 0) - markers.interpolate(should_refine) - - refined_mesh = mesh.refine_marked_elements(markers) - return refined_mesh - - rect1 = WorkPlane(Axes((0, 0, 0), n=Z, h=X)).Rectangle(1, 2).Face() - rect2 = WorkPlane(Axes((0, 1, 0), n=Z, h=X)).Rectangle(2, 1).Face() - L = rect1 + rect2 - - geo = OCCGeometry(L, dim=2) - ngmsh = geo.GenerateMesh(maxh=0.1) - mesh = Mesh(ngmsh) - - max_iterations = 10 - error_estimators = [] - dofs = [] - for i in range(max_iterations): - uh = solve_poisson(mesh) - (eta, error_est) = estimate_error(mesh, uh) - error_estimators.append(error_est) - dofs.append(uh.function_space().dim()) + if error_est < 0.05: + break mesh = adapt(mesh, eta) assert error_estimators[-1] < 0.06 From d1b8c4c1e47c9c4ac77ffd931aef3aace64d5a36 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Mar 2026 11:21:45 +0000 Subject: [PATCH 11/25] Fixes --- firedrake/mesh.py | 84 ++++++++++---------- firedrake/mg/netgen.py | 26 +++--- firedrake/netgen.py | 74 ++++++----------- tests/firedrake/multigrid/test_netgen_gmg.py | 25 +++--- 4 files changed, 96 insertions(+), 113 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 8ef5c6cc9c..b4f3a897e8 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2975,7 +2975,7 @@ 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=False): + def curve_field(self, order, permutation_tol=1e-8, cg_field=None): '''Return a function containing the curved coordinates of the mesh. This method requires that the mesh has been constructed from a @@ -2983,15 +2983,19 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=False): :arg order: the order of the curved mesh. :arg permutation_tol: tolerance used to construct the permutation of the reference element. - :arg cg_field: return a CG function field representing the mesh, rather than the - default DG field. + :arg cg_field: return a CG function field representing the mesh, as opposed to a DG field. ''' - from firedrake.netgen import find_permutation, netgen_to_plex_numbering, netgen_distribute - from firedrake.functionspace import VectorFunctionSpace, FunctionSpace - from firedrake.function import Function utils.check_netgen_installed() + + from firedrake.netgen import find_permutation, plex_to_netgen_numbering, netgen_distribute + from firedrake.functionspace import FunctionSpace + from firedrake.function import Function + + if cg_field is None: + cg_field = not self.coordinates.function_space().finat_element.is_dg() + # Check if the mesh is a surface mesh or two dimensional mesh if self.topological_dimension == 2: ng_element = self.netgen_mesh.Elements2D() @@ -2999,26 +3003,22 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=False): ng_element = self.netgen_mesh.Elements3D() ng_dimension = len(ng_element) - # Construct the mesh as a Firedrake function - if cg_field: - coords_space = VectorFunctionSpace(self, "CG", order) - else: - low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] - ufl_element = low_order_element.reconstruct(degree=order) - coords_space = VectorFunctionSpace(self, finat.ufl.BrokenElement(ufl_element)) + # Construct the coordinates as a Firedrake function + low_order_element = self.coordinates.function_space().ufl_element() + ufl_element = low_order_element.reconstruct(degree=order) + if not cg_field: + ufl_element = finat.ufl.BrokenElement(ufl_element) + coords_space = FunctionSpace(self, ufl_element) new_coordinates = Function(coords_space).interpolate(self.coordinates) # Compute reference points using fiat fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent - entity_ids = fiat_element.entity_dofs() nodes = fiat_element.dual_basis() ref = [] - for dim in entity_ids: - for entity in entity_ids[dim]: - for dof in entity_ids[dim][entity]: - # Assert singleton point for each node. - pt, = nodes[dof].get_point_dict().keys() - ref.append(pt) + for node in nodes: + # Assert singleton point for each node. + pt, = node.get_point_dict().keys() + ref.append(pt) reference_space_points = np.array(ref) # Construct numpy arrays for physical domain data @@ -3034,37 +3034,35 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=False): self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) curved = ng_element.NumPy()["curved"] - # Get numbering - DG0 = FunctionSpace(self, "DG", 0) - rstart, rend = DG0.dof_dset.layout_vec.getOwnershipRange() - num_cells = rend - rstart - _, iperm = netgen_to_plex_numbering(self) - iperm -= rstart - # Distribute curved cell data - own_curved = netgen_distribute(self, curved) + cell_node_map = new_coordinates.cell_node_map() + num_cells = cell_node_map.values.shape[0] + DG0 = FunctionSpace(self, "DG", 0) + own_curved = netgen_distribute(DG0, curved) own_curved = np.flatnonzero(own_curved[:num_cells]) - cell_node_map = new_coordinates.cell_node_map() + # Get numbering + iperm = plex_to_netgen_numbering(self) pyop2_index = cell_node_map.values[iperm[own_curved]] # Distribute coordinate data - own_curved_points = netgen_distribute(self, curved_space_points)[own_curved] - own_physical_points = netgen_distribute(self, physical_space_points)[own_curved] - - if any(own_curved): - # Find the correct coordinate permutation for each cell - permutation = find_permutation( - own_physical_points, - new_coordinates.dat.data[pyop2_index].real, - tol=permutation_tol, - ) - # Apply the permutation to each cell in turn - for ii, p in enumerate(own_curved_points): - own_curved_points[ii] = p[permutation[ii]] + broken_space = coords_space.broken_space() + own_curved_points = netgen_distribute(broken_space, curved_space_points)[own_curved] + own_physical_points = netgen_distribute(broken_space, physical_space_points)[own_curved] + + # Find the correct coordinate permutation for each cell + permutation = find_permutation( + own_physical_points, + new_coordinates.dat.data_ro_with_halos[pyop2_index].real, + tol=permutation_tol, + ) + self.comm.Barrier() + # Apply the permutation to each cell in turn + for i in range(own_curved_points.shape[0]): + own_curved_points[i, ...] = own_curved_points[i, permutation[i], ...] # Assign the curved coordinates to the dat - new_coordinates.dat.data[pyop2_index] = own_curved_points + new_coordinates.dat.data_wo_with_halos[pyop2_index] = own_curved_points return new_coordinates diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 516e94a956..e94abbf2af 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -41,14 +41,15 @@ def snapToNetgenDMPlex(ngmesh, petscPlex, comm): ''' logger.info(f"\t\t\t[{time.time()}]Snapping the DMPlex to NETGEN mesh") if len(ngmesh.Elements3D()) == 0: - ng_coelement = ngmesh.Elements1D + ng_coelement = ngmesh.Elements1D() else: - ng_coelement = ngmesh.Elements2D - if comm.rank == 0: - nodes_to_correct = ng_coelement().NumPy()["nodes"] + ng_coelement = ngmesh.Elements2D() + + nodes_to_correct = ng_coelement.NumPy()["nodes"] + sf = None + if sf is not None: nodes_to_correct = comm.bcast(nodes_to_correct, root=0) - else: - nodes_to_correct = comm.bcast(None, root=0) + logger.info(f"\t\t\t[{time.time()}]Point distributed") nodes_to_correct = trim_util(nodes_to_correct) nodes_to_correct_sorted = np.hstack(nodes_to_correct.reshape((-1, 1))) @@ -253,7 +254,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): snap = flags.get("snap_to", "geometry") snap_smoothing = flags.get("snap_smoothing", "hyperelastic") logger.info(f"\tSnap to {snap} using {snap_smoothing} smoothing (if snapping to coarse)") - cg = flags.get("cg", False) + cg = flags.get("cg", not mesh.coordinates.function_space().finat_element.is_dg()) nested = flags.get("nested", snap in ["coarse"]) permutation_tol = flags.get("permutation_tol", 1e-8) # Firedrake quantities @@ -289,6 +290,10 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm) cdm = rdm + # Snap the mesh to the Netgen mesh + if snap == "geometry": + snapToNetgenDMPlex(ngmesh, rdm, comm) + if optMoves: # Optimises the mesh, for example smoothing if ngmesh.dim == 2: @@ -298,10 +303,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): else: raise ValueError("Only 2D and 3D meshes can be optimised.") - # Snap the mesh to the Netgen mesh - if snap == "geometry": - snapToNetgenDMPlex(ngmesh, rdm, comm) - # We construct a Firedrake mesh from the DMPlex mesh parameters = {} if distribution_parameters is not None: @@ -333,7 +334,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): mesh = reconstruct_mesh(mesh, coordinates) elif snap == "coarse": mesh = snapToCoarse(meshes[0].coordinates, mesh, order[l+1], snap_smoothing, cg) - toc = time.time() logger.info(f"\t\t\tMeshed curved. Time taken: {toc-tic}") logger.info(f"\t\tLevel {l+1}: with {ngmesh.Coordinates().shape[0]}\ @@ -344,6 +344,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): # Populate the coarse to fine map coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes, lgmaps) + return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, 1, nested=nested) @@ -361,4 +362,5 @@ def reconstruct_mesh(mesh, *args, **kwargs): tmesh._did_reordering = mesh._did_reordering tmesh.netgen_mesh = mesh.netgen_mesh tmesh.netgen_flags = mesh.netgen_flags + tmesh.sfBC_orig = mesh.sfBC_orig return tmesh diff --git a/firedrake/netgen.py b/firedrake/netgen.py index f908edbe07..f7b852eaa0 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -29,70 +29,45 @@ class comp: Mesh = type(None) -def netgen_distribute(mesh, netgen_data): - from firedrake import FunctionSpace - # Create Netgen to Plex reordering +def netgen_distribute(V, netgen_data): + mesh = V.mesh() plex = mesh.topology_dm sf = mesh.sfBC_orig - perm, iperm = netgen_to_plex_numbering(mesh) - if sf is not None: - netgen_data = np.asarray(netgen_data) - dtype = netgen_data.dtype - dtype = mesh.comm.bcast(dtype, root=0) - - netgen_data = netgen_data.transpose() - shp = netgen_data.shape[:-1] - shp = mesh.comm.bcast(shp, root=0) - if mesh.comm.rank != 0: - netgen_data = np.empty((*shp, 0), dtype=dtype) - - M = FunctionSpace(mesh, "DG", 0) - marked = M.dof_dset.layout_vec.copy() - marked.set(0) + fstart, fend = plex.getHeightStratum(0) + netgen_data = np.asarray(netgen_data) + nshape = netgen_data.shape + dtype = netgen_data.dtype + shp = V.shape + if sf is not None: + section = V.dm.getDefaultSection() + marked = V.dof_dset.layout_vec sfBCInv = sf.createInverse() - section, marked0 = plex.distributeField(sfBCInv, mesh._cell_numbering, marked) + section0, marked0 = plex.distributeField(sfBCInv, section, marked) plex_data = None + + marked0.set(0) + d = netgen_data for i in np.ndindex(shp): - marked0[:netgen_data.shape[-1]] = netgen_data[i] - _, marked = plex.distributeField(sf, section, marked0) + di = d[(..., *i)].flatten() + marked0[:len(di)] = di + _, marked = plex.distributeField(sf, section0, marked0) arr = marked.getArray() if plex_data is None: - plex_data = np.empty(shp + arr.shape, dtype=dtype) - plex_data[i] = arr.astype(dtype) - - plex_data = plex_data.transpose() + plex_data = np.empty(arr.shape + shp, dtype=dtype) + plex_data[(..., *i)] = arr else: plex_data = netgen_data + plex_data = plex_data.reshape(-1, *nshape[1:]) return plex_data -def netgen_to_plex_numbering(mesh): - from firedrake import FunctionSpace - - sf = mesh.sfBC_orig +def plex_to_netgen_numbering(mesh): plex = mesh.topology_dm - cellNum = plex.getCellNumbering().indices - cellNum[cellNum < 0] = -cellNum[cellNum < 0]-1 fstart, fend = plex.getHeightStratum(0) cids = list(map(mesh._cell_numbering.getOffset, range(fstart, fend))) - - # Create Netgen to Plex reordering - M = FunctionSpace(mesh, "DG", 0) - marked = M.dof_dset.layout_vec.copy() - marked.set(0) - - cstart, cend = marked.getOwnershipRange() - iperm = cellNum[cids[:cend-cstart]] - marked.setValues(iperm, np.arange(cstart, cend)) - marked.assemble() - marked0 = marked - if sf is not None: - sfBCInv = sf.createInverse() - _, marked0 = plex.distributeField(sfBCInv, mesh._cell_numbering, marked) - - perm = marked0.getArray()[:M.dim()].astype(PETSc.IntType) - return perm, iperm + iperm = np.asarray(cids) + return iperm @PETSc.Log.EventDecorator() @@ -112,6 +87,9 @@ def find_permutation(points_a: npt.NDArray[np.inexact], points_b: npt.NDArray[np 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)] + if len(p) == 0: + return p + try: permutation = np.array(p, ndmin=2) except ValueError as e: diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 3e788318b0..8274f32012 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -2,7 +2,7 @@ from firedrake import * -@pytest.fixture(params=[2, 3]) +@pytest.fixture(params=[3, 2]) def ngmesh(request): dim = request.param if dim == 2: @@ -23,20 +23,25 @@ def ngmesh(request): @pytest.mark.skipcomplex @pytest.mark.skipnetgen -# @pytest.mark.parallel([1, 3]) -@pytest.mark.parametrize("netgen_degree", [1, 3, (1, 2, 3)]) +@pytest.mark.parallel([1, 3]) +@pytest.mark.parametrize("netgen_degree", [1, 3, (1, 2, 3)], ids=lambda degree: f"{degree=}") def test_netgen_mg(ngmesh, netgen_degree): dparams = {"overlap_type": (DistributedMeshOverlapType.VERTEX, 1)} - mesh = Mesh(ngmesh, distribution_parameters=dparams) - nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": netgen_degree}) - mesh = nh[-1] + base = Mesh(ngmesh, distribution_parameters=dparams) + mh = MeshHierarchy(base, 2, netgen_flags={"degree": netgen_degree}) + mesh = mh[-1] try: len(netgen_degree) except TypeError: - netgen_degree = (netgen_degree,)*len(nh) + netgen_degree = (netgen_degree,)*len(mh) - for m, deg in zip(nh, netgen_degree): - assert m.coordinates.function_space().ufl_element().degree() == deg + coords_space = base.coordinates.function_space() + assert coords_space.ufl_element().degree() == 1 + assert not coords_space.finat_element.is_dg() + for m, deg in zip(mh, netgen_degree): + coords_space = m.coordinates.function_space() + assert coords_space.ufl_element().degree() == deg + assert not coords_space.finat_element.is_dg() V = FunctionSpace(mesh, "CG", 3) u = TrialFunction(V) @@ -57,7 +62,7 @@ def test_netgen_mg(ngmesh, netgen_degree): solve(a == L, uh, bcs=bcs, solver_parameters={ "ksp_type": "cg", "ksp_norm_type": "natural", - "ksp_max_it": 12, + "ksp_max_it": 14, "ksp_rtol": rtol, "ksp_monitor": None, "pc_type": "mg", From 8f404e1cafce0de18885f1963a34e1455ee098e0 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Mar 2026 11:38:21 +0000 Subject: [PATCH 12/25] tidy --- firedrake/mesh.py | 12 +++++------- firedrake/mg/netgen.py | 15 ++++----------- firedrake/netgen.py | 23 ++++++++++------------- 3 files changed, 19 insertions(+), 31 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index b4f3a897e8..100725c257 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3004,11 +3004,10 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): ng_dimension = len(ng_element) # Construct the coordinates as a Firedrake function - low_order_element = self.coordinates.function_space().ufl_element() - ufl_element = low_order_element.reconstruct(degree=order) + coords_space = self.coordinates.function_space().reconstruct(degree=order) + broken_space = coords_space.broken_space() if not cg_field: - ufl_element = finat.ufl.BrokenElement(ufl_element) - coords_space = FunctionSpace(self, ufl_element) + coords_space = broken_space new_coordinates = Function(coords_space).interpolate(self.coordinates) # Compute reference points using fiat @@ -3042,11 +3041,10 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): own_curved = np.flatnonzero(own_curved[:num_cells]) # Get numbering - iperm = plex_to_netgen_numbering(self) - pyop2_index = cell_node_map.values[iperm[own_curved]] + cellNum = plex_to_netgen_numbering(self) + pyop2_index = cell_node_map.values[cellNum[own_curved]] # Distribute coordinate data - broken_space = coords_space.broken_space() own_curved_points = netgen_distribute(broken_space, curved_space_points)[own_curved] own_physical_points = netgen_distribute(broken_space, physical_space_points)[own_curved] diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index e94abbf2af..a788cffa16 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -45,12 +45,10 @@ def snapToNetgenDMPlex(ngmesh, petscPlex, comm): else: ng_coelement = ngmesh.Elements2D() + # When we create a netgen mesh from a refined plex, + # the netgen mesh represents the local submesh. + # Therefore, there is no need to distribute the netgen data nodes_to_correct = ng_coelement.NumPy()["nodes"] - sf = None - if sf is not None: - nodes_to_correct = comm.bcast(nodes_to_correct, root=0) - - logger.info(f"\t\t\t[{time.time()}]Point distributed") nodes_to_correct = trim_util(nodes_to_correct) nodes_to_correct_sorted = np.hstack(nodes_to_correct.reshape((-1, 1))) nodes_to_correct_sorted.sort() @@ -279,21 +277,18 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): no = impl.create_lgmap(cdm) o = impl.create_lgmap(mesh.topology_dm) lgmaps.append((no, o)) - mesh.topology_dm.setRefineLevel(0) - ngmesh = mesh.netgen_mesh meshes.append(mesh) + ngmesh = mesh.netgen_mesh for l in range(levs): # Straighten the mesh ngmesh.Curve(1) rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm) cdm = rdm - # Snap the mesh to the Netgen mesh if snap == "geometry": snapToNetgenDMPlex(ngmesh, rdm, comm) - if optMoves: # Optimises the mesh, for example smoothing if ngmesh.dim == 2: @@ -341,10 +336,8 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): and optimisation moves {optMoves}.") mesh.topology_dm.setRefineLevel(l+1) meshes.append(mesh) - # Populate the coarse to fine map coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes, lgmaps) - return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, 1, nested=nested) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index f7b852eaa0..4e64862b32 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -30,16 +30,14 @@ class comp: def netgen_distribute(V, netgen_data): + netgen_data = np.asarray(netgen_data) mesh = V.mesh() - plex = mesh.topology_dm sf = mesh.sfBC_orig - fstart, fend = plex.getHeightStratum(0) - - netgen_data = np.asarray(netgen_data) - nshape = netgen_data.shape - dtype = netgen_data.dtype - shp = V.shape if sf is not None: + plex = mesh.topology_dm + nshape = netgen_data.shape + dtype = netgen_data.dtype + section = V.dm.getDefaultSection() marked = V.dof_dset.layout_vec sfBCInv = sf.createInverse() @@ -48,26 +46,25 @@ def netgen_distribute(V, netgen_data): marked0.set(0) d = netgen_data - for i in np.ndindex(shp): + for i in np.ndindex(V.shape): di = d[(..., *i)].flatten() marked0[:len(di)] = di _, marked = plex.distributeField(sf, section0, marked0) arr = marked.getArray() if plex_data is None: - plex_data = np.empty(arr.shape + shp, dtype=dtype) + plex_data = np.empty(arr.shape + V.shape, dtype=dtype) plex_data[(..., *i)] = arr + plex_data = plex_data.reshape(-1, *nshape[1:]) else: plex_data = netgen_data - plex_data = plex_data.reshape(-1, *nshape[1:]) return plex_data def plex_to_netgen_numbering(mesh): plex = mesh.topology_dm fstart, fend = plex.getHeightStratum(0) - cids = list(map(mesh._cell_numbering.getOffset, range(fstart, fend))) - iperm = np.asarray(cids) - return iperm + cellNum = list(map(mesh._cell_numbering.getOffset, range(fstart, fend))) + return np.asarray(cellNum) @PETSc.Log.EventDecorator() From 925072712b0be190ff1a7e8f87a5ddec13a9591b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Mar 2026 11:43:13 +0000 Subject: [PATCH 13/25] tidy --- firedrake/mesh.py | 2 +- firedrake/mg/netgen.py | 12 ++++++------ firedrake/netgen.py | 2 +- firedrake/pointquery_utils.py | 7 ++----- 4 files changed, 10 insertions(+), 13 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 100725c257..cd99d24e58 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3057,7 +3057,7 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): self.comm.Barrier() # Apply the permutation to each cell in turn for i in range(own_curved_points.shape[0]): - own_curved_points[i, ...] = own_curved_points[i, permutation[i], ...] + own_curved_points[i] = own_curved_points[i, permutation[i]] # Assign the curved coordinates to the dat new_coordinates.dat.data_wo_with_halos[pyop2_index] = own_curved_points diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index a788cffa16..595181efad 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -241,25 +241,25 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): comm = mesh.comm # Parse netgen flags if not isinstance(flags, dict): - flags = {} + flags = mesh.netgen_flags order = flags.get("degree", 1) - logger.info(f"\tOrder of the hierarchy: {order}") if isinstance(order, int): order = [order]*(levs+1) + permutation_tol = flags.get("permutation_tol", 1e-8) refType = flags.get("refinement_type", "uniform") - logger.info(f"\tRefinement type: {refType}") optMoves = flags.get("optimisation_moves", False) snap = flags.get("snap_to", "geometry") snap_smoothing = flags.get("snap_smoothing", "hyperelastic") - logger.info(f"\tSnap to {snap} using {snap_smoothing} smoothing (if snapping to coarse)") cg = flags.get("cg", not mesh.coordinates.function_space().finat_element.is_dg()) nested = flags.get("nested", snap in ["coarse"]) - permutation_tol = flags.get("permutation_tol", 1e-8) + logger.info(f"\tOrder of the hierarchy: {order}") + logger.info(f"\tRefinement type: {refType}") + logger.info(f"\tSnap to {snap} using {snap_smoothing} smoothing (if snapping to coarse)") # Firedrake quantities meshes = [] lgmaps = [] # Curve the mesh - if mesh.coordinates.function_space().ufl_element().degree() != order[0]: + if order[0] != mesh.coordinates.function_space().ufl_element().degree(): coordinates = mesh.curve_field( order=order[0], permutation_tol=permutation_tol, diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 4e64862b32..ea7b3e0b68 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -135,7 +135,7 @@ class FiredrakeMesh: def __init__(self, mesh, netgen_flags=None, user_comm=fd.COMM_WORLD): self.comm = user_comm # Parsing netgen flags - if netgen_flags is None: + if not isinstance(netgen_flags, dict): netgen_flags = {} split2tets = netgen_flags.get("split_to_tets", False) split = netgen_flags.get("split", False) diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 4102d0f020..141eb3a939 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -128,10 +128,7 @@ def init_X(fiat_cell, parameters): @PETSc.Log.EventDecorator() -def to_reference_coords_newton_step(ufl_coordinate_element, parameters, - x0_dtype="double", - dX_dtype=ScalarType, - kernel_name="to_reference_coords_newton_step"): +def to_reference_coords_newton_step(ufl_coordinate_element, parameters, x0_dtype="double", dX_dtype=ScalarType): # Set up UFL form cell = ufl_coordinate_element.cell domain = ufl.Mesh(ufl_coordinate_element) @@ -194,7 +191,7 @@ def predicate(index): impero_c = impero_utils.compile_gem(assignments, ()) kernel, _ = tsfc.loopy.generate( impero_c, loopy_args, ScalarType, - kernel_name=kernel_name) + kernel_name="to_reference_coords_newton_step") return lp.generate_code_v2(kernel).device_code() From f1cd773d58786fa15222962f9600297b6cbe2d3d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Mar 2026 12:09:43 +0000 Subject: [PATCH 14/25] update demo --- demos/netgen/netgen_mesh.py.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 480220cbd0..2f4439187c 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -225,7 +225,7 @@ It is now time to define the solve, mark and refine loop that is at the heart of for i in range(max_iterations): - print("level {}".format(i)) + print(f"level {i}") lam, uh, V = Solve(msh, labels) mark = Mark(msh, uh, lam) msh = msh.refine_marked_elements(mark) @@ -332,8 +332,7 @@ As usual, we generate a mesh for the described geometry and use the Firedrake-Ne High-order Meshes ------------------ It is possible to construct high-order meshes for a geometry constructed in Netgen. -In order to do so we need to use the ``curve_field`` method of a Firedrake ``Mesh`` object generated from a Netgen mesh. -In particular, we need to pass the degree of the polynomial field we want to use to parametrise the coordinates of the domain to the ``curve_field`` method, which will return a ``Function`` constructed on a DG space for this purpose. :: +We can set the degree of the geometry via ``netgen_flags`` keyword argument of the ``Mesh`` constructor. :: from netgen.occ import WorkPlane, OCCGeometry import netgen @@ -347,7 +346,7 @@ In particular, we need to pass the degree of the polynomial field we want to use ngmesh = OCCGeometry(shape,dim=2).GenerateMesh(maxh=1.) else: ngmesh = netgen.libngpy._meshing.Mesh(2) - mesh = Mesh(Mesh(ngmesh, comm=COMM_WORLD).curve_field(4)) + mesh = Mesh(ngmesh, comm=COMM_WORLD, netgen_flags={"degree": 3}) VTKFile("output/MeshExample5.pvd").write(mesh) .. figure:: Example5.png @@ -366,7 +365,8 @@ We will now show how to solve the Poisson problem on a high-order mesh, of order ngmesh = OCCGeometry(shape,dim=3).GenerateMesh(maxh=1.) else: ngmesh = netgen.libngpy._meshing.Mesh(3) - mesh = Mesh(Mesh(ngmesh).curve_field(3)) + mesh = Mesh(ngmesh, netgen_flags={"degree": 3}) + # Solving the Poisson problem VTKFile("output/MeshExample6.pvd").write(mesh) x, y, z = SpatialCoordinate(mesh) @@ -405,7 +405,7 @@ It is also possible to construct high-order meshes using the ``SplineGeometry``, ngmesh = geo.GenerateMesh(maxh=0.2) else: ngmesh = netgen.libngpy._meshing.Mesh(2) - mesh = Mesh(Mesh(ngmesh,comm=COMM_WORLD).curve_field(2)) + mesh = Mesh(ngmesh, comm=COMM_WORLD, netgen_flags={"degree": 2}) VTKFile("output/MeshExample7.pvd").write(mesh) .. figure:: Example7.png From e58dd9cf1f017961f8a13947f539f7436ae514ed Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Mar 2026 21:32:01 +0000 Subject: [PATCH 15/25] review comments --- demos/netgen/netgen_mesh.py.rst | 26 +++++++------- firedrake/mesh.py | 39 +++++++++++---------- firedrake/mg/netgen.py | 1 - firedrake/netgen.py | 61 ++++++++++++++++++++------------- 4 files changed, 70 insertions(+), 57 deletions(-) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 2f4439187c..ea2ee5a38d 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -112,7 +112,7 @@ In code this becomes: :: Now we are ready to assemble the stiffness matrix for the problem. Since we want to enforce Dirichlet boundary conditions we construct a ``DirichletBC`` object and we use the ``GetRegionNames`` method from the Netgen mesh in order to map the label we have given when describing the geometry to the PETSc ``DMPLEX`` IDs. In particular if we look for the IDs of boundary element labeled either "line" or "curve" we would get:: - labels = [i+1 for i, name in enumerate(ngmsh.GetRegionNames(codim=1)) if name in ["line","curve"]] + labels = [i+1 for i, name in enumerate(ngmsh.GetRegionNames(codim=1)) if name in ["line", "curve"]] bc = DirichletBC(V, 0, labels) print(labels) @@ -151,7 +151,7 @@ Then a SLEPc Eigenvalue Problem Solver (``EPS``) is initialised and set up to us u = TrialFunction(V) v = TestFunction(V) a = inner(grad(u), grad(v))*dx - m = (u*v)*dx + m = inner(u, v)*dx uh = Function(V) bc = DirichletBC(V, 0, labels) A = assemble(a, bcs=bc) @@ -217,7 +217,7 @@ In order to do so we begin by computing the value of the indicator using a piece sum_marked_eta += sum(eta[new_marked]) marked += new_marked frac -= delfrac - markedVec0.getArray()[:] = 1.0*marked[:] + markedVec0.getArray()[:] = marked[:] sct(markedVec0, markedVec, mode=PETSc.Scatter.Mode.REVERSE) return mark @@ -343,7 +343,7 @@ We can set the degree of the geometry via ``netgen_flags`` keyword argument of t for i in range(6): wp.Line(0.6).Arc(0.4, 60) shape = wp.Face() - ngmesh = OCCGeometry(shape,dim=2).GenerateMesh(maxh=1.) + ngmesh = OCCGeometry(shape, dim=2).GenerateMesh(maxh=1.) else: ngmesh = netgen.libngpy._meshing.Mesh(2) mesh = Mesh(ngmesh, comm=COMM_WORLD, netgen_flags={"degree": 3}) @@ -362,16 +362,16 @@ We will now show how to solve the Poisson problem on a high-order mesh, of order if COMM_WORLD.rank == 0: shape = Sphere(Pnt(0,0,0), 1) - ngmesh = OCCGeometry(shape,dim=3).GenerateMesh(maxh=1.) + ngmesh = OCCGeometry(shape, dim=3).GenerateMesh(maxh=1.) else: ngmesh = netgen.libngpy._meshing.Mesh(3) - mesh = Mesh(ngmesh, netgen_flags={"degree": 3}) + mesh = Mesh(ngmesh, netgen_flags={"degree": 4}) # Solving the Poisson problem VTKFile("output/MeshExample6.pvd").write(mesh) x, y, z = SpatialCoordinate(mesh) V = FunctionSpace(mesh, "CG", 3) - f = Function(V).interpolate(1+0*x) + f = Function(V).interpolate(Constant(1)) u = TrialFunction(V) v = TestFunction(V) a = inner(grad(u), grad(v)) * dx @@ -397,12 +397,12 @@ It is also possible to construct high-order meshes using the ``SplineGeometry``, from mpi4py import MPI if COMM_WORLD.rank == 0: - geo = CSG2d() - circle = Circle(center=(1,1), radius=0.1, bc="curve").Maxh(0.01) - rect = Rectangle(pmin=(0,1), pmax=(1,2), - bottom="b", left="l", top="t", right="r") - geo.Add(rect-circle) - ngmesh = geo.GenerateMesh(maxh=0.2) + geo = CSG2d() + circle = Circle(center=(1,1), radius=0.1, bc="curve").Maxh(0.01) + rect = Rectangle(pmin=(0,1), pmax=(1,2), + bottom="b", left="l", top="t", right="r") + geo.Add(rect-circle) + ngmesh = geo.GenerateMesh(maxh=0.2) else: ngmesh = netgen.libngpy._meshing.Mesh(2) mesh = Mesh(ngmesh, comm=COMM_WORLD, netgen_flags={"degree": 2}) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 4ef210c00e..a625d8592a 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2984,12 +2984,12 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): :arg order: the order of the curved mesh. :arg permutation_tol: tolerance used to construct the permutation of the reference element. :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. ''' - utils.check_netgen_installed() - from firedrake.netgen import find_permutation, plex_to_netgen_numbering, netgen_distribute + from firedrake.netgen import find_permutation, netgen_distribute from firedrake.functionspace import FunctionSpace from firedrake.function import Function @@ -3013,24 +3013,24 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): # Compute reference points using fiat fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent nodes = fiat_element.dual_basis() - ref = [] + ref_pts = [] for node in nodes: # Assert singleton point for each node. pt, = node.get_point_dict().keys() - ref.append(pt) - reference_space_points = np.array(ref) + ref_pts.append(pt) + reference_points = np.array(ref_pts) # Construct numpy arrays for physical domain data - physical_space_points = np.zeros( - (ng_dimension, reference_space_points.shape[0], self.geometric_dimension) + physical_points = np.zeros( + (ng_dimension, reference_points.shape[0], self.geometric_dimension) ) - curved_space_points = np.zeros( - (ng_dimension, reference_space_points.shape[0], self.geometric_dimension) + curved_points = np.zeros( + (ng_dimension, reference_points.shape[0], self.geometric_dimension) ) - self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) + self.netgen_mesh.CalcElementMapping(reference_points, physical_points) # NOTE: This will segfault for CSG! self.netgen_mesh.Curve(order) - self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) + self.netgen_mesh.CalcElementMapping(reference_points, curved_points) curved = ng_element.NumPy()["curved"] # Distribute curved cell data @@ -3040,18 +3040,19 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): own_curved = netgen_distribute(DG0, curved) own_curved = np.flatnonzero(own_curved[:num_cells]) - # Get numbering - cellNum = plex_to_netgen_numbering(self) - pyop2_index = cell_node_map.values[cellNum[own_curved]] - # Distribute coordinate data - own_curved_points = netgen_distribute(broken_space, curved_space_points)[own_curved] - own_physical_points = netgen_distribute(broken_space, physical_space_points)[own_curved] + own_curved_points = netgen_distribute(broken_space, curved_points)[own_curved] + own_physical_points = netgen_distribute(broken_space, physical_points)[own_curved] + + # Get broken indices + cstart, cend = self.topology_dm.getHeightStratum(0) + cellNum = np.array(list(map(self._cell_numbering.getOffset, range(cstart, cend)))) + broken_indices = cell_node_map.values[cellNum[own_curved]] # Find the correct coordinate permutation for each cell permutation = find_permutation( own_physical_points, - new_coordinates.dat.data_ro_with_halos[pyop2_index].real, + new_coordinates.dat.data_ro_with_halos[broken_indices].real, tol=permutation_tol, ) self.comm.Barrier() @@ -3060,7 +3061,7 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): own_curved_points[i] = own_curved_points[i, permutation[i]] # Assign the curved coordinates to the dat - new_coordinates.dat.data_wo_with_halos[pyop2_index] = own_curved_points + new_coordinates.dat.data_wo_with_halos[broken_indices] = own_curved_points return new_coordinates diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 595181efad..c3194c5212 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -280,7 +280,6 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): mesh.topology_dm.setRefineLevel(0) meshes.append(mesh) ngmesh = mesh.netgen_mesh - for l in range(levs): # Straighten the mesh ngmesh.Curve(1) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index ea7b3e0b68..75e09453ae 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -5,10 +5,11 @@ ''' import numpy as np import numpy.typing as npt -from petsc4py import PETSc from scipy.spatial.distance import cdist -import firedrake as fd +from pyop2.mpi import COMM_WORLD +from firedrake.petsc import PETSc +import firedrake # Netgen and ngsPETSc are not available when the documentation is getting built # because they do not have ARM wheels. @@ -29,44 +30,56 @@ class comp: Mesh = type(None) -def netgen_distribute(V, netgen_data): +def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, + netgen_data: npt.NDArray[np.inexact]): + """ + Distribute data from the netgen layout into the DMPlex layout. + + Parameters + ---------- + V + The target function space defining the data DMPlex layout. + netgen_data + The data in the layout of the underlying netgen mesh. + + Returns + ------- + npt.NDArray[np.inexact] + The data in the target DMPlex layout. + + """ netgen_data = np.asarray(netgen_data) mesh = V.mesh() sf = mesh.sfBC_orig - if sf is not None: + if sf is None: + # This mesh was not redistributed at construction. + # This means that the underlying netgen mesh represents + # the local part of the mesh owned by this process. + # Therefore the netgen data is already distributed. + plex_data = netgen_data + else: plex = mesh.topology_dm nshape = netgen_data.shape dtype = netgen_data.dtype - section = V.dm.getDefaultSection() - marked = V.dof_dset.layout_vec sfBCInv = sf.createInverse() - section0, marked0 = plex.distributeField(sfBCInv, section, marked) + section = V.dm.getDefaultSection() + vec = V.dof_dset.layout_vec + section0, vec0 = plex.distributeField(sfBCInv, section, vec) + vec0.set(0) plex_data = None - - marked0.set(0) - d = netgen_data for i in np.ndindex(V.shape): - di = d[(..., *i)].flatten() - marked0[:len(di)] = di - _, marked = plex.distributeField(sf, section0, marked0) - arr = marked.getArray() + di = netgen_data[(..., *i)].flatten() + vec0[:len(di)] = di + _, vec = plex.distributeField(sf, section0, vec0) + arr = vec.getArray() if plex_data is None: plex_data = np.empty(arr.shape + V.shape, dtype=dtype) plex_data[(..., *i)] = arr plex_data = plex_data.reshape(-1, *nshape[1:]) - else: - plex_data = netgen_data return plex_data -def plex_to_netgen_numbering(mesh): - plex = mesh.topology_dm - fstart, fend = plex.getHeightStratum(0) - cellNum = list(map(mesh._cell_numbering.getOffset, range(fstart, fend))) - return np.asarray(cellNum) - - @PETSc.Log.EventDecorator() def find_permutation(points_a: npt.NDArray[np.inexact], points_b: npt.NDArray[np.inexact], tol: float = 1e-5): @@ -132,7 +145,7 @@ class FiredrakeMesh: :param netgen_flags: The dictionary of flags to be passed to ngsPETSc. :arg comm: the MPI communicator. ''' - def __init__(self, mesh, netgen_flags=None, user_comm=fd.COMM_WORLD): + def __init__(self, mesh, netgen_flags=None, user_comm=COMM_WORLD): self.comm = user_comm # Parsing netgen flags if not isinstance(netgen_flags, dict): From 546f11a59bb64d76925d08bff7a3618402653be0 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 18 Mar 2026 21:46:54 +0000 Subject: [PATCH 16/25] fix --- demos/netgen/netgen_mesh.py.rst | 4 ++-- firedrake/netgen.py | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index ea2ee5a38d..1e6f80965e 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -346,7 +346,7 @@ We can set the degree of the geometry via ``netgen_flags`` keyword argument of t ngmesh = OCCGeometry(shape, dim=2).GenerateMesh(maxh=1.) else: ngmesh = netgen.libngpy._meshing.Mesh(2) - mesh = Mesh(ngmesh, comm=COMM_WORLD, netgen_flags={"degree": 3}) + mesh = Mesh(ngmesh, comm=COMM_WORLD, netgen_flags={"degree": 4}) VTKFile("output/MeshExample5.pvd").write(mesh) .. figure:: Example5.png @@ -365,7 +365,7 @@ We will now show how to solve the Poisson problem on a high-order mesh, of order ngmesh = OCCGeometry(shape, dim=3).GenerateMesh(maxh=1.) else: ngmesh = netgen.libngpy._meshing.Mesh(3) - mesh = Mesh(ngmesh, netgen_flags={"degree": 4}) + mesh = Mesh(ngmesh, netgen_flags={"degree": 3}) # Solving the Poisson problem VTKFile("output/MeshExample6.pvd").write(mesh) diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 75e09453ae..3de208463d 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -38,7 +38,7 @@ def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, Parameters ---------- V - The target function space defining the data DMPlex layout. + The target function space defining the DMPlex layout. netgen_data The data in the layout of the underlying netgen mesh. @@ -97,6 +97,7 @@ def find_permutation(points_a: npt.NDArray[np.inexact], points_b: npt.NDArray[np 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)] + if len(p) == 0: return p @@ -145,7 +146,7 @@ class FiredrakeMesh: :param netgen_flags: The dictionary of flags to be passed to ngsPETSc. :arg comm: the MPI communicator. ''' - def __init__(self, mesh, netgen_flags=None, user_comm=COMM_WORLD): + def __init__(self, mesh, netgen_flags, user_comm=COMM_WORLD): self.comm = user_comm # Parsing netgen flags if not isinstance(netgen_flags, dict): From 8d917e7cc48734d327bf971751d04b9d2042b939 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Mar 2026 14:31:44 +0000 Subject: [PATCH 17/25] Remove branching --- firedrake/mesh.py | 41 +++++++++++------------ tests/firedrake/regression/test_netgen.py | 30 +++++++++++------ 2 files changed, 40 insertions(+), 31 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index a625d8592a..96d1d12276 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2936,8 +2936,9 @@ def refine_marked_elements(self, mark, netgen_flags=None): with mark.dat.vec as mvec: if self.sfBC_orig is None: - perm = list(map(self._cell_numbering.getOffset, range(mvec.getSize()))) - mark_np = mvec.getArray()[perm] + cstart, cend = self.topology_dm.getHeightStratum(0) + cellNum = list(map(self._cell_numbering.getOffset, range(cstart, cend))) + mark_np = mvec.getArray()[cellNum] else: sfBCInv = self.sfBC_orig.createInverse() _, marked0 = self.topology_dm.distributeField(sfBCInv, @@ -2945,29 +2946,27 @@ def refine_marked_elements(self, mark, netgen_flags=None): mvec) mark_np = marked0.getArray() + max_refs = 0 if mark_np.size == 0 else int(mark_np.max()) netgen_mesh = self.netgen_mesh.Copy() refine_faces = netgen_flags.get("refine_faces", False) - if self.comm.rank == 0: - max_refs = int(mark_np.max()) - for r in range(max_refs): - cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() - cells.NumPy()["refine"] = (mark_np[:len(cells)] > 0) - if not refine_faces and dim == 3: - netgen_mesh.Elements2D().NumPy()["refine"] = False - netgen_mesh.Refine(adaptive=True) - mark_np -= 1 - if r < max_refs - 1: - parents = netgen_mesh.parentelements if dim == 3 else netgen_mesh.parentsurfaceelements - parents = parents.NumPy().astype(int).flatten() - cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() - num_coarse_cells = len(cells) - num_fine_cells = parents.shape[0] - indices = np.arange(num_fine_cells, dtype=int) - fine_cells = indices > num_coarse_cells + for r in range(max_refs): + cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() + cells.NumPy()["refine"] = (mark_np[:len(cells)] > 0) + if not refine_faces and dim == 3: + netgen_mesh.Elements2D().NumPy()["refine"] = False + netgen_mesh.Refine(adaptive=True) + mark_np -= 1 + if r < max_refs - 1: + parents = netgen_mesh.parentelements if dim == 3 else netgen_mesh.parentsurfaceelements + parents = parents.NumPy().astype(int).flatten() + num_fine_cells = parents.shape[0] + num_coarse_cells = mark_np.size + indices = np.arange(num_fine_cells, dtype=int) + while (indices >= num_coarse_cells).any(): + fine_cells = (indices >= num_coarse_cells) indices[fine_cells] = parents[indices[fine_cells]] - mark_np = mark_np[indices] + mark_np = mark_np[indices] - self.comm.Barrier() return Mesh(netgen_mesh, reorder=self._did_reordering, distribution_parameters=self._distribution_parameters, diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 21a054b643..31fbac30e7 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -4,24 +4,34 @@ @pytest.mark.skipnetgen +@pytest.mark.parallel([1, 3]) def test_create_netgen_mesh_high_order(): from netgen.geom2d import Circle, CSG2d geo = CSG2d() - - circle = Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle") - geo.Add(circle) - + geo.Add(Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle")) ngmesh = geo.GenerateMesh(maxh=0.75) # Test that setting the degree in netgen_flags produces a high-order mesh - mesh = Mesh(ngmesh, netgen_flags={"degree": 3}) - assert mesh.coordinates.function_space().ufl_element().degree() == 3 + order = 3 + mesh1 = Mesh(ngmesh, netgen_flags={"degree": order}) + assert mesh1.coordinates.function_space().ufl_element().degree() == order + dim = mesh1.topological_dimension + DG0 = FunctionSpace(mesh1, "DG", 0) + markers = Function(DG0) + + # Test mesh refinement: 1 refinement + markers.assign(1) + mesh2 = mesh1.refine_marked_elements(markers) + assert FunctionSpace(mesh1, "DG", 0).dim() * 2**dim == FunctionSpace(mesh2, "DG", 0).dim() + # Test that refining a high-order mesh gives a high-order mesh + assert mesh2.coordinates.function_space().ufl_element().degree() == order + # Test mesh refinement: 2 refinements + markers.assign(2) + mesh3 = mesh1.refine_marked_elements(markers) + assert FunctionSpace(mesh1, "DG", 0).dim() * 4**dim == FunctionSpace(mesh3, "DG", 0).dim() # Test that refining a high-order mesh gives a high-order mesh - DG0 = FunctionSpace(mesh, "DG", 0) - markers = Function(DG0).assign(1) - mesh2 = mesh.refine_marked_elements(markers) - assert mesh2.coordinates.function_space().ufl_element().degree() == 3 + assert mesh3.coordinates.function_space().ufl_element().degree() == order def square_geometry(h): From 942ae517078a433197017cf4458fa3576baeb18a Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Mar 2026 22:37:46 +0000 Subject: [PATCH 18/25] Update demos/netgen/netgen_mesh.py.rst --- demos/netgen/netgen_mesh.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 1e6f80965e..f0f4864a16 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -332,7 +332,7 @@ As usual, we generate a mesh for the described geometry and use the Firedrake-Ne High-order Meshes ------------------ It is possible to construct high-order meshes for a geometry constructed in Netgen. -We can set the degree of the geometry via ``netgen_flags`` keyword argument of the ``Mesh`` constructor. :: +We can set the degree of the geometry via the ``netgen_flags`` keyword argument of the ``Mesh`` constructor. :: from netgen.occ import WorkPlane, OCCGeometry import netgen From d8fc40b3ed74f9267b7e43f5132d1b3e860c6be6 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Mar 2026 22:37:56 +0000 Subject: [PATCH 19/25] Update demos/netgen/netgen_mesh.py.rst --- demos/netgen/netgen_mesh.py.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index f0f4864a16..6d0f488c83 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -371,7 +371,7 @@ We will now show how to solve the Poisson problem on a high-order mesh, of order VTKFile("output/MeshExample6.pvd").write(mesh) x, y, z = SpatialCoordinate(mesh) V = FunctionSpace(mesh, "CG", 3) - f = Function(V).interpolate(Constant(1)) + f = Constant(1) u = TrialFunction(V) v = TestFunction(V) a = inner(grad(u), grad(v)) * dx From 2551b9d93028cff07c2a8bc514c74a5d924e0154 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Mar 2026 22:51:05 +0000 Subject: [PATCH 20/25] Update firedrake/mesh.py --- firedrake/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 96d1d12276..8f7ee0b6bc 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -3352,7 +3352,7 @@ def Mesh(meshfile, **kwargs): degree = netgen_flags.get("degree", 1) if degree != 1: permutation_tol = netgen_flags.get("permutation_tol", 1e-8) - cg = netgen_flags.get("cg", False) + cg = netgen_flags.get("cg", None) coordinates = mesh.curve_field( order=degree, permutation_tol=permutation_tol, From a6d3171cc4b587e8b3960cace2feec27dc01a093 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 19 Mar 2026 23:17:00 +0000 Subject: [PATCH 21/25] Update tests/firedrake/multigrid/test_netgen_gmg.py --- tests/firedrake/multigrid/test_netgen_gmg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 8274f32012..79f8d50fc3 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -51,8 +51,8 @@ def test_netgen_mg(ngmesh, netgen_degree): labels = [i+1 for i, name in enumerate(ngmesh.GetRegionNames(codim=1)) if name in ["surface"]] x = SpatialCoordinate(mesh) - uexact = dot(x, x) - bcs = DirichletBC(V, uexact, labels) + uexact = 1-dot(x, x) + bcs = DirichletBC(V, 0, labels) L = a(uexact, v) uh = Function(V) From 9a4dc3abbd01c5bf427cb66c02555da92a611250 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 20 Mar 2026 11:57:50 +0000 Subject: [PATCH 22/25] Fix --- firedrake/mesh.py | 5 +++- firedrake/mg/netgen.py | 26 +++++++++----------- tests/firedrake/multigrid/test_netgen_gmg.py | 9 +++---- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 8f7ee0b6bc..1ddb40c4b1 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2992,6 +2992,9 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): from firedrake.functionspace import FunctionSpace from firedrake.function import Function + if not hasattr(self, "netgen_mesh"): + raise ValueError("Cannot curve a mesh that has not been generated by netgen.") + if cg_field is None: cg_field = not self.coordinates.function_space().finat_element.is_dg() @@ -3027,7 +3030,7 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): (ng_dimension, reference_points.shape[0], self.geometric_dimension) ) self.netgen_mesh.CalcElementMapping(reference_points, physical_points) - # NOTE: This will segfault for CSG! + # 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"] diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index c3194c5212..13646ae29b 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -139,17 +139,15 @@ def ref_grad(u): dmhooks.push_appctx(dm, ctx) solver.solve() if not cg: - element = ho.function_space().ufl_element().sub_elements[0].reconstruct(degree=degree) - space = fd.VectorFunctionSpace(linear, fd.BrokenElement(element)) - ho = fd.Function(space).interpolate(ho) + ho = fd.Function(ho.function_space().broken_space()).interpolate(ho) else: raise NotImplementedError("Snapping to Netgen meshes is only implemented for 2D meshes.") - return fd.Mesh(ho, comm=linear.comm, distribution_parameters=linear._distribution_parameters) + return reconstruct_mesh(linear, ho) def uniformRefinementRoutine(ngmesh, cdm): ''' - Routing called inside of NetgenHierarchy to compute refined ngmesh and plex. + Routine called inside of NetgenHierarchy to compute refined ngmesh and plex. ''' # We refine the DMPlex mesh uniformly logger.info(f"\t\t\t[{time.time()}]Refining the plex") @@ -235,7 +233,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): the coarse mesh, otherwise, these options override the default. """ - if mesh.geometric_dimension > 3: + if mesh.geometric_dimension not in {2, 3}: raise NotImplementedError("Netgen hierachies are only implemented for 2D and 3D meshes.") logger.info(f"Creating a Netgen hierarchy with {levs} levels.") comm = mesh.comm @@ -280,7 +278,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): mesh.topology_dm.setRefineLevel(0) meshes.append(mesh) ngmesh = mesh.netgen_mesh - for l in range(levs): + for l in range(1, levs+1): # Straighten the mesh ngmesh.Curve(1) rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm) @@ -292,7 +290,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): # Optimises the mesh, for example smoothing if ngmesh.dim == 2: ngmesh.OptimizeMesh2d(MeshingParameters(optimize2d=optMoves)) - elif mesh.dim == 3: + elif ngmesh.dim == 3: ngmesh.OptimizeVolumeMesh(MeshingParameters(optimize3d=optMoves)) else: raise ValueError("Only 2D and 3D meshes can be optimised.") @@ -316,24 +314,24 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): lgmaps.append((no, o)) # Curve the mesh - if order[l+1] != mesh.coordinates.function_space().ufl_element().degree(): + if order[l] != mesh.coordinates.function_space().ufl_element().degree(): logger.info("\t\t\tCurving the mesh ...") tic = time.time() if snap == "geometry": coordinates = mesh.curve_field( - order=order[l+1], + order=order[l], permutation_tol=permutation_tol, cg_field=cg, ) mesh = reconstruct_mesh(mesh, coordinates) elif snap == "coarse": - mesh = snapToCoarse(meshes[0].coordinates, mesh, order[l+1], snap_smoothing, cg) + mesh = snapToCoarse(meshes[0].coordinates, mesh, order[l], snap_smoothing, cg) toc = time.time() logger.info(f"\t\t\tMeshed curved. Time taken: {toc-tic}") - logger.info(f"\t\tLevel {l+1}: with {ngmesh.Coordinates().shape[0]}\ - vertices, with order {order[l+1]}, snapping to {snap}\ + logger.info(f"\t\tLevel {l}: with {ngmesh.Coordinates().shape[0]}\ + vertices, with order {order[l]}, snapping to {snap}\ and optimisation moves {optMoves}.") - mesh.topology_dm.setRefineLevel(l+1) + mesh.topology_dm.setRefineLevel(l) meshes.append(mesh) # Populate the coarse to fine map coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes, lgmaps) diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 79f8d50fc3..bb05426739 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -56,14 +56,12 @@ def test_netgen_mg(ngmesh, netgen_degree): L = a(uexact, v) uh = Function(V) - rtol = 1E-8 uerr = uexact - uh - err0 = assemble(a(uerr, uerr)) solve(a == L, uh, bcs=bcs, solver_parameters={ "ksp_type": "cg", "ksp_norm_type": "natural", "ksp_max_it": 14, - "ksp_rtol": rtol, + "ksp_rtol": 1E-8, "ksp_monitor": None, "pc_type": "mg", "mg_levels_pc_type": "python", @@ -72,5 +70,6 @@ def test_netgen_mg(ngmesh, netgen_degree): "mg_coarse_pc_type": "lu", "mg_coarse_pc_factor_mat_solver_type": "mumps", }) - errf = assemble(a(uerr, uerr)) - assert errf < err0 * (2*rtol) + err = assemble(a(uerr, uerr)) ** 0.5 + expected = 6E-2 if netgen_degree[-1] == 1 else 3E-4 + assert err < expected From b4dc8222055e538f4b4f8c5aecb724deb6bede4f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 20 Mar 2026 14:37:35 +0000 Subject: [PATCH 23/25] fix docs --- firedrake/mesh.py | 2 +- firedrake/netgen.py | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 1ddb40c4b1..bf853cc8d1 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2983,7 +2983,7 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): :arg order: the order of the curved mesh. :arg permutation_tol: tolerance used to construct the permutation of the reference element. :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. + Defaults to the continuity of the coordinates of the original mesh. ''' utils.check_netgen_installed() diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 3de208463d..57b442b83a 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -4,7 +4,6 @@ This file was copied from ngsPETSc. ''' import numpy as np -import numpy.typing as npt from scipy.spatial.distance import cdist from pyop2.mpi import COMM_WORLD @@ -31,7 +30,7 @@ class comp: def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, - netgen_data: npt.NDArray[np.inexact]): + netgen_data: np.ndarray): """ Distribute data from the netgen layout into the DMPlex layout. @@ -44,7 +43,7 @@ def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, Returns ------- - npt.NDArray[np.inexact] + ``np.ndarray`` The data in the target DMPlex layout. """ @@ -81,7 +80,7 @@ def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, @PETSc.Log.EventDecorator() -def find_permutation(points_a: npt.NDArray[np.inexact], points_b: npt.NDArray[np.inexact], +def find_permutation(points_a: np.ndarray, points_b: np.ndarray, tol: float = 1e-5): """ Find all permutations between a list of two sets of points. From b82acbada1412d23231eed75ae218a32e878adeb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 21 Mar 2026 13:41:03 +0000 Subject: [PATCH 24/25] fix --- firedrake/mesh.py | 29 ++++++++++++++--------------- firedrake/mg/netgen.py | 10 +++++++--- 2 files changed, 21 insertions(+), 18 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index bf853cc8d1..07998afe48 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2930,10 +2930,9 @@ def refine_marked_elements(self, mark, netgen_flags=None): raise ValueError("Adaptive refinement requires a netgen mesh.") if netgen_flags is None: netgen_flags = self.netgen_flags - dim = self.geometric_dimension - if dim not in {2, 3}: + gdim = self.geometric_dimension + if gdim not in {2, 3}: raise NotImplementedError("No implementation for dimension other than 2 and 3.") - with mark.dat.vec as mvec: if self.sfBC_orig is None: cstart, cend = self.topology_dm.getHeightStratum(0) @@ -2941,27 +2940,28 @@ def refine_marked_elements(self, mark, netgen_flags=None): mark_np = mvec.getArray()[cellNum] else: sfBCInv = self.sfBC_orig.createInverse() - _, marked0 = self.topology_dm.distributeField(sfBCInv, - self._cell_numbering, - mvec) - mark_np = marked0.getArray() - + _, mvec0 = self.topology_dm.distributeField(sfBCInv, + self._cell_numbering, + mvec) + mark_np = mvec0.getArray() max_refs = 0 if mark_np.size == 0 else int(mark_np.max()) + # Create a copy of the netgen mesh netgen_mesh = self.netgen_mesh.Copy() refine_faces = netgen_flags.get("refine_faces", False) for r in range(max_refs): - cells = netgen_mesh.Elements3D() if dim == 3 else netgen_mesh.Elements2D() + cells = netgen_mesh.Elements3D() if gdim == 3 else netgen_mesh.Elements2D() cells.NumPy()["refine"] = (mark_np[:len(cells)] > 0) - if not refine_faces and dim == 3: - netgen_mesh.Elements2D().NumPy()["refine"] = False + if gdim == 3: + faces = netgen_mesh.Elements2D() + faces.NumPy()["refine"] = refine_faces netgen_mesh.Refine(adaptive=True) mark_np -= 1 if r < max_refs - 1: - parents = netgen_mesh.parentelements if dim == 3 else netgen_mesh.parentsurfaceelements - parents = parents.NumPy().astype(int).flatten() + parents = netgen_mesh.parentelements if gdim == 3 else netgen_mesh.parentsurfaceelements + parents = parents.NumPy()["i"] num_fine_cells = parents.shape[0] num_coarse_cells = mark_np.size - indices = np.arange(num_fine_cells, dtype=int) + indices = np.arange(num_fine_cells, dtype=PETSc.IntType) while (indices >= num_coarse_cells).any(): fine_cells = (indices >= num_coarse_cells) indices[fine_cells] = parents[indices[fine_cells]] @@ -2987,7 +2987,6 @@ def curve_field(self, order, permutation_tol=1e-8, cg_field=None): ''' utils.check_netgen_installed() - from firedrake.netgen import find_permutation, netgen_distribute from firedrake.functionspace import FunctionSpace from firedrake.function import Function diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index 13646ae29b..c3f9979033 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -40,9 +40,13 @@ def snapToNetgenDMPlex(ngmesh, petscPlex, comm): This function snaps the coordinates of a DMPlex mesh to the coordinates of a Netgen mesh. ''' logger.info(f"\t\t\t[{time.time()}]Snapping the DMPlex to NETGEN mesh") - if len(ngmesh.Elements3D()) == 0: + + tdim = petscPlex.getDimension() + if tdim == 1: + ng_coelement = ngmesh.Elements0D() + elif tdim == 2: ng_coelement = ngmesh.Elements1D() - else: + elif tdim == 3: ng_coelement = ngmesh.Elements2D() # When we create a netgen mesh from a refined plex, @@ -57,7 +61,7 @@ def snapToNetgenDMPlex(ngmesh, petscPlex, comm): tic = time.time() ngCoordinates = ngmesh.Coordinates() petscCoordinates = petscPlex.getCoordinatesLocal().getArray() - petscCoordinates = petscCoordinates.reshape(-1, ngmesh.dim) + petscCoordinates = petscCoordinates.reshape(-1, petscPlex.getCoordinateDim()) petscCoordinates[nodes_to_correct_index] = ngCoordinates[nodes_to_correct_index] petscPlexCoordinates = petscPlex.getCoordinatesLocal() petscPlexCoordinates.setArray(petscCoordinates.reshape((-1, 1))) From 2b116d5b9dc213460157fd6a397605b4be241fb1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 21 Mar 2026 15:34:08 +0000 Subject: [PATCH 25/25] fixes --- firedrake/mg/netgen.py | 34 ++++++++++++++-------------------- 1 file changed, 14 insertions(+), 20 deletions(-) diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py index c3f9979033..aaa44ace85 100644 --- a/firedrake/mg/netgen.py +++ b/firedrake/mg/netgen.py @@ -35,27 +35,25 @@ def trim_util(T): return T -def snapToNetgenDMPlex(ngmesh, petscPlex, comm): +def snapToNetgenDMPlex(ngmesh, petscPlex): ''' This function snaps the coordinates of a DMPlex mesh to the coordinates of a Netgen mesh. ''' logger.info(f"\t\t\t[{time.time()}]Snapping the DMPlex to NETGEN mesh") - tdim = petscPlex.getDimension() - if tdim == 1: + gdim = petscPlex.getCoordinateDim() + if gdim == 1: ng_coelement = ngmesh.Elements0D() - elif tdim == 2: + elif gdim == 2: ng_coelement = ngmesh.Elements1D() - elif tdim == 3: + elif gdim == 3: ng_coelement = ngmesh.Elements2D() - # When we create a netgen mesh from a refined plex, # the netgen mesh represents the local submesh. # Therefore, there is no need to distribute the netgen data nodes_to_correct = ng_coelement.NumPy()["nodes"] nodes_to_correct = trim_util(nodes_to_correct) - nodes_to_correct_sorted = np.hstack(nodes_to_correct.reshape((-1, 1))) - nodes_to_correct_sorted.sort() + nodes_to_correct_sorted = nodes_to_correct.flatten() nodes_to_correct_index = np.unique(nodes_to_correct_sorted) logger.info(f"\t\t\t[{time.time()}]Nodes have been corrected") tic = time.time() @@ -64,7 +62,7 @@ def snapToNetgenDMPlex(ngmesh, petscPlex, comm): petscCoordinates = petscCoordinates.reshape(-1, petscPlex.getCoordinateDim()) petscCoordinates[nodes_to_correct_index] = ngCoordinates[nodes_to_correct_index] petscPlexCoordinates = petscPlex.getCoordinatesLocal() - petscPlexCoordinates.setArray(petscCoordinates.reshape((-1, 1))) + petscPlexCoordinates.setArray(petscCoordinates.flatten()) petscPlex.setCoordinatesLocal(petscPlexCoordinates) toc = time.time() logger.info(f"\t\t\tSnap the DMPlex to NETGEN mesh. Time taken: {toc - tic} seconds") @@ -237,10 +235,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): the coarse mesh, otherwise, these options override the default. """ - if mesh.geometric_dimension not in {2, 3}: - raise NotImplementedError("Netgen hierachies are only implemented for 2D and 3D meshes.") - logger.info(f"Creating a Netgen hierarchy with {levs} levels.") - comm = mesh.comm + tdim = mesh.topological_dimension # Parse netgen flags if not isinstance(flags, dict): flags = mesh.netgen_flags @@ -254,6 +249,7 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): snap_smoothing = flags.get("snap_smoothing", "hyperelastic") cg = flags.get("cg", not mesh.coordinates.function_space().finat_element.is_dg()) nested = flags.get("nested", snap in ["coarse"]) + logger.info(f"Creating a Netgen hierarchy with {levs} levels.") logger.info(f"\tOrder of the hierarchy: {order}") logger.info(f"\tRefinement type: {refType}") logger.info(f"\tSnap to {snap} using {snap_smoothing} smoothing (if snapping to coarse)") @@ -268,10 +264,8 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): cg_field=cg, ) mesh = reconstruct_mesh(mesh, coordinates) - # Make a plex (cdm) without overlap. dm_cell_type, = mesh.dm_cell_types - tdim = mesh.topology_dm.getDimension() cdm = dmcommon.submesh_create(mesh.topology_dm, tdim, "celltype", dm_cell_type, True) cdm.removeLabel("pyop2_core") cdm.removeLabel("pyop2_owned") @@ -287,17 +281,17 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): ngmesh.Curve(1) rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm) cdm = rdm - # Snap the mesh to the Netgen mesh - if snap == "geometry": - snapToNetgenDMPlex(ngmesh, rdm, comm) if optMoves: # Optimises the mesh, for example smoothing - if ngmesh.dim == 2: + if tdim == 2: ngmesh.OptimizeMesh2d(MeshingParameters(optimize2d=optMoves)) - elif ngmesh.dim == 3: + elif tdim == 3: ngmesh.OptimizeVolumeMesh(MeshingParameters(optimize3d=optMoves)) else: raise ValueError("Only 2D and 3D meshes can be optimised.") + # Snap the mesh to the Netgen mesh + if snap == "geometry": + snapToNetgenDMPlex(ngmesh, rdm) # We construct a Firedrake mesh from the DMPlex mesh parameters = {}