diff --git a/.github/actions/install/action.yml b/.github/actions/install/action.yml index 54c1411d7f..97a50c5d26 100644 --- a/.github/actions/install/action.yml +++ b/.github/actions/install/action.yml @@ -144,7 +144,7 @@ runs: : # because they rely on non-PyPI versions of petsc4py. pip install --no-build-isolation --no-deps \ "$PETSC_DIR"/"$PETSC_ARCH"/externalpackages/git.slepc/src/binding/slepc4py - pip install --no-deps git+https://github.com/NGSolve/ngsPETSc.git netgen-mesher netgen-occt + pip install --no-deps git+https://github.com/NGSolve/ngsPETSc.git@pbrubeck/netgen-plex netgen-mesher netgen-occt : # We have to pass '--no-build-isolation' to use a custom petsc4py EXTRA_PIP_FLAGS='--no-build-isolation' diff --git a/firedrake/adapt.py b/firedrake/adapt.py index 071ca8231a..767dd4cbed 100644 --- a/firedrake/adapt.py +++ b/firedrake/adapt.py @@ -8,7 +8,6 @@ from firedrake.function import Function from firedrake.functionspace import FunctionSpace from firedrake.mesh import Mesh, DISTRIBUTION_PARAMETERS_NOOP -from firedrake.netgen import _transfer_high_order_coordinates # PETSc's DMAdaptFlag value requesting refinement, for the adapt label. @@ -63,10 +62,8 @@ def _copy_adaptive_refinement_metadata(source_mesh, target_mesh): target_mesh._distribution_parameters = dict(source_mesh._distribution_parameters) target_mesh._did_reordering = source_mesh._did_reordering target_mesh._tolerance = source_mesh.tolerance - if hasattr(source_mesh, "netgen_mesh") and not hasattr(target_mesh, "netgen_mesh"): - target_mesh.netgen_mesh = source_mesh.netgen_mesh - if hasattr(source_mesh, "netgen_flags") and not hasattr(target_mesh, "netgen_flags"): - target_mesh.netgen_flags = source_mesh.netgen_flags + if target_mesh._geometry_source is None: + target_mesh._geometry_source = source_mesh._geometry_source def refine_marked_elements(mesh, cell_marker): @@ -106,6 +103,8 @@ def refine_marked_elements(mesh, cell_marker): try: for ref in range(num_refinements): new_dm = _adapt_marked_cells(current_mesh, current_mark) + if mesh._geometry_source is not None: + mesh._geometry_source.snap(new_dm) current_mesh = Mesh( new_dm, dim=mesh.geometric_dimension, @@ -131,10 +130,9 @@ def refine_marked_elements(mesh, cell_marker): coarse_dm.removeLabel(PARENT_LABEL) final_mesh = current_mesh - if hasattr(mesh, "netgen_mesh"): + if mesh._geometry_source is not None: order = mesh.coordinates.function_space().ufl_element().degree() - if order > 1: - final_mesh = _transfer_high_order_coordinates(mesh, final_mesh, order) + final_mesh = mesh._geometry_source.recurve(final_mesh, order) final_mesh.topology_dm.removeLabel(PARENT_LABEL) final_mesh.adaptive_parent = mesh diff --git a/firedrake/cython/dmcommon.pyx b/firedrake/cython/dmcommon.pyx index 54774b456f..3d5c72a039 100644 --- a/firedrake/cython/dmcommon.pyx +++ b/firedrake/cython/dmcommon.pyx @@ -2108,6 +2108,144 @@ def reordered_coords(PETSc.DM dm, PETSc.Section global_numbering, shape, referen return coords +@cython.boundscheck(False) +@cython.wraparound(False) +def set_cell_coordinates(PETSc.DM dm, + np.ndarray[PetscScalar, ndim=3, mode="c"] values): + """Set coordinate closures from cellwise values. + + Parameters + ---------- + dm : PETSc.DMPlex + The DMPlex whose coordinate vector is populated. + values : numpy.ndarray + Array of shape ``(num_cells, num_nodes, coordinate_dim)`` in + the coordinate ``PetscFE`` closure ordering. + """ + cdef: + PETSc.Section section = dm.getCoordinateSection() + PETSc.Vec coordinates = dm.getCoordinatesLocal() + PetscInt cStart, cEnd, c + PetscInt closure_size + PetscScalar *closure = NULL + PetscInt expected_closure_size = values.shape[1] * values.shape[2] + + get_height_stratum(dm.dm, 0, &cStart, &cEnd) + if values.shape[0] != cEnd - cStart: + raise ValueError( + f"Expected coordinate data for {cEnd - cStart} cells, " + f"got {values.shape[0]}" + ) + if cStart < cEnd: + CHKERR(DMPlexVecGetClosure( + dm.dm, + section.sec, + coordinates.vec, + cStart, + &closure_size, + &closure, + )) + CHKERR(DMPlexVecRestoreClosure( + dm.dm, + section.sec, + coordinates.vec, + cStart, + &closure_size, + &closure, + )) + if closure_size != expected_closure_size: + raise ValueError( + f"Coordinate closure has size {closure_size}, " + f"expected {expected_closure_size}" + ) + for c in range(cStart, cEnd): + CHKERR(DMPlexVecSetClosure( + dm.dm, + section.sec, + coordinates.vec, + c, + &values[c - cStart, 0, 0], + PETSC_INSERT_VALUES, + )) + dm.setCoordinatesLocal(coordinates) + + +@cython.boundscheck(False) +@cython.wraparound(False) +def reordered_coords_high_order(PETSc.DM dm, + PETSc.Section firedrake_section, + shape): + """Return high-order DMPlex coordinates in a Firedrake layout. + + The DMPlex coordinate discretization and Firedrake coordinate element + must assign the same number and ordering of nodes to each topological + entity. + + Parameters + ---------- + dm : PETSc.DMPlex + The DMPlex containing high-order coordinates. + firedrake_section : PETSc.Section + Scalar section of the matching Firedrake coordinate space. + shape : tuple + Output shape ``(num_coordinate_nodes, coordinate_dim)``. + """ + cdef: + PETSc.Section coordinate_section = dm.getCoordinateSection() + PETSc.Vec coordinate_vector = dm.getCoordinatesLocal() + const PetscScalar *dm_coordinates + PetscInt pStart, pEnd, qStart, qEnd, p + PetscInt firedrake_dof, coordinate_dof + PetscInt firedrake_offset, coordinate_offset + PetscInt i, j, gdim = shape[1], total_dof = 0 + np.ndarray coords = np.empty(shape, dtype=ScalarType) + + pStart, pEnd = firedrake_section.getChart() + qStart, qEnd = coordinate_section.getChart() + if (pStart, pEnd) != (qStart, qEnd): + raise ValueError( + "DMPlex and Firedrake coordinate sections have different charts: " + f"{(qStart, qEnd)} != {(pStart, pEnd)}" + ) + + CHKERR(VecGetArrayRead(coordinate_vector.vec, &dm_coordinates)) + try: + for p in range(pStart, pEnd): + CHKERR(PetscSectionGetDof( + firedrake_section.sec, p, &firedrake_dof + )) + CHKERR(PetscSectionGetDof( + coordinate_section.sec, p, &coordinate_dof + )) + if coordinate_dof != gdim * firedrake_dof: + raise ValueError( + f"Coordinate sections disagree at DMPlex point {p}: " + f"{coordinate_dof} != {gdim} * {firedrake_dof}" + ) + if firedrake_dof == 0: + continue + CHKERR(PetscSectionGetOffset( + firedrake_section.sec, p, &firedrake_offset + )) + CHKERR(PetscSectionGetOffset( + coordinate_section.sec, p, &coordinate_offset + )) + for i in range(firedrake_dof): + for j in range(gdim): + coords[firedrake_offset + i, j] = \ + dm_coordinates[coordinate_offset + gdim * i + j] + total_dof += firedrake_dof + finally: + CHKERR(VecRestoreArrayRead(coordinate_vector.vec, &dm_coordinates)) + + if total_dof != shape[0]: + raise ValueError( + f"Coordinate section contains {total_dof} nodes, " + f"expected {shape[0]}" + ) + return coords + + def _get_expanded_dm_dg_coords(dm: PETSc.DM, ndofs: np.ndarray): """Return the DM DG coordinates expanded to the full closure size. diff --git a/firedrake/cython/petschdr.pxi b/firedrake/cython/petschdr.pxi index 42ac97e24d..12a402c168 100644 --- a/firedrake/cython/petschdr.pxi +++ b/firedrake/cython/petschdr.pxi @@ -25,6 +25,8 @@ cdef extern from "petsc.h": ctypedef enum PetscErrorCode: PETSC_SUCCESS PETSC_ERR_LIB + ctypedef enum PetscInsertMode "InsertMode": + PETSC_INSERT_VALUES "INSERT_VALUES" cdef extern from "petscsys.h" nogil: PetscErrorCode PetscMalloc1(PetscInt,void*) @@ -82,6 +84,9 @@ cdef extern from "petscdmplex.h" nogil: PetscErrorCode DMPlexSetCellType(PETSc.PetscDM,PetscInt,PetscDMPolytopeType) PetscErrorCode DMPlexGetCellType(PETSc.PetscDM,PetscInt,PetscDMPolytopeType*) + PetscErrorCode DMPlexVecGetClosure(PETSc.PetscDM,PETSc.PetscSection,PETSc.PetscVec,PetscInt,PetscInt*,PetscScalar**) + PetscErrorCode DMPlexVecRestoreClosure(PETSc.PetscDM,PETSc.PetscSection,PETSc.PetscVec,PetscInt,PetscInt*,PetscScalar**) + PetscErrorCode DMPlexVecSetClosure(PETSc.PetscDM,PETSc.PetscSection,PETSc.PetscVec,PetscInt,PetscScalar[],PetscInsertMode) cdef extern from "petscdmlabel.h" nogil: struct _n_DMLabel diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 881cdfe523..9fc81c825f 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -544,7 +544,6 @@ 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) @@ -1147,7 +1146,6 @@ 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 @@ -2395,6 +2393,7 @@ def __init__(self, coordinates): self.extruded = isinstance(topology, ExtrudedMeshTopology) self.variable_layers = self.extruded and topology.variable_layers self._base_mesh = None # this is set by extruded meshes in a later step + self._geometry_source = None # these are set by firedrake.adapt.refine_marked_elements self.adaptive_parent = None self.adaptive_cell_maps = None @@ -2977,104 +2976,43 @@ def refine_marked_elements(self, mark): return refine_marked_elements(self, mark) @PETSc.Log.EventDecorator() - def curve_field(self, order, permutation_tol=None, cg_field=None): - '''Return a function containing the curved coordinates of the mesh. + def curve_field(self, + order: int, + permutation_tol: float = 1.e-8, + cg_field: bool | None = None) -> object: + """Return a function containing re-evaluated mesh coordinates. - This method requires that the mesh has been constructed from a - netgen mesh. + This method requires a mesh construction source which supports + coordinate re-evaluation, such as Netgen. - :arg order: the order of the curved mesh. - :arg permutation_tol: ignored. - :arg cg_field: return a CG function field representing the mesh, as opposed to a DG field. - Defaults to the continuity of the coordinates of the original mesh. - - ''' - 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.") - if permutation_tol is not None: - warnings.warn( - "permutation_tol is no longer required to obtain the curved coordinates. " - "This kwarg will be removed in a future release.", - FutureWarning, - ) + Parameters + ---------- + order + Polynomial degree of the re-evaluated coordinates. + permutation_tol + Retained for compatibility. DMPlex coordinate elements do not use + geometric permutation searches. + cg_field + Whether to return continuous coordinates. Defaults to the + continuity of the current coordinates. - if cg_field is None: - cg_field = not self.coordinates.function_space().finat_element.is_dg() + Returns + ------- + Function + The re-evaluated coordinate field. - # Check if the mesh is a surface mesh or two dimensional mesh - if self.topological_dimension == 2: - ng_element = self.netgen_mesh.Elements2D() - else: - ng_element = self.netgen_mesh.Elements3D() - ng_dimension = len(ng_element) - - # 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 - nodes = fiat_element.dual_basis() - ref_pts = [] - entity_ids = fiat_element.entity_dofs() - for dim in sorted(entity_ids): - for entity in sorted(entity_ids[dim]): - for i in entity_ids[dim][entity]: - # Assert singleton point for each node. - pt, = nodes[i].get_point_dict().keys() - ref_pts.append(pt) - reference_points = np.array(ref_pts) - - # Construct numpy arrays for physical domain data - 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.Curve(1) - self.netgen_mesh.CalcElementMapping(reference_points, physical_points) - self.netgen_mesh.Curve(order) - self.netgen_mesh.CalcElementMapping(reference_points, curved_points) - curved = ng_element.NumPy()["curved"] - - # 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]) - - # 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] - - # 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[broken_indices].real, + """ + if self._geometry_source is None: + raise ValueError( + "This mesh has no geometry source capable of re-evaluating " + "its coordinates." + ) + return self._geometry_source.curve_field( + self, + order, + permutation_tol=permutation_tol, + cg_field=cg_field, ) - 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_wo_with_halos[broken_indices] = own_curved_points - return new_coordinates @PETSc.Log.EventDecorator() @@ -3117,6 +3055,8 @@ def make_mesh_from_coordinates(coordinates, name, tolerance=0.5): mesh._tolerance = tolerance mesh._did_reordering = orig_mesh._did_reordering mesh._distribution_parameters = orig_mesh._distribution_parameters + if isinstance(orig_mesh, MeshGeometry): + mesh._geometry_source = orig_mesh._geometry_source return mesh @@ -3375,7 +3315,11 @@ def Mesh(meshfile, **kwargs): utils._init() - from_netgen = netgen and isinstance(meshfile, netgen.libngpy._meshing.Mesh) + from_netgen = user_comm.bcast( + bool(netgen and isinstance(meshfile, netgen.libngpy._meshing.Mesh)) + if user_comm.rank == 0 else None, + root=0, + ) # We don't need to worry about using a user comm in these cases as # they all immediately call a petsc4py which in turn uses a PETSc @@ -3386,12 +3330,13 @@ def Mesh(meshfile, **kwargs): if MPI.Comm.Compare(user_comm, plex.comm.tompi4py()) not in {MPI.CONGRUENT, MPI.IDENT}: raise ValueError("Communicator used to create `plex` must be at least congruent to the communicator used to create the mesh") elif from_netgen: - from firedrake.netgen import FiredrakeMesh + from firedrake.netgen import NetgenGeometry petsctools.cite("Betteridge2024") - netgen_flags = kwargs.get("netgen_flags", {"quad": False, "transform": None, "purify_to_tets": False}) - netgen_firedrake_mesh = FiredrakeMesh(meshfile, netgen_flags, user_comm) - plex = netgen_firedrake_mesh.meshMap.petscPlex + netgen_flags = kwargs.get("netgen_flags", {}) + geometry_source = NetgenGeometry(meshfile, netgen_flags, user_comm) + netgen_flags = geometry_source.options + plex = geometry_source.plex plex.setName(_generate_default_mesh_topology_name(name)) else: @@ -3422,36 +3367,17 @@ def Mesh(meshfile, **kwargs): permutation_name=kwargs.get("permutation_name"), submesh_parent=submesh_parent.topology if submesh_parent else None, comm=user_comm) - mesh = make_mesh_from_mesh_topology(topology, name) + if from_netgen and netgen_flags.get("degree", 1) != 1: + from firedrake.netgen import _mesh_from_coordinate_dm + + mesh = _mesh_from_coordinate_dm( + topology, name, netgen_flags["degree"] + ) + else: + mesh = make_mesh_from_mesh_topology(topology, name) 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", None) - 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.sfBC = mesh.sfBC - temp.sfBC_orig = mesh.sfBC_orig - temp._distribution_parameters = mesh._distribution_parameters - temp._did_reordering = mesh._did_reordering - mesh = temp + mesh._geometry_source = geometry_source mesh.submesh_parent = submesh_parent mesh._tolerance = tolerance diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py index 8c1a6c237a..8523478f75 100644 --- a/firedrake/mg/mesh.py +++ b/firedrake/mg/mesh.py @@ -201,9 +201,11 @@ def MeshHierarchy(mesh, refinement_levels=0, the number of refinements for each level in the hierarchy. Adaptive refinement only supports one refinement per level. netgen_flags : bool, dict - either a bool or a dictionary containing options for Netgen. - If not False the hierachy is constructed using ngsPETSc, if - None hierarchy constructed in a standard manner. + options for the Netgen geometry attached to ``mesh``. The only + option read here is ``degree``: either an integer, or a sequence + with one entry per mesh in the hierarchy, giving the polynomial + degree of the coordinates of each level. It defaults to the degree + of the coordinates of ``mesh``. distribution_parameters : dict options controlling mesh distribution, see :py:func:`.Mesh` for details. If ``None``, use the same distribution @@ -229,17 +231,37 @@ def MeshHierarchy(mesh, refinement_levels=0, HierarchyBase The mesh hierarchy. - """ + Notes + ----- + If ``mesh`` carries a geometry source, such as the Netgen geometry of a + mesh built from a Netgen mesh, then it is attached to every level, and + the coordinates of each level are re-evaluated on that geometry. - if (isinstance(netgen_flags, bool) and netgen_flags) or isinstance(netgen_flags, dict): + """ + if netgen_flags is True or isinstance(netgen_flags, dict): utils.check_netgen_installed() - from firedrake.mg.netgen import NetgenHierarchy + if mesh._geometry_source is None: + raise RuntimeError("Cannot pass netgen_flags to a mesh that has not " + "been generated by Netgen.") + else: + netgen_flags = {} + + # The coordinate degree of each mesh in the hierarchy. + nmeshes = refinement_levels*refinements_per_level + 1 + coarse_degree = mesh.coordinates.function_space().ufl_element().degree() + degree = netgen_flags.get("degree", coarse_degree) + if isinstance(degree, int): + degree = (degree,)*nmeshes + degree = tuple(degree) + if len(degree) != nmeshes: + raise ValueError(f"Expecting one coordinate degree per mesh in the " + f"hierarchy ({nmeshes}), got {len(degree)}") + + geometry_source = mesh._geometry_source + if geometry_source is not None: + if degree[0] != coarse_degree: + mesh = firedrake.Mesh(mesh.curve_field(degree[0])) - if hasattr(mesh, "netgen_mesh"): - return NetgenHierarchy(mesh, refinement_levels, flags=netgen_flags) - else: - raise RuntimeError("Cannot create a NetgenHierarchy from a mesh that has not been generated by\ - Netgen.") if callbacks is not None: before, after = callbacks else: @@ -251,6 +273,11 @@ def MeshHierarchy(mesh, refinement_levels=0, cdm = mesh.topology_dm if refinement_levels > 0: cdm = make_unoverlapped_dm(cdm) + if geometry_source is not None and degree[0] > 1: + # Refinement, and conversion back to Netgen, both require the + # straight-sided cells of a linear coordinate discretization. + from firedrake.netgen import _linearize_coordinate_dm + _linearize_coordinate_dm(cdm) cdm.setRefinementUniform(True) dms = [cdm] for i in range(refinement_levels*refinements_per_level): @@ -259,6 +286,8 @@ def MeshHierarchy(mesh, refinement_levels=0, rdm = cdm.refine() if i % refinements_per_level == 0: after(rdm, i) + if geometry_source is not None: + geometry_source.snap(rdm) # Fix up coords if refining embedded circle or sphere if hasattr(mesh, '_radius'): # FIXME, really we need some CAD-like representation @@ -281,7 +310,7 @@ def MeshHierarchy(mesh, refinement_levels=0, parameters["partition"] = False meshes = [mesh] - for rdm in dms[1:]: + for rdm, order in zip(dms[1:], degree[1:]): fmesh = mesh_builder( rdm, dim=mesh.geometric_dimension, @@ -289,6 +318,8 @@ def MeshHierarchy(mesh, refinement_levels=0, reorder=reorder, comm=mesh.comm, ) + if geometry_source is not None: + fmesh = geometry_source.recurve(fmesh, order) meshes.append(fmesh) # Build local-to-global maps and coarse/fine cell maps between diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py deleted file mode 100644 index a957e50c7d..0000000000 --- a/firedrake/mg/netgen.py +++ /dev/null @@ -1,356 +0,0 @@ -""" -This file was copied from ngsPETSc. -""" -import time -from fractions import Fraction - -import numpy as np -import ufl -from packaging.version import Version -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 - -# Netgen and ngsPETSc are not available when the documentation is getting built -# because they do not have ARM wheels. -try: - from netgen.meshing import MeshingParameters - from ngsPETSc.plex import MeshMapping -except ImportError: - pass - - -def trim_util(T): - """ - Trim zeros from a connectivity array T. - """ - if Version(np.__version__) >= Version("2.2"): - T = np.trim_zeros(T, "b", axis=1).astype(PETSc.IntType) - 1 - else: - T = (np.array([list(np.trim_zeros(a, "b")) for a in list(T)], dtype=PETSc.IntType) - 1) - return T - - -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") - - 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 = 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, petscPlex.getCoordinateDim()) - petscCoordinates[nodes_to_correct_index] = ngCoordinates[nodes_to_correct_index] - petscPlexCoordinates = petscPlex.getCoordinatesLocal() - 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") - - -def snapToCoarse(coarse, linear, degree, snap_smoothing, cg): - ''' - This function snaps the coordinates of a DMPlex mesh to the coordinates of a Netgen mesh. - ''' - dim = linear.geometric_dimension - if dim == 2: - space = fd.VectorFunctionSpace(linear, "CG", degree) - ho = fd.assemble(fd.interpolate(coarse, space)) - if snap_smoothing == "hyperelastic": - # Hyperelastic Smoothing - bcs = [fd.DirichletBC(space, ho, "on_boundary")] - quad_degree = 2*(degree+1)-1 - d = linear.topological_dimension - Q = fd.TensorFunctionSpace(linear, "DG", degree=0) - Jinv = ufl.JacobianInverse(linear) - hinv = fd.Function(Q) - hinv.interpolate(Jinv) - G = ufl.Jacobian(linear) * hinv - ijac = 1/abs(ufl.det(G)) - - def ref_grad(u): - return ufl.dot(ufl.grad(u), G) - - params = { - "snes_type": "newtonls", - "snes_linesearch_type": "l2", - "snes_max_it": 50, - "snes_rtol": 1E-8, - "snes_atol": 1E-8, - "snes_ksp_ew": True, - "snes_ksp_ew_rtol0": 1E-2, - "snes_ksp_ew_rtol_max": 1E-2, - } - params["mat_type"] = "aij" - coarse = { - "ksp_type": "preonly", - "pc_type": "lu", - "pc_mat_factor_type": "mumps", - } - gmg = { - "pc_type": "mg", - "mg_coarse": coarse, - "mg_levels": { - "ksp_max_it": 2, - "ksp_type": "chebyshev", - "pc_type": "jacobi", - }, - } - l = fd.mg.utils.get_level(linear)[1] - pc = gmg if l else coarse - params.update(pc) - ksp = { - "ksp_rtol": 1E-8, - "ksp_atol": 0, - "ksp_type": "minres", - "ksp_norm_type": "preconditioned", - } - params.update(ksp) - u = ho - F = ref_grad(u) - J = ufl.det(F) - psi = (1/2) * (ufl.inner(F, F)-d - ufl.ln(J**2)) - U = (psi * ijac)*fd.dx(degree=quad_degree) - dU = ufl.derivative(U, u, fd.TestFunction(space)) - problem = fd.NonlinearVariationalProblem(dU, u, bcs) - solver = fd.NonlinearVariationalSolver(problem, solver_parameters=params) - solver.set_transfer_manager(None) - ctx = solver._ctx - for c in problem.F.coefficients(): - dm = c.function_space().dm - dmhooks.push_appctx(dm, ctx) - solver.solve() - if not cg: - ho = fd.Function(ho.function_space().broken_space()).interpolate(ho) - else: - raise NotImplementedError("Snapping to Netgen meshes is only implemented for 2D meshes.") - return reconstruct_mesh(linear, ho) - - -def uniformRefinementRoutine(ngmesh, cdm): - ''' - 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") - cdm.setRefinementUniform(True) - rdm = cdm.refine() - rdm.removeLabel("pyop2_core") - rdm.removeLabel("pyop2_owned") - rdm.removeLabel("pyop2_ghost") - logger.info(f"\t\t\t[{time.time()}]Mapping the mesh to Netgen mesh") - tic = time.time() - mapping = MeshMapping(rdm, geo=ngmesh) - toc = time.time() - logger.info(f"\t\t\t[{time.time()}]Mapped the mesh to Netgen. Time taken: {toc-tic}") - return (rdm, mapping.ngMesh) - - -def uniformMapRoutine(meshes, lgmaps): - ''' - This function computes the coarse to fine and fine to coarse maps - for a uniform mesh hierarchy. - ''' - refinements_per_level = 1 - coarse_to_fine_cells = [] - fine_to_coarse_cells = [None] - for (coarse, fine), (clgmaps, flgmaps) in zip( - zip(meshes[:-1], meshes[1:]), - zip(lgmaps[:-1], lgmaps[1:]) - ): - c2f, f2c = impl.coarse_to_fine_cells(coarse, fine, clgmaps, flgmaps) - coarse_to_fine_cells.append(c2f) - fine_to_coarse_cells.append(f2c) - - coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), c2f) - for i, c2f in enumerate(coarse_to_fine_cells)) - fine_to_coarse_cells = dict((Fraction(i, refinements_per_level), f2c) - for i, f2c in enumerate(fine_to_coarse_cells)) - return (coarse_to_fine_cells, fine_to_coarse_cells) - - -def alfeldRefinementRoutine(ngmesh, cdm): - ''' - Routing called inside of NetgenHierarchy to compute refined ngmesh and plex. - ''' - # We refine the netgen mesh alfeld - ngmesh.SplitAlfeld() - # We refine the DMPlex mesh alfeld - tr = PETSc.DMPlexTransform().create(comm=PETSc.COMM_WORLD) - tr.setType(PETSc.DMPlexTransformType.REFINEREGULAR) - tr.setDM(cdm) - tr.setUp() - rdm = tr.apply(cdm) - return (rdm, ngmesh) - - -def alfeldMapRoutine(meshes): - ''' - This function computes the coarse to fine and fine to coarse maps - for a alfeld mesh hierarchy. - ''' - raise NotImplementedError("Alfeld refinement is not implemented yet.") - - -refinementTypes = {"uniform": (uniformRefinementRoutine, uniformMapRoutine), - "Alfeld": (alfeldRefinementRoutine, alfeldMapRoutine)} - - -def NetgenHierarchy(mesh, levs, flags, distribution_parameters=None): - """Create a Firedrake mesh hierarchy from Netgen/NGSolve meshes. - - :arg mesh: the Netgen/NGSolve mesh - :arg levs: the number of levels in the hierarchy - :arg flags: either a bool or a dictionary containing options for Netgen. - If not False the hierachy is constructed using ngsPETSc, if None hierarchy - constructed in a standard manner. Netgen flags includes: - - - degree, either an integer denoting the degree of curvature of all levels of - the mesh or a list of levs+1 integers denoting the degree of curvature of - each level of the mesh. - - tol, geometric tolerance adopted in snapToNetgenDMPlex. - - refinement_type, the refinment type to be used: uniform (default), Alfeld - :kwarg distribution_parameters: a dict of options controlling mesh distribution. - If ``None``, use the same distribution parameters as were used to distribute - the coarse mesh, otherwise, these options override the default. - - """ - tdim = mesh.topological_dimension - # Parse netgen flags - if not isinstance(flags, dict): - flags = mesh.netgen_flags - order = flags.get("degree", 1) - if isinstance(order, int): - order = [order]*(levs+1) - permutation_tol = flags.get("permutation_tol", None) - refType = flags.get("refinement_type", "uniform") - optMoves = flags.get("optimisation_moves", False) - snap = flags.get("snap_to", "geometry") - 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)") - # Firedrake quantities - meshes = [] - lgmaps = [] - # Curve the mesh - if order[0] != mesh.coordinates.function_space().ufl_element().degree(): - coordinates = mesh.curve_field( - order=order[0], - permutation_tol=permutation_tol, - cg_field=cg, - ) - mesh = reconstruct_mesh(mesh, coordinates) - # Make a plex (cdm) without overlap. - cdm = dmcommon.submesh_create(mesh.topology_dm, tdim, "depth", tdim, True) - cdm.removeLabel("pyop2_core") - cdm.removeLabel("pyop2_owned") - cdm.removeLabel("pyop2_ghost") - no = impl.create_lgmap(cdm) - o = impl.create_lgmap(mesh.topology_dm) - lgmaps.append((no, o)) - mesh.topology_dm.setRefineLevel(0) - meshes.append(mesh) - base_ngmesh = mesh.netgen_mesh - comm = mesh.comm - for l in range(1, levs+1): - rdm, ngmesh = refinementTypes[refType][0](base_ngmesh, cdm) - # `fd.Mesh` mutates `rdm` in place (e.g. adding overlap), so clone - # it first to keep an unoverlapped dm for the next refinement. - cdm = rdm.clone() - if optMoves: - # Optimises the mesh, for example smoothing - if tdim == 2: - ngmesh.OptimizeMesh2d(MeshingParameters(optimize2d=optMoves)) - 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, - comm=comm) - 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] != 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], - permutation_tol=permutation_tol, - cg_field=cg, - ) - mesh = reconstruct_mesh(mesh, coordinates) - elif snap == "coarse": - 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}: with {ngmesh.Coordinates().shape[0]}\ - vertices, with order {order[l]}, snapping to {snap}\ - and optimisation moves {optMoves}.") - 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 = mesh.sfBC - tmesh.sfBC_orig = mesh.sfBC_orig - return tmesh diff --git a/firedrake/netgen.py b/firedrake/netgen.py index c0c0369779..85d65cfe06 100644 --- a/firedrake/netgen.py +++ b/firedrake/netgen.py @@ -1,13 +1,12 @@ -''' -This module contains all the functions related to wrapping NGSolve meshes to Firedrake +"""Conversion between Netgen meshes and Firedrake meshes.""" +import pickle +from functools import cached_property -This file was copied from ngsPETSc. -''' import numpy as np -from scipy.spatial.distance import cdist -from pyop2.mpi import COMM_WORLD +from pyop2.mpi import COMM_WORLD, MPI from firedrake.petsc import PETSc +from firedrake.utils import ScalarType import firedrake # Netgen and ngsPETSc are not available when the documentation is getting built @@ -26,107 +25,245 @@ class ngs: "dummy class" class comp: "dummy class" - Mesh = type(None) + Mesh = type("_MissingNGSolveMesh", (), {}) -def netgen_distribute(V: firedrake.functionspaceimpl.WithGeometryBase, - netgen_data: np.ndarray): +def _coordinate_finite_element(plex: PETSc.DMPlex, + degree: int, + is_simplex: bool) -> PETSc.FE: + """Create the PETSc finite element used for DMPlex coordinates. + + Parameters + ---------- + plex + The DMPlex whose dimension and coordinate dimension determine the + element's reference cell and value shape. + degree + Polynomial degree of the coordinate element. + is_simplex + Whether ``plex`` is a simplex mesh. + + Returns + ------- + PETSc.FE + A Lagrange finite element of the requested degree. """ - Distribute data from the netgen layout into the DMPlex layout. + prefix = "netgen_coordinate_" + key = f"{prefix}petscspace_degree" + options = PETSc.Options() + previous = options.getString(key) if options.hasName(key) else None + options[key] = degree + try: + element = PETSc.FE().createDefault( + plex.getDimension(), + plex.getCoordinateDim(), + is_simplex, + prefix=prefix, + comm=plex.comm, + ) + finally: + if previous is None: + del options[key] + else: + options[key] = previous + return element + + +def _coordinate_reference_points(element: PETSc.FE, + tdim: int, + gdim: int) -> np.ndarray: + """Return PETSc coordinate nodes in Netgen reference coordinates.""" + dual = element.getDualSpace() + dimension = dual.getDimension() + if dimension % gdim: + raise ValueError( + f"Coordinate dual dimension {dimension} is not divisible by " + f"the geometric dimension {gdim}" + ) + + points = [] + for node in range(0, dimension, gdim): + point, _ = dual.getFunctional(node).getData() + if point.size != tdim: + raise ValueError("Netgen coordinates require point-evaluation nodes.") + for component in range(1, gdim): + component_point, _ = dual.getFunctional(node + component).getData() + if not np.array_equal(point, component_point): + raise ValueError( + "Coordinate components use different dual evaluation points." + ) + points.append(point) + + # Express the dual nodes in barycentric coordinates relative to PETSc's + # oriented reference-cell closure. In particular, PETSc's tetrahedron + # does not order the coordinate axes in the same way as its triangle. + reference_dm = dual.getDM() + reference_cell = reference_dm.getHeightStratum(0)[0] + vertices = reference_dm.getVecClosure( + reference_dm.getCoordinateSection(), + reference_dm.getCoordinatesLocal(), + reference_cell, + ).reshape(-1, tdim) + points = np.asarray(points) + coordinates = np.linalg.solve( + (vertices[1:] - vertices[0]).T, + (points - vertices[0]).T, + ).T + barycentric = np.column_stack((1 - coordinates.sum(axis=1), coordinates)) + + # Netgen's element transformation associates its reference vertices + # with the cyclically shifted connectivity order (last, first, ..., + # penultimate). + netgen_vertex_permutation = np.roll(np.arange(tdim + 1), 1) + return barycentric[:, netgen_vertex_permutation][:, 1:] + + +def _is_simplex(plex: PETSc.DMPlex) -> bool: + """Return whether every cell of a DMPlex, on every rank, is a simplex. Parameters ---------- - V - The target function space defining the DMPlex layout. - netgen_data - The data in the layout of the underlying netgen mesh. + plex + The DMPlex to test. Returns ------- - ``np.ndarray`` - The data in the target DMPlex layout. + bool + Whether the mesh is a simplex mesh. Ranks holding no cells abstain. + """ + cell_start, cell_end = plex.getHeightStratum(0) + return plex.comm.tompi4py().allreduce( + cell_start == cell_end or plex.isSimplex(), + op=MPI.LAND, + ) + +def _set_netgen_coordinates(plex: PETSc.DMPlex, + ngmesh: object, + degree: int, + *, + root_only: bool = True, + reset: bool = False) -> None: + """Attach Netgen-evaluated high-order coordinates to a DMPlex. + + Parameters + ---------- + plex + The DMPlex to curve. Its coordinate discretization is replaced by a + Lagrange element of the requested degree. + ngmesh + The Netgen mesh whose element mapping evaluates the coordinates. + degree + Polynomial degree of the new coordinates. + root_only + Whether only rank 0 holds a nonempty ``ngmesh`` and ``plex``, as is + the case before the mesh is distributed. Other ranks then contribute + no cell coordinates. + reset + Whether to first reset ``plex`` to a linear coordinate + discretization, needed when it already carries a stale + higher-order one. """ - 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 + from firedrake.cython import dmcommon + + if not isinstance(degree, int) or degree < 1: + raise ValueError("The Netgen coordinate degree must be a positive integer.") + comm = plex.comm.tompi4py() + is_simplex = _is_simplex(plex) + if not is_simplex: + raise NotImplementedError( + "High-order Netgen coordinates currently require a simplex mesh." + ) + + tdim = plex.getDimension() + gdim = plex.getCoordinateDim() + if reset: + # Refined DMs can carry the parent coordinate PetscFE alongside a + # linear coordinate section. Restore a consistent linear coordinate + # DM before installing the requested element. + linear_element = _coordinate_finite_element(plex, 1, is_simplex) + plex.setCoordinateDisc(linear_element, False, False) + element = _coordinate_finite_element(plex, degree, is_simplex) + points = _coordinate_reference_points(element, tdim, gdim) + plex.setCoordinateDisc(element, False, True) + + if not root_only or comm.rank == 0: + elements = { + 1: ngmesh.Elements1D, + 2: ngmesh.Elements2D, + 3: ngmesh.Elements3D, + }[tdim]() + values = np.empty((len(elements), len(points), gdim), dtype=np.float64) + ngmesh.Curve(degree) + ngmesh.CalcElementMapping(points, values) + values = np.asarray(values, dtype=ScalarType) 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: np.ndarray, points_b: np.ndarray): - """ Find all permutations between a list of two sets of points. - - Given two numpy arrays of shape (ncells, npoints, dim) containing - floating point coordinates for each cell, determine each index - permutation that takes `points_a` to `points_b`. Ie: - ``` - permutation = find_permutation(points_a, points_b) - assert np.allclose(points_a[permutation], points_b, rtol=0, atol=tol) - ``` - """ - if points_a.shape != points_b.shape: - raise ValueError("`points_a` and `points_b` must have the same shape.") - - # Match reference points instead of physical points to ensure scale invariance - dim = points_a.shape[-1] - vids = list(range(dim+1)) - # Infer the affine mapping (A, b) from the image of the vertices (first dim+1 dofs) - bs = points_a[:, vids[:1], :] - As = points_a[:, vids[1:], :] - As -= bs - Ainvs = np.linalg.inv(As) - # x_phys = A * x_ref + b <==> x_ref = inv(A) * (x_phys - b) - # Multiply inv(A) from the right, since the data is row-major - ref_points_a = np.matmul(points_a - bs, Ainvs) - ref_points_b = np.matmul(points_b - bs, Ainvs) - - p = [np.argmin(cdist(a, b), axis=0) for a, b in zip(ref_points_a, ref_points_b)] - - if len(p) == 0: - return p + values = np.empty((0, len(points), gdim), dtype=ScalarType) + dmcommon.set_cell_coordinates(plex, values) - try: - permutation = np.array(p, ndmin=2) - except ValueError as e: - raise ValueError( - "It was not possible to find a permutation for every cell" - " within the provided tolerance" - ) from e - if permutation.shape != points_a.shape[0:2]: - raise ValueError( - "It was not possible to find a permutation for every cell" - " within the provided tolerance" +def _linearize_coordinate_dm(plex: PETSc.DMPlex) -> None: + """Replace a simplex coordinate discretization by its linear interpolant. + + Parameters + ---------- + plex + The DMPlex whose coordinate discretization is replaced in place by + a degree-1 Lagrange element. + """ + is_simplex = _is_simplex(plex) + if not is_simplex: + raise NotImplementedError( + "Linearizing Netgen coordinates currently requires a simplex mesh." ) + element = _coordinate_finite_element(plex, 1, is_simplex) + plex.setCoordinateDisc(element, False, True) - return permutation + +def _mesh_from_coordinate_dm(topology: object, + name: str, + degree: int) -> object: + """Construct a Firedrake mesh from a high-order DMPlex coordinate field. + + Parameters + ---------- + topology + The mesh topology whose DMPlex already carries a degree-``degree`` + coordinate discretization. + name + Name of the returned mesh. + degree + Polynomial degree of the coordinate field to read from ``topology``. + + Returns + ------- + firedrake.mesh.MeshGeometry + A mesh built from the DMPlex's coordinates, reordered into + Firedrake's coordinate function space layout. + """ + import finat.ufl + from firedrake.function import CoordinatelessFunction + from firedrake.functionspace import FunctionSpace + from firedrake.cython import dmcommon + from firedrake.mesh import make_mesh_from_coordinates + + dm = topology.topology_dm + element = finat.ufl.VectorElement( + "Lagrange", + topology.ufl_cell(), + degree, + dim=dm.getCoordinateDim(), + ) + function_space = FunctionSpace(topology, element) + section = function_space.dm.getDefaultSection() + shape = (section.getStorageSize(), dm.getCoordinateDim()) + values = dmcommon.reordered_coords_high_order(dm, section, shape) + coordinates = CoordinatelessFunction( + function_space, + val=values, + name=f"{name}_coordinates", + ) + return make_mesh_from_coordinates(coordinates, name) def _transfer_high_order_coordinates(coarse_mesh, fine_mesh, order): @@ -185,45 +322,201 @@ def splitToQuads(plex, dim, comm): "Powell-Sabin": lambda x: x.SplitPowellSabin()} -class FiredrakeMesh: - ''' - This class creates a Firedrake mesh from Netgen/NGSolve meshes. +class NetgenGeometry: + """A Netgen geometry source associated with a Firedrake mesh. - :arg mesh: the mesh object, it can be either a Netgen/NGSolve mesh or a PETSc DMPlex - :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=COMM_WORLD): + Parameters + ---------- + mesh + A Netgen or NGSolve mesh. + options + Netgen construction and coordinate options. + user_comm + The communicator on which to build the DMPlex. + plex + An existing DMPlex represented by ``mesh``. If provided, topology + conversion is skipped. + """ + + def __init__(self, + mesh: object, + options: dict | None, + user_comm=COMM_WORLD, + plex: PETSc.DMPlex | None = None) -> None: self.comm = user_comm - # Parsing netgen flags - if not isinstance(netgen_flags, dict): - netgen_flags = {} - split2tets = netgen_flags.get("split_to_tets", False) - split = netgen_flags.get("split", False) - quad = netgen_flags.get("quad", False) - optMoves = netgen_flags.get("optimisation_moves", False) - # Checking the mesh format - if isinstance(mesh, (ngs.comp.Mesh, ngm.Mesh)): - if split2tets: - mesh = mesh.Split2Tets() - if split: - # Split mesh this includes Alfeld and Powell-Sabin - splitTypes[split](mesh) - if optMoves: - # Optimises the mesh, for example smoothing - if mesh.dim == 2: - mesh.OptimizeMesh2d(MeshingParameters(optimize2d=optMoves)) - elif mesh.dim == 3: - mesh.OptimizeVolumeMesh(MeshingParameters(optimize3d=optMoves)) - else: - raise ValueError("Only 2D and 3D meshes can be optimised.") - # We create the plex from the netgen mesh - self.meshMap = MeshMapping(mesh, comm=self.comm) - # We apply the DMPLEX transform + self.options = dict(options) if isinstance(options, dict) else {} + self.mesh = mesh + self.mesh_mapping = None + self.plex = plex + self._mesh_is_replicated = plex is not None + + if plex is not None: + return + self._mesh_is_replicated = self.comm.allreduce( + isinstance(mesh, (ngs.comp.Mesh, ngm.Mesh)), + op=MPI.LAND, + ) + if isinstance(mesh, ngs.comp.Mesh): + mesh = mesh.ngmesh + is_netgen = self.comm.bcast( + isinstance(mesh, ngm.Mesh) if self.comm.rank == 0 else None, + root=0, + ) + if is_netgen: + split2tets = self.options.get("split_to_tets", False) + split = self.options.get("split", False) + quad = self.options.get("quad", False) + opt_moves = self.options.get("optimisation_moves", False) + degree = self.options.get("degree", 1) + if split2tets or split or quad or opt_moves: + self._mesh_is_replicated = False + if self.comm.rank == 0: + if split2tets: + mesh = mesh.Split2Tets() + if split: + # Split mesh this includes Alfeld and Powell-Sabin + splitTypes[split](mesh) + if opt_moves: + # Optimises the mesh, for example smoothing + if mesh.dim == 2: + mesh.OptimizeMesh2d(MeshingParameters(optimize2d=opt_moves)) + elif mesh.dim == 3: + mesh.OptimizeVolumeMesh(MeshingParameters(optimize3d=opt_moves)) + else: + raise ValueError("Only 2D and 3D meshes can be optimised.") + self.mesh = mesh + self.mesh_mapping = MeshMapping(mesh, comm=self.comm) + self.plex = self.mesh_mapping.petscPlex + self.mesh = self.mesh_mapping.ngMesh if quad: - newplex = splitToQuads(self.meshMap.petscPlex, mesh.dim, comm=self.comm) - self.meshMap = MeshMapping(newplex) - elif isinstance(mesh, PETSc.DMPlex): - self.meshMap = MeshMapping(mesh) + dim = self.comm.bcast(mesh.dim if self.comm.rank == 0 else None, root=0) + self.plex = splitToQuads(self.plex, dim, comm=self.comm) + self.mesh_mapping = MeshMapping(self.plex) + self.mesh = self.mesh_mapping.ngMesh + if degree > 1: + # The linear coordinates are already those of the Netgen mesh. + _set_netgen_coordinates(self.plex, self.mesh, degree) else: raise ValueError("Mesh format not recognised.") + + @cached_property + def _local_mesh(self) -> object: + """Return a rank-local copy of the Netgen mesh when it is needed.""" + if self._mesh_is_replicated: + return self.mesh + serialized_mesh = self.comm.bcast( + pickle.dumps(self.mesh) if self.comm.rank == 0 else None, + root=0, + ) + return self.mesh if self.comm.rank == 0 else pickle.loads(serialized_mesh) + + def curve_field(self, + mesh: object, + order: int, + **kwargs: object) -> object: + """Return re-evaluated Netgen coordinates. + + Parameters + ---------- + mesh + The Firedrake mesh whose topology and construction options are + reused. + order + Polynomial degree of the new coordinates. + **kwargs + Compatibility options accepted by + :meth:`firedrake.mesh.MeshGeometry.curve_field`. + + Returns + ------- + firedrake.Function + The coordinate field of a rebuilt mesh. + """ + if kwargs.get("cg_field", True) is False: + raise NotImplementedError( + "cg_field=False is not supported: DMPlex high-order Netgen " + "coordinates are always continuous." + ) + options = { + key: value for key, value in self.options.items() + if key not in { + "split_to_tets", "split", "quad", "optimisation_moves", + "degree", "permutation_tol", "cg", + } + } + options["degree"] = order + rebuilt = firedrake.Mesh( + self.mesh, + name=mesh.name, + comm=mesh.comm, + reorder=mesh._did_reordering, + distribution_parameters=mesh._distribution_parameters, + tolerance=mesh.tolerance, + netgen_flags=options, + ) + return rebuilt.coordinates + + def snap(self, plex: PETSc.DMPlex) -> None: + """Snap the vertices of a derived DMPlex onto the geometry. + + Converting a DMPlex back to Netgen projects its boundary onto the + geometry, so the points of the resulting Netgen mesh are the snapped + vertices, in the order they were read out of ``plex``. A refined + DMPlex must be snapped before it is refined again, so that each level + is a subdivision of the level below it. + + Parameters + ---------- + plex + A DMPlex derived from this geometry, with linear coordinates, + which are snapped in place. + """ + if not _is_simplex(plex): + # ngsPETSc only converts simplex DMPlexes back to Netgen. + return + ngmesh = createNetgenMesh(plex, self._local_mesh) + coordinates = plex.getCoordinatesLocal() + coordinates.array[:] = ngmesh.Coordinates().reshape(-1) + plex.setCoordinatesLocal(coordinates) + + def recurve(self, fine_mesh: object, order: int) -> object: + """Re-evaluate coordinates on a derived mesh. + + The derived mesh is converted back to Netgen and its coordinates are + evaluated from the Netgen element mapping, which follows the geometry. + + Parameters + ---------- + fine_mesh + A refined Firedrake mesh with linear, and already snapped, + DMPlex coordinates. + order + Polynomial degree of the new coordinates. + + Returns + ------- + firedrake.mesh.MeshGeometry + A mesh using the re-evaluated coordinates. + """ + if order == 1 or not fine_mesh.ufl_cell().is_simplex: + # Snapping has already placed the vertices on the geometry. + fine_mesh._geometry_source = self + return fine_mesh + + dm = fine_mesh.topology_dm + fresh_mesh = createNetgenMesh(dm, self._local_mesh) + _set_netgen_coordinates( + dm, fresh_mesh, order, root_only=False, reset=True + ) + curved_mesh = _mesh_from_coordinate_dm( + fine_mesh.topology, fine_mesh.name, order + ) + curved_mesh._geometry_source = type(self)( + fresh_mesh, self.options, fine_mesh.comm, plex=dm + ) + curved_mesh._distribution_parameters = dict( + fine_mesh._distribution_parameters + ) + curved_mesh._did_reordering = fine_mesh._did_reordering + curved_mesh._tolerance = fine_mesh.tolerance + return curved_mesh diff --git a/tests/firedrake/multigrid/test_adaptive_multigrid.py b/tests/firedrake/multigrid/test_adaptive_multigrid.py index e333c3be74..8a144b45de 100644 --- a/tests/firedrake/multigrid/test_adaptive_multigrid.py +++ b/tests/firedrake/multigrid/test_adaptive_multigrid.py @@ -55,7 +55,6 @@ def coarse_mesh(request): raise NotImplementedError(f"Unrecognized mesher {mesher}") -@pytest.fixture def mh(coarse_mesh): return corner_adaptive_hierarchy(coarse_mesh, nlevels=2) @@ -275,15 +274,15 @@ def test_adapt_preserves_mesh_metadata(degree): markers = Function(FunctionSpace(mesh, "DG", 0)).assign(1) refined = mesh.refine_marked_elements(markers) - assert refined.netgen_mesh is not None - assert refined.netgen_flags == mesh.netgen_flags + assert refined._geometry_source is not None + assert refined._geometry_source.options == mesh._geometry_source.options assert refined._distribution_parameters == mesh._distribution_parameters assert refined.tolerance == mesh.tolerance assert refined.coordinates.function_space().ufl_element().degree() == degree markers = Function(FunctionSpace(refined, "DG", 0)).assign(1) twice_refined = refined.refine_marked_elements(markers) - assert twice_refined.netgen_flags == mesh.netgen_flags + assert twice_refined._geometry_source.options == mesh._geometry_source.options assert twice_refined._distribution_parameters == mesh._distribution_parameters assert twice_refined.coordinates.function_space().ufl_element().degree() == degree @@ -292,7 +291,7 @@ def test_adapt_preserves_mesh_metadata(degree): @pytest.mark.parametrize("refine", [1, 2]) def test_adapt_after_uniform_refinement(coarse_mesh, refine): """A hierarchy built by uniform refinement can be adaptively refined.""" - netgen_flags = {} if hasattr(coarse_mesh, "netgen_mesh") else None + netgen_flags = {} if coarse_mesh._geometry_source is not None else None mh = MeshHierarchy(coarse_mesh, refine, netgen_flags=netgen_flags) _assert_adapt_after_uniform_refinement(mh) diff --git a/tests/firedrake/multigrid/test_netgen_gmg.py b/tests/firedrake/multigrid/test_netgen_gmg.py index 6216907800..64cf03efc8 100644 --- a/tests/firedrake/multigrid/test_netgen_gmg.py +++ b/tests/firedrake/multigrid/test_netgen_gmg.py @@ -6,6 +6,8 @@ @pytest.fixture(params=[(2, "occ"), (2, "spline"), (2, "csg"), (3, "occ"), (3, "csg")], ids=lambda val: "-".join(map(str, val))) def ngmesh(request): + if COMM_WORLD.rank != 0: + return None dim, geo_type = request.param maxh = 0.75 if dim == 2: @@ -67,14 +69,21 @@ def test_netgen_mg(ngmesh, netgen_degree): assert not coords_space.finat_element.is_dg() errors = [] + if COMM_WORLD.rank == 0: + labels = [ + i + 1 for i, name in enumerate( + ngmesh.GetRegionNames(codim=1) + ) if name == "surface" + ] + else: + labels = None + labels = COMM_WORLD.bcast(labels, root=0) for mesh in mh[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 ["surface"]] - x = SpatialCoordinate(mesh) uexact = 1-dot(x, x) bcs = DirichletBC(V, 0, labels) diff --git a/tests/firedrake/regression/test_netgen.py b/tests/firedrake/regression/test_netgen.py index 64abf21c31..aafeebb4a6 100644 --- a/tests/firedrake/regression/test_netgen.py +++ b/tests/firedrake/regression/test_netgen.py @@ -7,14 +7,25 @@ @pytest.mark.parallel([1, 2]) def test_netgen_csg_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) + if COMM_WORLD.rank == 0: + geo = CSG2d() + geo.Add(Circle(center=(0, 0), radius=1.0, mat="mat1", bc="circle")) + ngmesh = geo.GenerateMesh(maxh=0.75) + else: + ngmesh = None # 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 + coordinate_fe, _ = mesh1.topology_dm.getCoordinateDM().getField(0) + assert coordinate_fe.getBasisSpace().getDegree() == (order, order) + coordinate_section = mesh1.topology_dm.getCoordinateSection() + edge_start, _ = mesh1.topology_dm.getDepthStratum(1) + assert coordinate_section.getDof(edge_start) == ( + mesh1.geometric_dimension * (order - 1) + ) + assert abs(assemble(1 * dx(domain=mesh1)) - np.pi) < 2.e-4 dim = mesh1.topological_dimension DG0 = FunctionSpace(mesh1, "DG", 0) markers = Function(DG0) @@ -25,6 +36,7 @@ def test_netgen_csg_mesh_high_order(): 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 + assert abs(assemble(1 * dx(domain=mesh2)) - np.pi) < 2.e-4 # Test mesh refinement: 2 refinements markers.assign(2) @@ -32,6 +44,7 @@ def test_netgen_csg_mesh_high_order(): 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 + assert abs(assemble(1 * dx(domain=mesh3)) - np.pi) < 2.e-4 def square_geometry(h, L=np.pi): @@ -215,7 +228,6 @@ def test_netgen_csg_manifold(): from netgen.csg import CSGeometry, Pnt, Sphere from netgen.meshing import MeshingParameters from netgen.meshing import MeshingStep - import netgen comm = COMM_WORLD if comm.rank == 0: @@ -224,11 +236,12 @@ def test_netgen_csg_manifold(): mp = MeshingParameters(maxh=0.05, perfstepsend=MeshingStep.MESHSURFACE) ngmesh = geo.GenerateMesh(mp=mp) else: - ngmesh = netgen.libngpy._meshing.Mesh(3) + ngmesh = None - msh = Mesh(ngmesh) + msh = Mesh(ngmesh, netgen_flags={"degree": 2}) assert msh.topological_dimension == 2 assert msh.geometric_dimension == 3 + assert msh.coordinates.ufl_element().degree() == 2 V = FunctionSpace(msh, "CG", 3) f = assemble(interpolate(Constant(1), V)) @@ -270,7 +283,6 @@ def Curve(t): @pytest.mark.parallel([1, 2]) def test_netgen_csg_high_order_integral(): from netgen.csg import CSGeometry, Pnt, Sphere - import netgen comm = COMM_WORLD if comm.rank == 0: @@ -278,7 +290,7 @@ def test_netgen_csg_high_order_integral(): geo.Add(Sphere(Pnt(0, 0, 0), 1).bc("sphere")) ngmesh = geo.GenerateMesh(maxh=0.7) else: - ngmesh = netgen.libngpy._meshing.Mesh(3) + ngmesh = None homsh = Mesh(ngmesh, netgen_flags={"degree": 2}) V = FunctionSpace(homsh, "CG", 2)