diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index affaf0fc3f..fe51c69082 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -186,7 +186,7 @@ jobs: : # 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/connorjward/ngsPETSc.git@connorjward/move-firedrake-code netgen-mesher netgen-occt : # We have to pass '--no-build-isolation' to use a custom petsc4py EXTRA_BUILD_ARGS='--no-isolation' diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 87757a2e54..bc94e08144 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np import ctypes import os @@ -2834,6 +2836,170 @@ def mark_entities(self, f, label_value, label_name=None): """ self.topology.mark_entities(f.topological, label_value, label_name) + # NOTE: Used to be curve_field and return coordinates (deprecate this) + # TODO: 'snap_to_netgen_mesh'? + @PETSc.Log.EventDecorator() + def curve(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=False) -> MeshGeometry: + ''' + This method returns a curved mesh as a Firedrake function. + + :arg order: the order of the curved mesh. + :arg permutation_tol: tolerance used to construct the permutation of the reference element. + :arg location_tol: tolerance used to locate the cell a point belongs to. + :arg cg_field: return a CG function field representing the mesh, rather than the + default DG field. + + ''' + import ngsPETSc + + # 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 + else: + ng_element = self.netgen_mesh.Elements3D + ng_dimension = len(ng_element()) + geom_dim = self.geometric_dimension() + + # Construct the mesh as a Firedrake function + if cg_field: + firedrake_space = fd.VectorFunctionSpace(self, "CG", order) + 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(interpolate(self.coordinates, firedrake_space)) + + # 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) + 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) + ) + + # Broadcast curved cell point data + self.comm.Bcast(physical_space_points, root=0) + self.comm.Bcast(curved_space_points, root=0) + cell_node_map = new_coordinates.cell_node_map() + + # Select only the points in curved cells + barycentres = np.average(physical_space_points, axis=1) + ng_index = [*map(lambda x: self.locate_cell(x, tolerance=location_tol), barycentres)] + + # 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]) + + # Find the correct coordinate permutation for each cell + # NB: Coordinates must be cast to real when running Firedrake in complex mode + permutation = find_permutation( + physical_space_points, + new_coordinates.dat.data[pyop2_index].reshape( + physical_space_points.shape + ).astype(np.float64, copy=False), + tol=permutation_tol + ) + + # Apply the permutation to each cell in turn + for ii, p in enumerate(curved_space_points): + curved_space_points[ii] = p[permutation[ii]] + + # Assign the curved coordinates to the dat + new_coordinates.dat.data[pyop2_index] = curved_space_points.reshape(-1, geom_dim) + + return new_coordinates + + def refine_marked_elements(self, mark, netgen_flags={}): + ''' + This method is used to refine a mesh based on a marking function + which is a Firedrake DG0 function. + + :arg mark: the marking function which is a Firedrake DG0 function. + :arg netgen_flags: the dictionary of flags to be passed to ngsPETSc. + It includes the option: + - refine_faces, which is a boolean specifyiong if you want to refine faces. + + ''' + 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 not in {2, 3}: + raise NotImplementedError( + "No implementation for dimension other than 2 and 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 #pylint: disable=C3001 + _, 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 Mesh(self.netgen_mesh, distribution_parameters=DistParams, comm=self.comm) + else: + return Mesh(netgen.libngpy._meshing.Mesh(dim), + distribution_parameters=DistParams, comm=self.comm) + @PETSc.Log.EventDecorator() def make_mesh_from_coordinates(coordinates, name, tolerance=0.5): @@ -3109,11 +3275,9 @@ 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 netgen and isinstance(meshfile, netgen.libngpy._meshing.Mesh): - netgen_firedrake_mesh.createFromTopology(topology, name=name, comm=user_comm) - mesh = netgen_firedrake_mesh.firedrakeMesh - else: - mesh = make_mesh_from_mesh_topology(topology, name) + mesh.netgen_mesh = netgen_firedrake_mesh mesh.submesh_parent = submesh_parent mesh._tolerance = tolerance return mesh diff --git a/firedrake/mg/mesh.py b/firedrake/mg/mesh.py index 791775054d..f23b40a192 100644 --- a/firedrake/mg/mesh.py +++ b/firedrake/mg/mesh.py @@ -9,6 +9,7 @@ import firedrake.cython.dmcommon as dmcommon from firedrake.utils import cached_property from firedrake.cython import mgimpl as impl +from .netgen import NetgenHierarchy from .utils import set_level __all__ = ("HierarchyBase", "MeshHierarchy", "ExtrudedMeshHierarchy", "NonNestedHierarchy", @@ -118,14 +119,8 @@ def MeshHierarchy(mesh, refinement_levels, A :py:class:`HierarchyBase` object representing the mesh hierarchy. """ - + # TODO: This is such a weird init procedure if (isinstance(netgen_flags, bool) and netgen_flags) or isinstance(netgen_flags, dict): - try: - from ngsPETSc import NetgenHierarchy - except ImportError: - raise ImportError("Unable to import netgen and ngsPETSc. Please ensure that netgen and ngsPETSc " - "are installed and available to Firedrake (see " - "https://www.firedrakeproject.org/install.html#netgen).") if hasattr(mesh, "netgen_mesh"): return NetgenHierarchy(mesh, refinement_levels, flags=netgen_flags) else: diff --git a/firedrake/mg/netgen.py b/firedrake/mg/netgen.py new file mode 100644 index 0000000000..f76fd24fa5 --- /dev/null +++ b/firedrake/mg/netgen.py @@ -0,0 +1,280 @@ +from fractions import Fraction + +import numpy as np +import ufl +from petsc4py import PETSc + +import firedrake as fd +from firedrake.cython import mgimpl as impl, dmcommon +from firedrake import dmhooks + + +__all__ = () + + +def snapToNetgenDMPlex(ngmesh, petscPlex): + ''' + This function snaps the coordinates of a DMPlex mesh to the coordinates of a Netgen mesh. + ''' + if petscPlex.getDimension() == 2: + ngCoordinates = ngmesh.Coordinates() + petscCoordinates = petscPlex.getCoordinatesLocal().getArray().reshape(-1, ngmesh.dim) + for i, pt in enumerate(petscCoordinates): + j = np.argmin(np.sum((ngCoordinates - pt)**2, axis=1)) + petscCoordinates[i] = ngCoordinates[j] + petscPlexCoordinates = petscPlex.getCoordinatesLocal() + petscPlexCoordinates.setArray(petscPlexCoordinates) + petscPlex.setCoordinatesLocal(petscPlexCoordinates) + else: + raise NotImplementedError("Snapping to Netgen meshes is only implemented for 2D meshes.") + + +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(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: + 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) + 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) + + +def uniformRefinementRoutine(ngmesh, cdm): + ''' + Routing called inside of NetgenHierarchy to compute refined ngmesh and plex. + ''' + #We refine the netgen mesh uniformly + ngmesh.Refine(adaptive=False) + #We refine the DMPlex mesh uniformly + cdm.setRefinementUniform(True) + rdm = cdm.refine() + rdm.removeLabel("pyop2_core") + rdm.removeLabel("pyop2_owned") + rdm.removeLabel("pyop2_ghost") + return (rdm, 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): + ''' + This function creates a Firedrake mesh hierarchy from Netgen/NGSolve meshes. + + :arg mesh: the Netgen/NGSolve mesh + :arg levs: the number of levels in the hierarchy + :arg netgen_flags: either a bool or a dictionray 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 + ''' + try: + from netgen.meshing import MeshingParameters + except ImportError: + raise ImportError("Unable to import netgen and ngsPETSc. Please ensure that netgen and ngsPETSc " + "are installed and available to Firedrake (see " + "https://www.firedrakeproject.org/install.html#netgen).") + + if mesh.geometric_dimension() == 3: + raise NotImplementedError("Netgen hierachies are only implemented for 2D meshes.") + comm = mesh.comm + #Parsing netgen flags + if not isinstance(flags, dict): + flags = {} + order = flags.get("degree", 1) + if isinstance(order, int): + order= [order]*(levs+1) + permutation_tol = flags.get("tol", 1e-8) + 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", False) + nested = flags.get("nested", snap in ["coarse"]) + #Firedrake quoantities + meshes = [] + lgmaps = [] + params = {"partition": False} + #We curve the mesh + if order[0]>1: + ho_field = mesh.curve_field( + order=order[0], + permutation_tol=permutation_tol, + cg_field=cg + ) + temp = fd.Mesh(ho_field,distribution_parameters=params, comm=comm) + temp.netgen_mesh = mesh.netgen_mesh + temp._tolerance = mesh.tolerance + mesh = temp + # 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") + 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 += [mesh] + ngmesh = mesh.netgen_mesh + for l in range(levs): + #Streightening the mesh + ngmesh.Curve(1) + rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm) + cdm = rdm + #We snap the mesh to the Netgen mesh + if snap == "geometry": + snapToNetgenDMPlex(ngmesh, rdm) + #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=params, 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: + ngmesh.OptimizeMesh2d(MeshingParameters(optimize2d=optMoves)) + elif mesh.dim == 3: + ngmesh.OptimizeVolumeMesh(MeshingParameters(optimize3d=optMoves)) + else: + raise ValueError("Only 2D and 3D meshes can be optimised.") + mesh.netgen_mesh = ngmesh + #We curve the mesh + if order[l+1] > 1: + if snap == "geometry": + mesh = fd.Mesh( + mesh.curve_field(order=order[l+1], permutation_tol=permutation_tol), + distribution_parameters=params, + comm=comm + ) + elif snap == "coarse": + mesh = snapToCoarse(ho_field, mesh, order[l+1], snap_smoothing, cg) + mesh.topology_dm.setRefineLevel(1 + l) + meshes += [mesh] + #We populate the coarse to fine map + coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes, lgmaps) + return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, + 1, nested=nested) diff --git a/firedrake/utils.py b/firedrake/utils.py index 901e2694ae..164029bffa 100644 --- a/firedrake/utils.py +++ b/firedrake/utils.py @@ -168,3 +168,15 @@ def wrapper(*args, **kwargs): return fn(*args, **kwargs) return wrapper return decorator + +# issue: sometimes not just import ngsPETSc here +# +# def import_ngsPETSc(): +# try: +# import ngsPETSc +# except ImportError: +# raise ImportError("Unable to import netgen and ngsPETSc. Please ensure that netgen and ngsPETSc " +# "are installed and available to Firedrake (see " +# "https://www.firedrakeproject.org/install.html#netgen).") +# else: +# return ngsPETSc