From df466b2df8279d19f5a4fb98d51786cc93f83175 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Tue, 6 Feb 2024 22:15:22 +0000 Subject: [PATCH] splits option in Netgen falgs Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake.py | 100 ++++++++++++++++++++++-------------- 1 file changed, 62 insertions(+), 38 deletions(-) diff --git a/ngsPETSc/utils/firedrake.py b/ngsPETSc/utils/firedrake.py index ead8021..7462eab 100644 --- a/ngsPETSc/utils/firedrake.py +++ b/ngsPETSc/utils/firedrake.py @@ -29,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 @@ -165,6 +174,9 @@ def curveField(self, order, tol=1e-8): 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. @@ -175,15 +187,21 @@ class FiredrakeMesh: ''' 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.") + 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=self.comm) @@ -248,15 +266,6 @@ def snapToNetgenDMPlex(ngmesh, petscPlex): else: raise NotImplementedError("Snapping to Netgen meshes is only implemented for 2D meshes.") -def flagsUtils(flags, option, default): - ''' - utility fuction used to parse Netgen flag options - ''' - try: - return flags[option] - except KeyError: - return default - def uniformRefinementRoutine(ngmesh, cdm): ''' Routing called inside of NetgenHierarchy to compute refined ngmesh and plex. @@ -271,6 +280,33 @@ def uniformRefinementRoutine(ngmesh, cdm): 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. @@ -285,8 +321,15 @@ def alfeldRefinementRoutine(ngmesh, cdm): rdm = tr.apply(cdm) return (rdm, ngmesh) -refinementTypes = {"uniform": uniformRefinementRoutine, - "Alfeld": alfeldRefinementRoutine} +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): ''' @@ -317,7 +360,6 @@ def NetgenHierarchy(mesh, levs, flags): refType = flagsUtils(flags, "refinement_type", "uniform") #Firedrake quoantities meshes = [] - refinements_per_level = 1 coarse_to_fine_cells = [] fine_to_coarse_cells = [None] params = {"partition": False} @@ -333,7 +375,7 @@ def NetgenHierarchy(mesh, levs, flags): for l in range(levs): #Streightening the mesh ngmesh.Curve(1) - rdm, ngmesh = refinementTypes[refType](ngmesh, cdm) + rdm, ngmesh = refinementTypes[refType][0](ngmesh, cdm) cdm = rdm #We snap the mesh to the Netgen mesh snapToNetgenDMPlex(ngmesh, rdm) @@ -347,24 +389,6 @@ def NetgenHierarchy(mesh, levs, flags): distribution_parameters=params, comm=comm) meshes += [mesh] #We populate the coarse to fine map - 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)) + coarse_to_fine_cells, fine_to_coarse_cells = refinementTypes[refType][1](meshes) return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, - refinements_per_level, nested=False) + 1, nested=False)