From 941486a6fe84be7001c376e34dad4bb262a817bc Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Sat, 17 Feb 2024 14:45:38 +0000 Subject: [PATCH] Added mesh smoothing Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake.py | 53 +++++++++++++++++++++++-------------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/ngsPETSc/utils/firedrake.py b/ngsPETSc/utils/firedrake.py index 7462eab..26495df 100644 --- a/ngsPETSc/utils/firedrake.py +++ b/ngsPETSc/utils/firedrake.py @@ -18,6 +18,7 @@ import netgen import netgen.meshing as ngm +from netgen.meshing import MeshingParameters try: import ngsolve as ngs except ImportError: @@ -174,6 +175,20 @@ def curveField(self, order, tol=1e-8): newFunctionCoordinates.sub(2).dat.data[datIdx] = curvedPhysPts[i][j][2] return newFunctionCoordinates +def splitToQuads(plex, dim, comm): + ''' + This method splits a Netgen mesh to quads, using a PETSc transform. + ''' + if dim == 2: + transform = PETSc.DMPlexTransform().create(comm=comm) + transform.setType(PETSc.DMPlexTransformType.REFINETOBOX) + transform.setDM(plex) + transform.setUp() + newplex = transform.apply(plex) + return newplex + else: + raise RuntimeError("Splitting to quads is only possible for 2D meshes.") + splitTypes = {"Alfeld": lambda x: x.SplitAlfeld(), "Powell-Sabin": lambda x: x.SplitPowellSabin()} @@ -192,35 +207,33 @@ def __init__(self, mesh, netgen_flags, user_comm=fd.COMM_WORLD): netgen_flags = {} split2tets = flagsUtils(netgen_flags, "split_to_tets", False) split = flagsUtils(netgen_flags, "split", False) + quad = flagsUtils(netgen_flags, "quad", False) + optMoves = flagsUtils(netgen_flags, "optimisation_moves", False) #Checking the mesh format if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)): if split2tets: + #Splits 2 tets 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)) + print("Optimising mesh") + 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 + if quad: + newplex = splitToQuads(self.meshMap.petscPlex, mesh.dim, comm=self.comm) + self.meshMap = MeshMapping(newplex) 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=self.comm) - transform.setType(PETSc.DMPlexTransformType.REFINETOBOX) - transform.setDM(self.meshMap.petscPlex) - transform.setUp() - newplex = transform.apply(self.meshMap.petscPlex) - self.meshMap = MeshMapping(newplex) - except KeyError: - warnings.warn("No quad flag found, mesh will not be quadrilateralised.") - try: - if netgen_flags["transform"] is not None: - transform = netgen_flags["transform"] - transform.setDM(self.meshMap.petscPlex) - transform.setUp() - newplex = transform.apply(self.meshMap.petscPlex) - self.meshMap = MeshMapping(newplex) - except KeyError: - warnings.warn("No PETSc transform found, mesh will not be transformed.") def createFromTopology(self, topology, name): '''