From e0afb7ef30511643bf6d7717f7cdac2744822146 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Sun, 7 Jan 2024 22:55:41 +0000 Subject: [PATCH] First attempt with adaptive mesh refinement Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake.py | 74 +++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/ngsPETSc/utils/firedrake.py b/ngsPETSc/utils/firedrake.py index 5bc065c..14d4f35 100644 --- a/ngsPETSc/utils/firedrake.py +++ b/ngsPETSc/utils/firedrake.py @@ -11,6 +11,8 @@ import warnings import numpy as np from petsc4py import PETSc +from fractions import Fraction +from copy import deepcopy import netgen import netgen.meshing as ngm @@ -216,3 +218,75 @@ def createFromTopology(self, topology, name): self.firedrakeMesh.comm = self.comm setattr(fd.MeshGeometry, "refine_marked_elements", refineMarkedElements) setattr(fd.MeshGeometry, "curve_field", curveField) + +def NetgenHierarchy(ngmesh, levs, order=1, comm=fd.COMM_WORLD): + ''' + 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 + ''' + meshes = [] + refinements_per_level = 1 + coarse_to_fine_cells = [] + fine_to_coarse_cells = [None] + #We construct the unrefined linear mesh + mesh = fd.Mesh(ngmesh, comm=comm) + if mesh.comm.size > 1 and mesh._grown_halos: + raise RuntimeError("Cannot refine parallel overlapped meshes ") + #Computing useful quantities + fd.FunctionSpace(mesh, "DG", 0) + coarse_index = mesh._cell_numbering.getOffset + if ngmesh.dim == 2: + number_coarse_cells = len(ngmesh.Elements2D()) + elif ngmesh.dim == 3: + pass + #Curving linear mesh + mesh = fd.Mesh(mesh.curve_field(order=order)) + meshes += [mesh] + for lv in range(levs): + #We refine the netgen mesh uniformly + ngmesh.Refine(adaptive=True) + #We construct the refined linear mesh and curve it + mesh = fd.Mesh(ngmesh, comm=comm) + fd.FunctionSpace(mesh, "DG", 0) + fine_index = lambda x: x + if ngmesh.dim == 2: + number_fine_cells = len(ngmesh.Elements2D()) + elif ngmesh.dim == 3: + pass + mesh = fd.Mesh(mesh.curve_field(order=order)) + meshes += [mesh] + #We populate the coarse to fine map + fine_cells_per_coarse_cell = number_fine_cells // number_coarse_cells + c2f = np.zeros((number_coarse_cells,fine_cells_per_coarse_cell),dtype=np.int32) + f2c = np.zeros((number_fine_cells,1),dtype=np.int32) + c2f_counter = [0]*number_coarse_cells + if ngmesh.dim == 2: + for i in range(len(ngmesh.parentsurfaceelements)): + if i < number_coarse_cells: + f2c[fine_index(i)][0] = coarse_index(i) + c2f[coarse_index(i)][0] = fine_index(i) + c2f_counter[coarse_index(i)] += 1 + else: + k = (ngmesh.parentsurfaceelements[i]-1)%number_coarse_cells + + c2f[coarse_index(k)][c2f_counter[coarse_index(k)]] = fine_index(i) + c2f_counter[coarse_index(k)] += 1 + f2c[fine_index(i)][0] = coarse_index(k) + coarse_to_fine_cells.append(c2f) + fine_to_coarse_cells.append(f2c) + coarse_index = fine_index + number_coarse_cells = number_fine_cells + + 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 meshes, coarse_to_fine_cells, fine_to_coarse_cells + return fd.HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells, + refinements_per_level, nested=False) + + +