Skip to content

Tidy up plex.py #106

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pbrubeck wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Tidy up plex.py #106

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
179373d
Removed JIT (I've not seen significant improvement in performance).
Mar 12, 2026
9c625f0
Tidy up plex.py
pbrubeck Mar 12, 2026
314ee67
remove unnecesary vertex permutation for positive Jacobian
pbrubeck Mar 18, 2026
cf66a98
cleanup
pbrubeck Mar 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
222 changes: 59 additions & 163 deletions ngsPETSc/plex.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,199 +16,100 @@
class comp:
"dummy class"
Mesh = type(None)
try:
from numba import jit
numba = True
except: # pylint: disable=bare-except
jit = None
numba = False

FACE_SETS_LABEL = "Face Sets"
CELL_SETS_LABEL = "Cell Sets"
EDGE_SETS_LABEL = "Edge Sets"

if numba:
@jit(nopython=True, cache=True)
def build2DCells(coordinates, sIndicies, cell_indicies, start_end):
"""
JIT Method used to construct 2D cells
"""
cStart = start_end[0]
cEnd = start_end[1]
cells = np.zeros((cell_indicies.shape[0], 3))
for i in range(cStart,cEnd):
sIndex = sIndicies[i]
if len(sIndex)==3:
cell = list(set(cell_indicies[i]))
A = np.zeros((2,2))
A[0,0] = (coordinates[cell[1]]-coordinates[cell[0]])[0]
A[0,1] = (coordinates[cell[1]]-coordinates[cell[0]])[1]
A[1,0] = (coordinates[cell[2]]-coordinates[cell[1]])[0]
A[1,1] = (coordinates[cell[2]]-coordinates[cell[1]])[1]
if np.linalg.det(A) > 0:
cells[i] = cell
else:
cells[i] = np.array([cell[0], cell[2], cell[1]])
else:
raise RuntimeError("We only support triangles.")
return cells

@jit(nopython=True, cache=True)
def build3DCells(coordinates, sIndicies, cell_indicies, start_end):
"""
JIT Method used to construct 3D cells
"""
cStart = start_end[0]
cEnd = start_end[1]
cells = np.zeros((cell_indicies.shape[0], 4))
for i in range(cStart,cEnd):
sIndex = sIndicies[i]
if len(sIndex)==4:
cell = list(set(cell_indicies[i]))
A = np.zeros((3,3))
A[0,0] = (coordinates[cell[1]]-coordinates[cell[0]])[0]
A[0,1] = (coordinates[cell[1]]-coordinates[cell[0]])[1]
A[0,2] = (coordinates[cell[1]]-coordinates[cell[0]])[2]
A[1,0] = (coordinates[cell[2]]-coordinates[cell[1]])[0]
A[1,1] = (coordinates[cell[2]]-coordinates[cell[1]])[1]
A[1,2] = (coordinates[cell[2]]-coordinates[cell[1]])[2]
A[2,0] = (coordinates[cell[3]]-coordinates[cell[2]])[0]
A[2,1] = (coordinates[cell[3]]-coordinates[cell[2]])[1]
A[2,2] = (coordinates[cell[3]]-coordinates[cell[2]])[2]
if np.linalg.det(A) > 0:
cells[i] = cell
else:
cells[i] = np.array([cell[0], cell[1], cell[2], cell[2]])
else:
raise RuntimeError("We only support tets.")
return cells
else:
def build2DCells(coordinates, sIndicies, cell_indicies, start_end):
"""
Method used to construct 2D cells
"""
cStart, cEnd = start_end
cells = np.zeros((cell_indicies.shape[0], 3))
for i in range(cStart,cEnd):
sIndex = sIndicies[i]
if len(sIndex)==3:
cell = list(set(cell_indicies[i]))
A = np.zeros((2,2))
for i in range(2):
A[i, :] = coordinates[cell[i+1]] - coordinates[cell[i]]
if np.linalg.det(A) > 0:
cells[i] = cell
else:
cells[i] = np.array([cell[0], cell[2], cell[1]])
else:
raise RuntimeError("We only support triangles.")
return cells

def build3DCells(coordinates, sIndicies, cell_indicies, start_end):
"""
Method used to construct 3D cells
"""
cStart, cEnd = start_end
cells = np.zeros((cell_indicies.shape[0], 4))
for i in range(cStart,cEnd):
sIndex = sIndicies[i]
if len(sIndex)==4:
cell = list(set(cell_indicies[i]))
A = np.zeros((3,3))
for i in range(3):
A[i, :] = coordinates[cell[i+1]] - coordinates[cell[i]]
if np.linalg.det(A) > 0:
cells[i] = cell
else:
cells[i] = np.array([cell[0], cell[1], cell[3], cell[2]])
else:
raise RuntimeError("We only support tets.")
return cells

def buildCells(coordinates, plex):

Check failure on line 25 in ngsPETSc/plex.py

View workflow job for this annotation

GitHub Actions / lint

W0613 

ngsPETSc/plex.py:25:15: W0613 Unused argument 'coordinates'
"""
Return a numpy.array with the vertices of each cell

:arg coordinates: vertices coordinates
:arg plex: PETSc DMPlex

"""
cStart, cEnd = plex.getHeightStratum(0)
vStart, vEnd = plex.getDepthStratum(0)

cells = [[v-vStart for v in plex.getAdjacency(c) if vStart <= v < vEnd]
for c in range(cStart, cEnd)]
cells = np.array(cells)
return cells


def create2DNetgenMesh(ngMesh, coordinates, plex, geoInfo):
"""
Method used to generate 2D NetgenMeshes

:arg ngMesh: the netgen mesh to be populated
:arg coordinates: vertices coordinates
:arg coordinates: vertices coordinates
:arg plex: PETSc DMPlex
:arg geoInfo: geometric information assosciated with the Netgen mesh

"""
ngMesh.AddPoints(coordinates)
cStart,cEnd = plex.getHeightStratum(0)
vStart, _ = plex.getHeightStratum(2)
# Outside of jitted loop we put all calls to plex
sIndicies = [plex.getCone(i) for i in range(cStart,cEnd)]
cells_indicies = np.vstack([np.hstack([plex.getCone(sIndex[k])-vStart
for k in range(len(sIndex))]) for sIndex in sIndicies])
cells = buildCells(coordinates, plex)
vStart, _ = plex.getDepthStratum(0)

ngMesh.Add(ngm.FaceDescriptor(bc=1))
cells = build2DCells(coordinates, sIndicies,
cells_indicies, (cStart, cEnd))
if cells.ndim == 2:
ngMesh.AddElements(dim=2, index=1, data=cells, base=0)
for bcLabel in range(1,plex.getLabelSize(FACE_SETS_LABEL)+1):
if plex.getStratumSize("Face Sets",bcLabel) == 0:
for bcLabel in range(1, plex.getLabelSize(FACE_SETS_LABEL)+1):
if plex.getStratumSize("Face Sets", bcLabel) == 0:
continue
bcIndices = plex.getStratumIS("Face Sets",bcLabel).indices
bcIndices = plex.getStratumIS("Face Sets", bcLabel).indices
for j in bcIndices:
bcIndex = plex.getCone(j)-vStart
bcIndex = plex.getCone(j) - vStart + 1
if len(bcIndex) == 2:
edge = ngm.Element1D([v+1 for v in bcIndex],
edge = ngm.Element1D(list(bcIndex),
index=bcLabel,
edgenr=bcLabel-1)
ngMesh.Add(edge, project_geominfo=geoInfo)


def create3DNetgenMesh(ngMesh, coordinates, plex, geoInfo):
"""
Method used to generate 3D NetgenMeshes

:arg ngMesh: the netgen mesh to be populated
:arg coordinates: vertices coordinates
:arg coordinates: vertices coordinates
:arg plex: PETSc DMPlex
:arg geoInfo: geometric information assosciated with the Netgen mesh

"""
ngMesh.AddPoints(coordinates)
cStart, cEnd = plex.getHeightStratum(0)
vStart, _ = plex.getHeightStratum(3)
# Outside of jitted loop we put all calls to plex
sIndicies = [plex.getCone(i) for i in range(cStart,cEnd)]
f1Indicies = np.array([plex.getCone(s[0]) for s in sIndicies])
f2Indicies = np.array([plex.getCone(s[1]) for s in sIndicies])
fIndicies = np.hstack([f1Indicies,f2Indicies])

cells_indicies = np.vstack([np.hstack([plex.getCone(sIndex[k])-vStart
for k in range(len(sIndex))]) for sIndex in fIndicies])
cells = buildCells(coordinates, plex)
vStart, _ = plex.getDepthStratum(0)

ngMesh.Add(ngm.FaceDescriptor(bc=1))
ngMesh.Add(ngm.FaceDescriptor(bc=plex.getLabelSize(FACE_SETS_LABEL)+1))
cells = build3DCells(coordinates, sIndicies,
cells_indicies, (cStart, cEnd))
if cells.ndim == 2:
ngMesh.AddElements(dim=3, index=plex.getLabelSize(FACE_SETS_LABEL)+1,
data=cells, base=0)
for bcLabel in range(1,plex.getLabelSize(FACE_SETS_LABEL)+1):
data=cells, base=0)
for bcLabel in range(1, plex.getLabelSize(FACE_SETS_LABEL)+1):
faces = []
if plex.getStratumSize("Face Sets",bcLabel) == 0:
if plex.getStratumSize("Face Sets", bcLabel) == 0:
continue
bcIndices = plex.getStratumIS("Face Sets",bcLabel).indices
bcIndices = plex.getStratumIS("Face Sets", bcLabel).indices
for j in bcIndices:
sIndex = plex.getCone(j)
if len(sIndex)==3:
S = list(itertools.chain.from_iterable([
list(plex.getCone(sIndex[k])-vStart) for k in range(len(sIndex))]))
face = list(dict.fromkeys(S))
A = np.array([coordinates[face[1]]-coordinates[face[0]],
coordinates[face[2]]-coordinates[face[1]],
coordinates[face[0]]-coordinates[face[2]]])
if len(sIndex) == 3:
S = dict.fromkeys(itertools.chain.from_iterable(
plex.getCone(p) - vStart
for p in sIndex)
)
face = list(S)
A = coordinates[face][[1, 2, 0]] - coordinates[face]
eig = np.linalg.eig(A)[0].real
if eig[1]*eig[2] > 0:
faces = faces + [face]
else:
faces = faces + [[face[0],face[2],face[1]]]
if eig[1] * eig[2] < 0:
face = [face[0], face[2], face[1]]
faces.append(face)

ngMesh.Add(ngm.FaceDescriptor(bc=bcLabel, surfnr=bcLabel))
ngMesh.AddElements(dim=2, index=bcLabel,
data=np.asarray(faces,dtype=np.int32), base=0,
data=np.asarray(faces, dtype=np.int32), base=0,
project_geometry = geoInfo)


Expand Down Expand Up @@ -240,7 +141,7 @@
:arg plex: the PETSc DMPlex to be converted in NGSolve mesh object

'''
if plex.getDimension() not in [2,3]:
if plex.getDimension() not in {2, 3}:
raise NotImplementedError(f"Not implemented for dimension {plex.getDimension()}.")
nv = plex.getDepthStratum(0)[1] - plex.getDepthStratum(0)[0]
coordinates = plex.getCoordinatesLocal().getArray()
Expand Down Expand Up @@ -291,7 +192,6 @@
T = trim_util(T)

plex = PETSc.DMPlex().createFromCellList(dim, T, V, comm=comm)
plex.setName(self.name)
vStart, _ = plex.getDepthStratum(0)
if surfMesh:
for e in self.ngMesh.Elements1D():
Expand All @@ -304,20 +204,17 @@
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))
self.petscPlex = plex
else:
plex = PETSc.DMPlex().createFromCellList(3,
np.zeros((0, 4), dtype=np.int32),
np.zeros((0, 3), dtype=np.double),
comm=comm)
self.petscPlex = plex
np.zeros((0, 4), dtype=np.int32),
np.zeros((0, 3), dtype=np.double),
comm=comm)
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
T = np.array([list(np.trim_zeros(a, 'b')) for a in T]) - 1
plex = PETSc.DMPlex().createFromCellList(2, T, V, comm=comm)
plex.setName(self.name)
vStart, _ = plex.getDepthStratum(0) # vertices
for e in self.ngMesh.Elements1D():
join = plex.getJoin([vStart+v.nr-1 for v in e.vertices])
Expand All @@ -326,11 +223,10 @@
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))

self.petscPlex = plex
else:
plex = PETSc.DMPlex().createFromCellList(2,
np.zeros((0, 3), dtype=np.int32),
np.zeros((0, 2), dtype=np.double),
comm=comm)
self.petscPlex = plex
np.zeros((0, 3), dtype=np.int32),
np.zeros((0, 2), dtype=np.double),
comm=comm)
plex.setName(self.name)
self.petscPlex = plex
Loading