diff --git a/.github/workflows/Docker.yml b/.github/workflows/Docker.yml index 29370b7..50683ef 100644 --- a/.github/workflows/Docker.yml +++ b/.github/workflows/Docker.yml @@ -4,6 +4,8 @@ on: schedule: - cron: '30 10 7,14,21,28 * *' + workflow_dispatch: + push: paths: - Dockerfile @@ -14,8 +16,9 @@ jobs: timeout-minutes: 360 steps: - name: Checkout code - uses: actions/checkout@v2 + uses: actions/checkout@v4 - name: Login to Docker Hub + if: github.repository == 'NGSolve/ngsPETSc' uses: docker/login-action@v1 with: username: ${{ secrets.DOCKERHUB_USERNAME }} @@ -24,6 +27,7 @@ jobs: run: | docker build -t ngspetsc:latest . - name: Push Docker image + if: github.repository == 'NGSolve/ngsPETSc' run: | docker tag ngspetsc:latest urzerbinati/ngspetsc:latest docker push urzerbinati/ngspetsc:latest diff --git a/.github/workflows/ngsPETSc.yml b/.github/workflows/ngsPETSc.yml index cacac83..281b445 100644 --- a/.github/workflows/ngsPETSc.yml +++ b/.github/workflows/ngsPETSc.yml @@ -1,8 +1,8 @@ # .github/workflows/app.yaml name: ngsPETSc tests -on: +on: push: - branches-ignore: + branches-ignore: - 'no-ci/*' schedule: @@ -16,7 +16,7 @@ jobs: steps: - name: Check out repository code - uses: actions/checkout@v2 + uses: actions/checkout@v4 - name: Install ngsPETSc run: | @@ -51,20 +51,20 @@ jobs: fenicsx: runs-on: ubuntu-latest - container: dolfinx/dolfinx:latest - timeout-minutes: 30 + container: dolfinx/dolfinx:nightly + timeout-minutes: 50 steps: - name: Check out repository code - uses: actions/checkout@v2 + uses: actions/checkout@v4 - name: Install Netgen and ngsPETSc run: | - pip3 install netgen-mesher \ - && pip3 install pylint \ - && export PYTHONPATH=$PYTHONPATH:/usr/local/lib/python3.10/site-packages \ - && echo "PYTHONPATH=$PYTHONPATH" >> $GITHUB_ENV \ - && NGSPETSC_NO_INSTALL_REQUIRED=ON pip install . + pip install netgen-mesher + pip install pylint + export PYTHONPATH=$PYTHONPATH:/usr/local/lib/python3.10/site-packages + echo "PYTHONPATH=$PYTHONPATH" >> $GITHUB_ENV + NGSPETSC_NO_INSTALL_REQUIRED=ON pip install . - name: Check formatting run: | @@ -73,4 +73,4 @@ jobs: - name: Run test suite in serial run: | - pytest -v tests/test_fenicsx.py \ No newline at end of file + pytest -v tests/test_fenicsx.py diff --git a/Makefile b/Makefile index 2f822d3..a93250e 100644 --- a/Makefile +++ b/Makefile @@ -1,9 +1,9 @@ MPI_EXEC = ${PETSC_DIR}/${PETSC_ARCH}/bin/mpiexec lint: - pylint --disable=C0412,C0103,C0415,C0321,E0401,E1101,E0611,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,I1101,W0201,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc - pylint --disable=C0412,C0103,C0415,C0321,C3001,E0401,E1101,E0611,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,I1101,W0201,W0406,W0212,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc/utils + pylint --disable=C0412,C0103,C0415,C0321,E0401,E1101,E0611,R1728,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,I1101,W0201,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc + pylint --disable=C0412,C0103,C0415,C0321,C3001,E0401,E1101,E0611,R1728,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,I1101,W0201,W0406,W0212,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc/utils lint_test: - pylint --disable=C0412,C0103,C0415,C0321,E0401,E1101,E0611,R1736,R0401,R0914,R0801,R0902,R1702,R0913,R0903,R0205,R0912,R0915,I1101,W0201,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase tests + pylint --disable=C0412,C0103,C0415,C0321,E0401,E1101,E0611,R1728,R1736,R0401,R0914,R0801,R0902,R1702,R0913,R0903,R0205,R0912,R0915,I1101,W0201,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase tests test: pytest tests/test_env.py pytest tests/test_vec.py diff --git a/docs/src/Tutorials.rst b/docs/src/Tutorials.rst index 8c44c77..4881a5f 100644 --- a/docs/src/Tutorials.rst +++ b/docs/src/Tutorials.rst @@ -9,3 +9,4 @@ How to use ngsPETSc notebooks/Meshes notebooks/Krylov Solver and Preconditioners notebooks/Eigenvalue Problems + Firedrake-Netgen interface via ngsPETSc \ No newline at end of file diff --git a/ngsPETSc/plex.py b/ngsPETSc/plex.py index 9882eb6..fd50d06 100644 --- a/ngsPETSc/plex.py +++ b/ngsPETSc/plex.py @@ -33,8 +33,9 @@ class MeshMapping: ''' - def __init__(self, mesh=None, name="Default"): + def __init__(self, mesh=None, comm=MPI.COMM_WORLD, name="Default"): self.name = name + self.comm = comm if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)): self.createPETScDMPlex(mesh) elif isinstance(mesh,PETSc.DMPlex): @@ -157,7 +158,7 @@ def createPETScDMPlex(self, mesh): self.ngMesh = mesh.ngmesh else: self.ngMesh = mesh - comm = MPI.COMM_WORLD + comm = self.comm if self.ngMesh.dim == 3: if comm.rank == 0: V = self.ngMesh.Coordinates() @@ -168,7 +169,7 @@ def createPETScDMPlex(self, mesh): surfMesh, dim = True, 2 T = self.ngMesh.Elements2D().NumPy()["nodes"] T = np.array([list(np.trim_zeros(a, 'b')) for a in list(T)])-1 - plex = PETSc.DMPlex().createFromCellList(dim, T, V) + plex = PETSc.DMPlex().createFromCellList(dim, T, V, comm=comm) plex.setName(self.name) vStart, _ = plex.getDepthStratum(0) if surfMesh: @@ -186,14 +187,15 @@ def createPETScDMPlex(self, mesh): else: plex = PETSc.DMPlex().createFromCellList(3, np.zeros((0, 4), dtype=np.int32), - np.zeros((0, 3), dtype=np.double)) + np.zeros((0, 3), dtype=np.double), + comm=comm) self.petscPlex = plex elif self.ngMesh.dim == 2: if comm.rank == 0: V = self.ngMesh.Coordinates() T = self.ngMesh.Elements2D().NumPy()["nodes"] T = np.array([list(np.trim_zeros(a, 'b')) for a in list(T)])-1 - plex = PETSc.DMPlex().createFromCellList(2, T, V) + plex = PETSc.DMPlex().createFromCellList(2, T, V, comm=comm) plex.setName(self.name) vStart, _ = plex.getDepthStratum(0) # vertices for e in self.ngMesh.Elements1D(): @@ -208,5 +210,6 @@ def createPETScDMPlex(self, mesh): else: plex = PETSc.DMPlex().createFromCellList(2, np.zeros((0, 3), dtype=np.int32), - np.zeros((0, 2), dtype=np.double)) + np.zeros((0, 2), dtype=np.double), + comm=comm) self.petscPlex = plex diff --git a/ngsPETSc/utils/fenicsx.py b/ngsPETSc/utils/fenicsx.py index a16c99e..91283fc 100644 --- a/ngsPETSc/utils/fenicsx.py +++ b/ngsPETSc/utils/fenicsx.py @@ -34,7 +34,7 @@ def model_to_mesh(self, hmax: float, gdim: int = 2, dolfinx.cpp.mesh.MeshTags_int32,dolfinx.cpp.mesh.MeshTags_int32]: """Given a NetGen model, take all physical entities of the highest topological dimension and create the corresponding DOLFINx mesh. - + This function only works in serial, at the moment. Args: @@ -73,9 +73,10 @@ def model_to_mesh(self, hmax: float, gdim: int = 2, V = ngmesh.Coordinates() T = ngmesh.Elements3D().NumPy()["nodes"] T = np.array([list(np.trim_zeros(a, 'b')) for a in list(T)])-1 - ufl_domain = dolfinx.io.gmshio.ufl_mesh(_ngs_to_cells[(gdim,T.shape[1])],gdim) + ufl_domain = dolfinx.io.gmshio.ufl_mesh( + _ngs_to_cells[(gdim,T.shape[1])], gdim, dolfinx.default_real_type) cell_perm = dolfinx.cpp.io.perm_gmsh(dolfinx.cpp.mesh.to_type(str(ufl_domain.ufl_cell())), T.shape[1]) - T = T[:, cell_perm] + T = np.ascontiguousarray(T[:, cell_perm]) mesh = dolfinx.mesh.create_mesh(self.comm, T, V, ufl_domain, partitioner) return mesh diff --git a/ngsPETSc/utils/firedrake.py b/ngsPETSc/utils/firedrake.py index 24cc3f3..7462eab 100644 --- a/ngsPETSc/utils/firedrake.py +++ b/ngsPETSc/utils/firedrake.py @@ -5,11 +5,13 @@ ''' try: import firedrake as fd - import ufl + from firedrake.logging import warning + from firedrake.cython import mgimpl as impl + from firedrake.__future__ import interpolate except ImportError: fd = None - ufl = None +from fractions import Fraction import warnings import numpy as np from petsc4py import PETSc @@ -27,6 +29,15 @@ class comp: from ngsPETSc import MeshMapping +def flagsUtils(flags, option, default): + ''' + utility fuction used to parse Netgen flag options + ''' + try: + return flags[option] + except KeyError: + return default + def refineMarkedElements(self, mark): ''' This method is used to refine a mesh based on a marking function @@ -54,21 +65,41 @@ def refineMarkedElements(self, mark): self.netgen_mesh.Refine(adaptive=True) return fd.Mesh(self.netgen_mesh) return fd.Mesh(netgen.libngpy._meshing.Mesh(2)) + + elif self.geometric_dimension() == 3: + with mark.dat.vec as marked: + marked0 = marked + getIdx = self._cell_numbering.getOffset + if self.sfBCInv is not None: + getIdx = lambda x: x + _, marked0 = self.topology_dm.distributeField(self.sfBCInv, + self._cell_numbering, + marked) + if self.comm.Get_rank() == 0: + mark = marked0.getArray() + for i, el in enumerate(self.netgen_mesh.Elements3D()): + if mark[getIdx(i)]: + el.refine = True + else: + el.refine = False + self.netgen_mesh.Refine(adaptive=True) + return fd.Mesh(self.netgen_mesh) + return fd.Mesh(netgen.libngpy._meshing.Mesh(3)) else: - raise NotImplementedError("No implementation for dimension other than 2.") + raise NotImplementedError("No implementation for dimension other than 2 and 3.") -def curveField(self, order): +def curveField(self, order, tol=1e-8): ''' - This method returns a curved mesh as a Firedrake funciton. + This method returns a curved mesh as a Firedrake function. :arg order: the order of the curved mesh ''' - low_order_element = self.coordinates.function_space().ufl_element().sub_elements()[0] + low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] element = low_order_element.reconstruct(degree=order) - space = fd.VectorFunctionSpace(self, ufl.BrokenElement(element)) - newFunctionCoordinates = fd.interpolate(self.coordinates, space) - + space = fd.VectorFunctionSpace(self, fd.BrokenElement(element)) + newFunctionCoordinates = fd.assemble(interpolate(self.coordinates, space)) + self.netgen_mesh = self.comm.bcast(self.netgen_mesh, root=0) #Computing reference points using fiat fiat_element = newFunctionCoordinates.function_space().finat_element.fiat_equivalent entity_ids = fiat_element.entity_dofs() @@ -80,11 +111,8 @@ def curveField(self, order): # Assert singleton point for each node. pt, = nodes[dof].get_point_dict().keys() refPts.append(pt) - V = newFunctionCoordinates.dat.data - getIdx = self._cell_numbering.getOffset refPts = np.array(refPts) - rnd = lambda x: round(x, 8) if self.geometric_dimension() == 2: #Mapping to the physical domain physPts = np.ndarray((len(self.netgen_mesh.Elements2D()), @@ -97,17 +125,23 @@ def curveField(self, order): self.netgen_mesh.CalcElementMapping(refPts, curvedPhysPts) cellMap = newFunctionCoordinates.cell_node_map() for i, el in enumerate(self.netgen_mesh.Elements2D()): + #Inefficent code but runs only on curved elements if el.curved: - pts = [tuple(map(rnd, pts)) - for pts in physPts[i][0:refPts.shape[0]]] - dofMap = {k: v for v, k in enumerate(pts)} - p = [dofMap[tuple(map(rnd, pts))] - for pts in V[cellMap.values[getIdx(i)]][0:refPts.shape[0]]] - curvedPhysPts[i] = curvedPhysPts[i][p] - for j, datIdx in enumerate(cellMap.values[getIdx(i)][0:refPts.shape[0]]): - newFunctionCoordinates.sub(0).dat.data[datIdx] = curvedPhysPts[i][j][0] - newFunctionCoordinates.sub(1).dat.data[datIdx] = curvedPhysPts[i][j][1] - + pts = physPts[i][0:refPts.shape[0]] + bary = sum([np.array(pts[i]) for i in range(len(pts))])/len(pts) + Idx = self.locate_cell(bary) + isInMesh = (0<=Idx tol: + fd.logging.warning("Not able to curve Firedrake element {}".format(Idx)) + else: + for j, datIdx in enumerate(cellMap.values[Idx][0:refPts.shape[0]]): + newFunctionCoordinates.sub(0).dat.data[datIdx] = curvedPhysPts[i][j][0] + newFunctionCoordinates.sub(1).dat.data[datIdx] = curvedPhysPts[i][j][1] if self.geometric_dimension() == 3: #Mapping to the physical domain physPts = np.ndarray((len(self.netgen_mesh.Elements3D()), @@ -120,19 +154,29 @@ def curveField(self, order): self.netgen_mesh.CalcElementMapping(refPts, curvedPhysPts) cellMap = newFunctionCoordinates.cell_node_map() for i, el in enumerate(self.netgen_mesh.Elements3D()): + #Inefficent code but runs only on curved elements if el.curved: - pts = [tuple(map(rnd, pts)) - for pts in physPts[i][0:refPts.shape[0]]] - dofMap = {k: v for v, k in enumerate(pts)} - p = [dofMap[tuple(map(rnd, pts))] - for pts in V[cellMap.values[getIdx(i)]][0:refPts.shape[0]]] - curvedPhysPts[i] = curvedPhysPts[i][p] - for j, datIdx in enumerate(cellMap.values[getIdx(i)][0:refPts.shape[0]]): - newFunctionCoordinates.sub(0).dat.data[datIdx] = curvedPhysPts[i][j][0] - newFunctionCoordinates.sub(1).dat.data[datIdx] = curvedPhysPts[i][j][1] - newFunctionCoordinates.sub(2).dat.data[datIdx] = curvedPhysPts[i][j][2] + pts = physPts[i][0:refPts.shape[0]] + bary = sum([np.array(pts[i]) for i in range(len(pts))])/len(pts) + Idx = self.locate_cell(bary) + isInMesh = (0<=Idx tol: + warning("Not able to curve element {}, residual is: {}".format(Idx, res)) + else: + for j, datIdx in enumerate(cellMap.values[Idx][0:refPts.shape[0]]): + newFunctionCoordinates.sub(0).dat.data[datIdx] = curvedPhysPts[i][j][0] + newFunctionCoordinates.sub(1).dat.data[datIdx] = curvedPhysPts[i][j][1] + newFunctionCoordinates.sub(2).dat.data[datIdx] = curvedPhysPts[i][j][2] return newFunctionCoordinates +splitTypes = {"Alfeld": lambda x: x.SplitAlfeld(), + "Powell-Sabin": lambda x: x.SplitPowellSabin()} + class FiredrakeMesh: ''' This class creates a Firedrake mesh from Netgen/NGSolve meshes. @@ -141,20 +185,26 @@ 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=PETSc.COMM_WORLD): + def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD): self.comm = user_comm + #Parsing netgen flags + if not isinstance(netgen_flags, dict): + netgen_flags = {} + split2tets = flagsUtils(netgen_flags, "split_to_tets", False) + split = flagsUtils(netgen_flags, "split", False) + #Checking the mesh format if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)): - try: - if netgen_flags["purify_to_tets"]: - mesh.Split2Tets() - except KeyError: - warnings.warn("No purify_to_tets flag found, mesh will not be purified to tets.") - self.meshMap = MeshMapping(mesh) + if split2tets: + mesh = mesh.Split2Tets() + if split: + splitTypes[split](mesh) + self.meshMap = MeshMapping(mesh, comm=self.comm) else: raise ValueError("Mesh format not recognised.") + #We quadrilateralise the mesh or apply a transform if requested try: if netgen_flags["quad"]: - transform = PETSc.DMPlexTransform().create(comm=PETSc.COMM_WORLD) + transform = PETSc.DMPlexTransform().create(comm=self.comm) transform.setType(PETSc.DMPlexTransformType.REFINETOBOX) transform.setDM(self.meshMap.petscPlex) transform.setUp() @@ -184,17 +234,161 @@ def createFromTopology(self, topology, name): cell = topology.ufl_cell() geometric_dim = topology.topology_dm.getCoordinateDim() cell = cell.reconstruct(geometric_dimension=geometric_dim) - element = ufl.VectorElement("Lagrange", cell, 1) + element = fd.VectorElement("Lagrange", cell, 1) # Create mesh object self.firedrakeMesh = fd.MeshGeometry.__new__(fd.MeshGeometry, element) self.firedrakeMesh._init_topology(topology) self.firedrakeMesh.name = name # Adding Netgen mesh and inverse sfBC as attributes self.firedrakeMesh.netgen_mesh = self.meshMap.ngMesh - if self.comm.Get_size() > 1: + if self.firedrakeMesh.sfBC is not None: self.firedrakeMesh.sfBCInv = self.firedrakeMesh.sfBC.createInverse() else: self.firedrakeMesh.sfBCInv = None - self.firedrakeMesh.comm = self.comm + #Generating ngs to Firedrake cell index map + #Adding refine_marked_elements and curve_field methods setattr(fd.MeshGeometry, "refine_marked_elements", refineMarkedElements) setattr(fd.MeshGeometry, "curve_field", curveField) + +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 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): + ''' + This function computes the coarse to fine and fine to coarse maps + for a uniform mesh hierarchy. + ''' + refinements_per_level = 1 + lgmaps = [] + for i, m in enumerate(meshes): + no = impl.create_lgmap(m.topology_dm) + m.init() + o = impl.create_lgmap(m.topology_dm) + m.topology_dm.setRefineLevel(i) + lgmaps.append((no, o)) + 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 tollerance adopted in snapToNetgenDMPlex. + -refinement_type, the refinment type to be used: uniform (default), Alfeld + ''' + if mesh.geometric_dimension() == 3: + raise NotImplementedError("Netgen hierachies are only implemented for 2D meshes.") + ngmesh = mesh.netgen_mesh + comm = mesh.comm + #Parsing netgen flags + if not isinstance(flags, dict): + flags = {} + order = flagsUtils(flags, "degree", 1) + if isinstance(order, int): + order= [order]*(levs+1) + tol = flagsUtils(flags, "tol", 1e-8) + refType = flagsUtils(flags, "refinement_type", "uniform") + #Firedrake quoantities + meshes = [] + coarse_to_fine_cells = [] + fine_to_coarse_cells = [None] + params = {"partition": False} + #We construct the unrefined linear mesh + if mesh.comm.size > 1 and mesh._grown_halos: + raise RuntimeError("Cannot refine parallel overlapped meshes ") + #We curve the mesh + if order[0]>1: + mesh = fd.Mesh(mesh.curve_field(order=order[0], tol=tol), + distribution_parameters=params, comm=comm) + meshes += [mesh] + cdm = meshes[-1].topology_dm + 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 + snapToNetgenDMPlex(ngmesh, rdm) + #We construct a Firedrake mesh from the DMPlex mesh + mesh = fd.Mesh(rdm, dim=meshes[-1].ufl_cell().geometric_dimension(), reorder=False, + distribution_parameters=params, comm=comm) + mesh.netgen_mesh = ngmesh + #We curve the mesh + if order[l+1] > 1: + mesh = fd.Mesh(mesh.curve_field(order=order[l+1], tol=1e-8), + distribution_parameters=params, comm=comm) + meshes += [mesh] + #We populate the coarse to fine map + coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes) + return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, + 1, nested=False) diff --git a/setup.py b/setup.py index a73e4c4..ae8902b 100644 --- a/setup.py +++ b/setup.py @@ -9,11 +9,13 @@ install_requires = [ 'petsc4py', 'mpi4py', + 'numpy', 'pytest', #For testing 'pylint', #For formatting ] else: install_requires=[ + 'netgen', 'ngsolve', 'petsc4py', 'mpi4py', @@ -38,7 +40,7 @@ "Operating System :: OS Independent", ], packages=["ngsPETSc", "ngsPETSc.utils"], - python_requires = ">=3.10", + python_requires = ">=3.9", install_requires=install_requires ) diff --git a/tests/.gitignore b/tests/.gitignore new file mode 100644 index 0000000..9e68810 --- /dev/null +++ b/tests/.gitignore @@ -0,0 +1 @@ +XDMF diff --git a/tests/test_fenicsx.py b/tests/test_fenicsx.py index f195147..3374db0 100644 --- a/tests/test_fenicsx.py +++ b/tests/test_fenicsx.py @@ -42,8 +42,8 @@ def test_poisson_netgen(): geo = SplineGeometry() geo.AddRectangle((0,0),(np.pi,np.pi)) geoModel = ngfx.GeometricModel(geo, MPI.COMM_WORLD) - msh =geoModel.model_to_mesh(hmax=0.1) - V = fem.FunctionSpace(msh, ("Lagrange", 2)) + msh = geoModel.model_to_mesh(hmax=0.1) + V = fem.functionspace(msh, ("Lagrange", 2)) #pylint: disable=E1120 facetsLR = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1), marker=lambda x: np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], np.pi)))