From 8ab35989d32934bee75380fc05d6ed3eb66c1647 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 24 Apr 2024 13:01:08 +0100 Subject: [PATCH 1/3] Relabeling netgen_flags Signed-off-by: Umberto Zerbinati --- ngsPETSc/plex.py | 24 +++++++++++++++++------- ngsPETSc/utils/firedrake.py | 5 +++-- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/ngsPETSc/plex.py b/ngsPETSc/plex.py index fd50d06..23a41a3 100644 --- a/ngsPETSc/plex.py +++ b/ngsPETSc/plex.py @@ -33,11 +33,11 @@ class MeshMapping: ''' - def __init__(self, mesh=None, comm=MPI.COMM_WORLD, name="Default"): + def __init__(self, mesh=None, comm=MPI.COMM_WORLD, name="default", labels={}): #pylint: disable=W0102 self.name = name self.comm = comm if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)): - self.createPETScDMPlex(mesh) + self.createPETScDMPlex(mesh, labels=labels) elif isinstance(mesh,PETSc.DMPlex): self.createNGSMesh(mesh) else: @@ -146,13 +146,13 @@ def createNGSMesh(self, plex): else: raise NotImplementedError("No implementation for dimension greater than 3.") - def createPETScDMPlex(self, mesh): + def createPETScDMPlex(self, mesh, labels): ''' This function generate an PETSc DMPlex from a Netgen/NGSolve mesh object - :arg plex: the Netgen/NGSolve mesh object to be converted into a PETSc + :arg mesh: the Netgen/NGSolve mesh object to be converted into a PETSc DMPlex. - + :arg labels: Dictionary with the new labels for the sets. ''' if isinstance(mesh,ngs.comp.Mesh): self.ngMesh = mesh.ngmesh @@ -192,6 +192,14 @@ def createPETScDMPlex(self, mesh): self.petscPlex = plex elif self.ngMesh.dim == 2: if comm.rank == 0: + # Map Netgen string to ids for the labels + newLabels = [] + if len(labels) > 0: + for key in labels.keys(): + for i, s in enumerate(self.ngMesh.GetRegionNames(dim=1)): + if key == s: + newLabels = newLabels + [(i+1, labels[key])] + newLabels = dict(newLabels) 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 @@ -204,8 +212,10 @@ def createPETScDMPlex(self, mesh): if not (1 == self.ngMesh.Elements2D().NumPy()["index"]).all(): for e in self.ngMesh.Elements2D(): join = plex.getFullJoin([vStart+v.nr-1 for v in e.vertices]) - plex.setLabelValue(CELL_SETS_LABEL, join[0], int(e.index)) - + if e.index in newLabels: + plex.setLabelValue(CELL_SETS_LABEL, join[0], int(newLabels[e.index])) + else: + plex.setLabelValue(CELL_SETS_LABEL, join[0], int(e.index)) self.petscPlex = plex else: plex = PETSc.DMPlex().createFromCellList(2, diff --git a/ngsPETSc/utils/firedrake.py b/ngsPETSc/utils/firedrake.py index 2645e33..28bbc87 100644 --- a/ngsPETSc/utils/firedrake.py +++ b/ngsPETSc/utils/firedrake.py @@ -211,6 +211,7 @@ def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD): split = flagsUtils(netgen_flags, "split", False) quad = flagsUtils(netgen_flags, "quad", False) optMoves = flagsUtils(netgen_flags, "optimisation_moves", False) + labels = flagsUtils(netgen_flags, "labels", {}) #Checking the mesh format if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)): if split2tets: @@ -227,7 +228,7 @@ def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD): 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) + self.meshMap = MeshMapping(mesh, comm=self.comm, labels=labels) #We apply the DMPLEX transform if quad: newplex = splitToQuads(self.meshMap.petscPlex, mesh.dim, comm=self.comm) @@ -359,7 +360,7 @@ def NetgenHierarchy(mesh, levs, flags): -tol, geometric tollerance adopted in snapToNetgenDMPlex. -refinement_type, the refinment type to be used: uniform (default), Alfeld ''' - if mesh.geometric_dimension() == 3: + if mesh.dim == 3: raise NotImplementedError("Netgen hierachies are only implemented for 2D meshes.") ngmesh = mesh.netgen_mesh comm = mesh.comm From 371e894764f86884d7c5852ef83c8fae460dc1b2 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 24 Apr 2024 13:07:31 +0100 Subject: [PATCH 2/3] Working on 1D bc not 2d area ... Signed-off-by: Umberto Zerbinati --- ngsPETSc/plex.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/ngsPETSc/plex.py b/ngsPETSc/plex.py index 23a41a3..934bba5 100644 --- a/ngsPETSc/plex.py +++ b/ngsPETSc/plex.py @@ -208,14 +208,14 @@ def createPETScDMPlex(self, mesh, labels): vStart, _ = plex.getDepthStratum(0) # vertices for e in self.ngMesh.Elements1D(): join = plex.getJoin([vStart+v.nr-1 for v in e.vertices]) - plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index)) + if e.index in newLabels: + plex.setLabelValue(FACE_SETS_LABEL, join[0], int(newLabels[e.index])) + else: + plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index)) if not (1 == self.ngMesh.Elements2D().NumPy()["index"]).all(): for e in self.ngMesh.Elements2D(): join = plex.getFullJoin([vStart+v.nr-1 for v in e.vertices]) - if e.index in newLabels: - plex.setLabelValue(CELL_SETS_LABEL, join[0], int(newLabels[e.index])) - else: - plex.setLabelValue(CELL_SETS_LABEL, join[0], int(e.index)) + plex.setLabelValue(CELL_SETS_LABEL, join[0], int(e.index)) self.petscPlex = plex else: plex = PETSc.DMPlex().createFromCellList(2, From 22fc5221af8c9e1507300233db312027e0a7afa8 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Thu, 25 Apr 2024 15:26:55 +0100 Subject: [PATCH 3/3] Also work in 3D Signed-off-by: Umberto Zerbinati --- ngsPETSc/plex.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/ngsPETSc/plex.py b/ngsPETSc/plex.py index 934bba5..5f87cb1 100644 --- a/ngsPETSc/plex.py +++ b/ngsPETSc/plex.py @@ -159,6 +159,14 @@ def createPETScDMPlex(self, mesh, labels): else: self.ngMesh = mesh comm = self.comm + # Map Netgen string to ids for the labels + newLabels = [] + if len(labels) > 0: + for key in labels.keys(): + for i, s in enumerate(self.ngMesh.GetRegionNames(dim=self.ngMesh.dim-1)): + if key == s: + newLabels = newLabels + [(i+1, labels[key])] + newLabels = dict(newLabels) if self.ngMesh.dim == 3: if comm.rank == 0: V = self.ngMesh.Coordinates() @@ -179,7 +187,10 @@ def createPETScDMPlex(self, mesh, labels): else: for e in self.ngMesh.Elements2D(): join = plex.getFullJoin([vStart+v.nr-1 for v in e.vertices]) - plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index)) + if e.index in newLabels: + plex.setLabelValue(FACE_SETS_LABEL, join[0], int(newLabels[e.index])) + else: + plex.setLabelValue(FACE_SETS_LABEL, join[0], int(e.index)) for e in self.ngMesh.Elements1D(): join = plex.getJoin([vStart+v.nr-1 for v in e.vertices]) plex.setLabelValue(EDGE_SETS_LABEL, join[0], int(e.index)) @@ -192,14 +203,6 @@ def createPETScDMPlex(self, mesh, labels): self.petscPlex = plex elif self.ngMesh.dim == 2: if comm.rank == 0: - # Map Netgen string to ids for the labels - newLabels = [] - if len(labels) > 0: - for key in labels.keys(): - for i, s in enumerate(self.ngMesh.GetRegionNames(dim=1)): - if key == s: - newLabels = newLabels + [(i+1, labels[key])] - newLabels = dict(newLabels) 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