diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 480220cbd0..6d0f488c83 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 @@ -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 the ``netgen_flags`` keyword argument of the ``Mesh`` constructor. :: from netgen.occ import WorkPlane, OCCGeometry import netgen @@ -344,10 +343,10 @@ In particular, we need to pass the degree of the polynomial field we want to use 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(Mesh(ngmesh, comm=COMM_WORLD).curve_field(4)) + mesh = Mesh(ngmesh, comm=COMM_WORLD, netgen_flags={"degree": 4}) VTKFile("output/MeshExample5.pvd").write(mesh) .. figure:: Example5.png @@ -363,15 +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(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) V = FunctionSpace(mesh, "CG", 3) - f = Function(V).interpolate(1+0*x) + f = Constant(1) u = TrialFunction(V) v = TestFunction(V) a = inner(grad(u), grad(v)) * dx @@ -397,15 +397,15 @@ 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(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 diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 59531066c1..07998afe48 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 @@ -2922,47 +2924,57 @@ 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 not hasattr(self, "netgen_mesh"): + raise ValueError("Adaptive refinement requires a netgen mesh.") if netgen_flags is None: - netgen_flags = {} - DistParams = self._distribution_parameters - els = {2: self.netgen_mesh.Elements2D, 3: self.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 - if self.sfBC is not None: - sfBCInv = self.sfBC.createInverse() - getIdx = lambda x: x - _, marked0 = self.topology_dm.distributeField(sfBCInv, - self._cell_numbering, - marked) - if self.comm.Get_rank() == 0: - mark = marked0.getArray() - max_refs = np.max(mark) - for _ in range(int(max_refs)): - for i, el in enumerate(els[dim]()): - if mark[getIdx(i)] > 0: - el.refine = True - else: - el.refine = False - if not refine_faces and dim == 3: - self.netgen_mesh.Elements2D().NumPy()["refine"] = 0 - self.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.libngpy._meshing.Mesh(dim), - distribution_parameters=DistParams, comm=self.comm) - else: + netgen_flags = self.netgen_flags + 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) + cellNum = list(map(self._cell_numbering.getOffset, range(cstart, cend))) + mark_np = mvec.getArray()[cellNum] + else: + sfBCInv = self.sfBC_orig.createInverse() + _, 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 gdim == 3 else netgen_mesh.Elements2D() + cells.NumPy()["refine"] = (mark_np[:len(cells)] > 0) + 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 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=PETSc.IntType) + while (indices >= num_coarse_cells).any(): + 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, + 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): + 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 @@ -2970,115 +2982,87 @@ 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. + :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. ''' - import firedrake as fd - utils.check_netgen_installed() + from firedrake.netgen import find_permutation, netgen_distribute + 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.") - from firedrake.netgen import find_permutation + 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 len(self.netgen_mesh.Elements3D()) == 0: - ng_element = self.netgen_mesh.Elements2D + if self.topological_dimension == 2: + ng_element = self.netgen_mesh.Elements2D() else: - ng_element = self.netgen_mesh.Elements3D - ng_dimension = len(ng_element()) - geom_dim = self.geometric_dimension + ng_element = self.netgen_mesh.Elements3D() + ng_dimension = len(ng_element) - # Construct the mesh as a Firedrake function - if cg_field: - firedrake_space = fd.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)) + # Construct the coordinates as a Firedrake function + coords_space = self.coordinates.function_space().reconstruct(degree=order) + broken_space = coords_space.broken_space() + if not cg_field: + coords_space = broken_space + 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) - 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) - ) - 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"] - # 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) - ) + ref_pts = [] + for node in nodes: + # Assert singleton point for each node. + pt, = node.get_point_dict().keys() + ref_pts.append(pt) + reference_points = np.array(ref_pts) + + # Construct numpy arrays for physical domain data + physical_points = np.zeros( + (ng_dimension, reference_points.shape[0], self.geometric_dimension) + ) + curved_points = np.zeros( + (ng_dimension, reference_points.shape[0], self.geometric_dimension) + ) + self.netgen_mesh.CalcElementMapping(reference_points, physical_points) + # NOTE: This will segfault for MeshHierarchy on a netgen CSG geometry + self.netgen_mesh.Curve(order) + self.netgen_mesh.CalcElementMapping(reference_points, curved_points) + curved = ng_element.NumPy()["curved"] - # Broadcast curved cell point data - self.comm.Bcast(physical_space_points, root=0) - self.comm.Bcast(curved_space_points, root=0) + # Distribute curved cell data 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]) - # 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)] + # Distribute coordinate data + own_curved_points = netgen_distribute(broken_space, curved_points)[own_curved] + own_physical_points = netgen_distribute(broken_space, physical_points)[own_curved] - # 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] - - # 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]) + # 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( - 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_ro_with_halos[broken_indices].real, + tol=permutation_tol, ) - + self.comm.Barrier() # Apply the permutation to each cell in turn - for ii, p in enumerate(curved_space_points): - curved_space_points[ii] = p[permutation[ii]] + 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] = curved_space_points.reshape(-1, geom_dim) - + new_coordinates.dat.data_wo_with_halos[broken_indices] = own_curved_points return new_coordinates @@ -3364,6 +3348,30 @@ 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", None) + coordinates = mesh.curve_field( + order=degree, + permutation_tol=permutation_tol, + cg_field=cg, + ) + # Do not redistribute the mesh + reorder_noop = None + temp = Mesh(coordinates, + reorder=reorder_noop, + perm_is=mesh._dm_renumbering, + distribution_parameters=DISTRIBUTION_PARAMETERS_NOOP, + comm=mesh.comm) + temp.netgen_mesh = mesh.netgen_mesh + temp.netgen_flags = mesh.netgen_flags + temp._distribution_parameters = mesh._distribution_parameters + temp._did_reordering = mesh._did_reordering + mesh = temp mesh.submesh_parent = submesh_parent mesh._tolerance = tolerance 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 82acab25aa..aaa44ace85 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 @@ -34,33 +35,34 @@ 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") - if len(ngmesh.Elements3D()) == 0: - ng_coelement = ngmesh.Elements1D - else: - ng_coelement = ngmesh.Elements2D - if comm.rank == 0: - nodes_to_correct = ng_coelement().NumPy()["nodes"] - 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") + + gdim = petscPlex.getCoordinateDim() + if gdim == 1: + ng_coelement = ngmesh.Elements0D() + elif gdim == 2: + ng_coelement = ngmesh.Elements1D() + 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() 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))) + 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") @@ -139,17 +141,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,50 +235,37 @@ def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): the coarse mesh, otherwise, these options override the default. """ - if mesh.geometric_dimension > 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 = {} + 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) - 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) 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"]) + 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)") # Firedrake quantities meshes = [] lgmaps = [] - parameters = {} - if distribution_parameters is not None: - parameters.update(distribution_parameters) - else: - parameters.update(mesh._distribution_parameters) - parameters["partition"] = False # Curve the mesh - if order[0] > 1: - ho_field = mesh.curve_field( + if order[0] != mesh.coordinates.function_space().ufl_element().degree(): + coordinates = mesh.curve_field( order=order[0], permutation_tol=permutation_tol, - cg_field=cg + cg_field=cg, ) - temp = fd.Mesh(ho_field, distribution_parameters=parameters, comm=comm) - temp.netgen_mesh = mesh.netgen_mesh - temp._tolerance = mesh.tolerance - mesh = temp + 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") @@ -289,49 +276,79 @@ 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) 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: + if tdim == 2: ngmesh.OptimizeMesh2d(MeshingParameters(optimize2d=optMoves)) - elif mesh.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 = {} + 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, + reorder=False, + distribution_parameters=parameters, + tolerance=mesh.tolerance) mesh.netgen_mesh = ngmesh + mesh.netgen_flags = flags + + no = impl.create_lgmap(rdm) + o = impl.create_lgmap(mesh.topology_dm) + lgmaps.append((no, o)) + # Curve the mesh - if order[l+1] > 1: + if order[l] != mesh.coordinates.function_space().ufl_element().degree(): logger.info("\t\t\tCurving the mesh ...") tic = time.time() if snap == "geometry": - mesh = fd.Mesh( - mesh.curve_field(order=order[l+1], - location_tol=location_tol, - permutation_tol=permutation_tol), - distribution_parameters=parameters, - comm=comm) + coordinates = mesh.curve_field( + order=order[l], + permutation_tol=permutation_tol, + cg_field=cg, + ) + mesh = reconstruct_mesh(mesh, coordinates) elif snap == "coarse": - mesh = snapToCoarse(ho_field, 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(1 + l) + 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) return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, 1, nested=nested) + + +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 + tmesh.sfBC_orig = mesh.sfBC_orig + return tmesh diff --git a/firedrake/netgen.py b/firedrake/netgen.py index 1033b2f627..57b442b83a 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -4,11 +4,11 @@ This file was copied from ngsPETSc. ''' 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,8 +29,58 @@ class comp: Mesh = type(None) +def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, + netgen_data: np.ndarray): + """ + Distribute data from the netgen layout into the DMPlex layout. + + Parameters + ---------- + V + The target function space defining the DMPlex layout. + netgen_data + The data in the layout of the underlying netgen mesh. + + Returns + ------- + ``np.ndarray`` + The data in the target DMPlex layout. + + """ + netgen_data = np.asarray(netgen_data) + mesh = V.mesh() + sf = mesh.sfBC_orig + 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 + + sfBCInv = sf.createInverse() + section = V.dm.getDefaultSection() + vec = V.dof_dset.layout_vec + section0, vec0 = plex.distributeField(sfBCInv, section, vec) + vec0.set(0) + plex_data = None + for i in np.ndindex(V.shape): + 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:]) + return plex_data + + @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. @@ -46,6 +96,10 @@ 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: @@ -91,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, user_comm=fd.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): diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 2b4dcc6602..bb05426739 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -1,105 +1,75 @@ import pytest - 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=[3, 2]) +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 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 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(): - ngmesh = create_netgen_mesh_circle() - mesh = Mesh(ngmesh) - 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) - - -@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) - nh = MeshHierarchy(mesh, 2, netgen_flags={"degree": [1, 2, 3]}) - mesh = nh[-1] +@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)} + 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(mh) + + 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) 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(): - ngmesh = create_netgen_mesh_circle() - mesh = Mesh(ngmesh) - 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 = 1-dot(x, x) + bcs = DirichletBC(V, 0, labels) + L = a(uexact, v) + uh = Function(V) + + uerr = uexact - uh + solve(a == L, uh, bcs=bcs, solver_parameters={ + "ksp_type": "cg", + "ksp_norm_type": "natural", + "ksp_max_it": 14, + "ksp_rtol": 1E-8, + "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", + }) + err = assemble(a(uerr, uerr)) ** 0.5 + expected = 6E-2 if netgen_degree[-1] == 1 else 3E-4 + assert err < expected diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index dca547a9f0..31fbac30e7 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -3,6 +3,37 @@ import pytest +@pytest.mark.skipnetgen +@pytest.mark.parallel([1, 3]) +def test_create_netgen_mesh_high_order(): + from netgen.geom2d import Circle, CSG2d + geo = CSG2d() + 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 + 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 + assert mesh3.coordinates.function_space().ufl_element().degree() == order + + def square_geometry(h): from netgen.geom2d import SplineGeometry geo = SplineGeometry() @@ -190,44 +221,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)) @@ -236,6 +243,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 @@ -305,81 +313,33 @@ def adapt(mesh, eta): (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.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 + assert error_estimators[-1] < 0.06 - 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] +def test_netgen_manifold(): + from netgen.meshing import MeshingStep + 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 + b = 1.6 - theta = 0.5 - should_refine = conditional(gt(eta, theta*eta_max), 1, 0) - markers.interpolate(should_refine) + def Curve(t): + return Pnt(0, R+a*np.cos(t), b*np.sin(t)) - refined_mesh = mesh.refine_marked_elements(markers) - return refined_mesh + n = 100 + pnts = [Curve(2*np.pi*t/n) for t in range(n+1)] - 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 + spline = SplineApproximation(pnts) + f = Face(Wire(spline)) - geo = OCCGeometry(L, dim=2) - ngmsh = geo.GenerateMesh(maxh=0.1) - mesh = Mesh(ngmsh) + 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) - 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()) - mesh = adapt(mesh, eta) - assert error_estimators[-1] < 0.06 + assert mesh.topological_dimension == 2 + assert mesh.geometric_dimension == 3