Skip to content

Commit 2605268

Browse files
UZerbinatiUZerbinati
authored
3D Hierarchy support
3D Hierarchy support and speed up of mesh construction.
1 parent 36719a4 commit 2605268

File tree

5 files changed

+310
-163
lines changed

5 files changed

+310
-163
lines changed

‎Makefile‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8,11 +8,11 @@ else
88
endif
99

1010
lint:
11-
pylint ${PYLINT_FORMAT} --disable=C0412,C0103,C0415,C0321,C3001,E0401,E1101,E0611,R1728,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,R0917,I1101,W0102,W0201,W0212,W0406,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc
12-
pylint ${PYLINT_FORMAT} --disable=C0412,C0103,C0415,C0321,C3001,E0401,E1101,E0611,R1728,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,R0917,I1101,W0102,W0201,W0212,W0406,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc/utils
11+
pylint ${PYLINT_FORMAT} --disable=C0412,C0103,C0415,C0321,C3001,E0401,E1101,E0611,R1728,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,R0917,I1101,W0102,W0201,W0212,W0406,W1203,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc
12+
pylint ${PYLINT_FORMAT} --disable=C0412,C0103,C0415,C0321,C3001,E0401,E1101,E0611,R1728,R1736,R0401,R0801,R0902,R1702,R0913,R0914,R0903,R0205,R0912,R0915,R0917,I1101,W0102,W0201,W0212,W0406,W1203,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase ngsPETSc/utils
1313

1414
lint_test:
15-
pylint ${PYLINT_FORMAT} --disable=C0412,C0103,C0415,C0321,E0401,E1101,E0611,R1728,R1736,R0401,R0914,R0801,R0902,R1702,R0913,R0903,R0205,R0912,R0915,I1101,W0201,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase tests
15+
pylint ${PYLINT_FORMAT} --disable=C0412,C0103,C0415,C0321,E0401,E1101,E0611,R1728,R1736,R0401,R0914,R0801,R0902,R1702,R0913,R0903,R0205,R0912,R0915,I1101,W0201,W1203,C0209 --variable-naming-style=camelCase --class-naming-style=PascalCase --argument-naming-style=camelCase --attr-naming-style=camelCase tests
1616

1717
test:
1818
pytest tests

‎ngsPETSc/plex.py‎

Lines changed: 213 additions & 106 deletions
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,9 @@
44
'''
55
import itertools
66
import numpy as np
7-
from packaging.version import Version
87
from petsc4py import PETSc
9-
108
import netgen.meshing as ngm
9+
from ngsPETSc.utils.utils import trim_util
1110
try:
1211
import ngsolve as ngs
1312
except ImportError:
@@ -16,11 +15,213 @@ class ngs:
1615
class comp:
1716
"dummy class"
1817
Mesh = type(None)
18+
try:
19+
from numba import jit
20+
numba = True
21+
except: # pylint: disable=bare-except
22+
jit = None
23+
numba = False
1924

2025
FACE_SETS_LABEL = "Face Sets"
2126
CELL_SETS_LABEL = "Cell Sets"
2227
EDGE_SETS_LABEL = "Edge Sets"
2328

29+
if numba:
30+
@jit(nopython=True, cache=True)
31+
def build2DCells(coordinates, sIndicies, cell_indicies, start_end):
32+
"""
33+
JIT Method used to construct 2D cells
34+
"""
35+
cStart = start_end[0]
36+
cEnd = start_end[1]
37+
cells = np.zeros((cell_indicies.shape[0], 3))
38+
for i in range(cStart,cEnd):
39+
sIndex = sIndicies[i]
40+
if len(sIndex)==3:
41+
cell = list(set(cell_indicies[i]))
42+
A = np.zeros((2,2))
43+
A[0,0] = (coordinates[cell[1]]-coordinates[cell[0]])[0]
44+
A[0,1] = (coordinates[cell[1]]-coordinates[cell[0]])[1]
45+
A[1,0] = (coordinates[cell[2]]-coordinates[cell[1]])[0]
46+
A[1,1] = (coordinates[cell[2]]-coordinates[cell[1]])[1]
47+
if np.linalg.det(A) > 0:
48+
cells[i] = cell
49+
else:
50+
cells[i] = np.array([cell[0],cell[2],cell[1]])
51+
else:
52+
raise RuntimeError("We only support triangles.")
53+
return cells
54+
55+
@jit(nopython=True, cache=True)
56+
def build3DCells(coordinates, sIndicies, cell_indicies, start_end):
57+
"""
58+
JIT Method used to construct 3D cells
59+
"""
60+
cStart = start_end[0]
61+
cEnd = start_end[1]
62+
cells = np.zeros((cell_indicies.shape[0], 4))
63+
for i in range(cStart,cEnd):
64+
sIndex = sIndicies[i]
65+
if len(sIndex)==4:
66+
cell = list(set(cell_indicies[i]))
67+
A = np.zeros((3,3))
68+
A[0,0] = (coordinates[cell[1]]-coordinates[cell[0]])[0]
69+
A[0,1] = (coordinates[cell[1]]-coordinates[cell[0]])[1]
70+
A[0,2] = (coordinates[cell[1]]-coordinates[cell[0]])[2]
71+
A[1,0] = (coordinates[cell[2]]-coordinates[cell[1]])[0]
72+
A[1,1] = (coordinates[cell[2]]-coordinates[cell[1]])[1]
73+
A[1,2] = (coordinates[cell[2]]-coordinates[cell[1]])[2]
74+
A[2,0] = (coordinates[cell[3]]-coordinates[cell[2]])[0]
75+
A[2,1] = (coordinates[cell[3]]-coordinates[cell[2]])[1]
76+
A[2,2] = (coordinates[cell[3]]-coordinates[cell[2]])[2]
77+
if np.linalg.det(A) > 0:
78+
cells[i] = cell
79+
else:
80+
cells[i] = np.array([cell[0],cell[1],cell[2], cell[2]])
81+
else:
82+
raise RuntimeError("We only support tets.")
83+
return cells
84+
else:
85+
def build2DCells(coordinates, sIndicies, cell_indicies, start_end):
86+
"""
87+
Method used to construct 2D cells
88+
"""
89+
cStart = start_end[0]
90+
cEnd = start_end[1]
91+
cells = np.zeros((cell_indicies.shape[0], 3))
92+
for i in range(cStart,cEnd):
93+
sIndex = sIndicies[i]
94+
if len(sIndex)==3:
95+
cell = list(set(cell_indicies[i]))
96+
A = np.zeros((2,2))
97+
A[0,0] = (coordinates[cell[1]]-coordinates[cell[0]])[0]
98+
A[0,1] = (coordinates[cell[1]]-coordinates[cell[0]])[1]
99+
A[1,0] = (coordinates[cell[2]]-coordinates[cell[1]])[0]
100+
A[1,1] = (coordinates[cell[2]]-coordinates[cell[1]])[1]
101+
if np.linalg.det(A) > 0:
102+
cells[i] = cell
103+
else:
104+
cells[i] = np.array([cell[0],cell[2],cell[1]])
105+
else:
106+
raise RuntimeError("We only support triangles.")
107+
return cells
108+
109+
def build3DCells(coordinates, sIndicies, cell_indicies, start_end):
110+
"""
111+
Method used to construct 3D cells
112+
"""
113+
cStart = start_end[0]
114+
cEnd = start_end[1]
115+
cells = np.zeros((cell_indicies.shape[0], 4))
116+
for i in range(cStart,cEnd):
117+
sIndex = sIndicies[i]
118+
if len(sIndex)==4:
119+
cell = list(set(cell_indicies[i]))
120+
A = np.zeros((3,3))
121+
A[0,0] = (coordinates[cell[1]]-coordinates[cell[0]])[0]
122+
A[0,1] = (coordinates[cell[1]]-coordinates[cell[0]])[1]
123+
A[0,2] = (coordinates[cell[1]]-coordinates[cell[0]])[2]
124+
A[1,0] = (coordinates[cell[2]]-coordinates[cell[1]])[0]
125+
A[1,1] = (coordinates[cell[2]]-coordinates[cell[1]])[1]
126+
A[1,2] = (coordinates[cell[2]]-coordinates[cell[1]])[2]
127+
A[2,0] = (coordinates[cell[3]]-coordinates[cell[2]])[0]
128+
A[2,1] = (coordinates[cell[3]]-coordinates[cell[2]])[1]
129+
A[2,2] = (coordinates[cell[3]]-coordinates[cell[2]])[2]
130+
if np.linalg.det(A) > 0:
131+
cells[i] = cell
132+
else:
133+
cells[i] = np.array([cell[0],cell[1],cell[3], cell[2]])
134+
else:
135+
raise RuntimeError("We only support tets.")
136+
return cells
137+
138+
def create2DNetgenMesh(ngMesh, coordinates, plex, geoInfo):
139+
"""
140+
Method used to generate 2D NetgenMeshes
141+
142+
:arg ngMesh: the netgen mesh to be populated
143+
:arg coordinates: vertices coordinates
144+
:arg plex: PETSc DMPlex
145+
:arg geoInfo: geometric information assosciated with the Netgen mesh
146+
147+
"""
148+
ngMesh.AddPoints(coordinates)
149+
cStart,cEnd = plex.getHeightStratum(0)
150+
vStart, _ = plex.getHeightStratum(2)
151+
# Outside of jitted loop we put all calls to plex
152+
sIndicies = [plex.getCone(i) for i in range(cStart,cEnd)]
153+
cells_indicies = np.vstack([np.hstack([plex.getCone(sIndex[k])-vStart
154+
for k in range(len(sIndex))]) for sIndex in sIndicies])
155+
ngMesh.Add(ngm.FaceDescriptor(bc=1))
156+
cells = build2DCells(coordinates, sIndicies,
157+
cells_indicies, (cStart, cEnd))
158+
if cells.ndim == 2:
159+
ngMesh.AddElements(dim=2, index=1, data=cells, base=0)
160+
for bcLabel in range(1,plex.getLabelSize(FACE_SETS_LABEL)+1):
161+
if plex.getStratumSize("Face Sets",bcLabel) == 0:
162+
continue
163+
bcIndices = plex.getStratumIS("Face Sets",bcLabel).indices
164+
for j in bcIndices:
165+
bcIndex = plex.getCone(j)-vStart
166+
if len(bcIndex) == 2:
167+
edge = ngm.Element1D([v+1 for v in bcIndex],
168+
index=bcLabel,
169+
edgenr=bcLabel-1)
170+
ngMesh.Add(edge, project_geominfo=geoInfo)
171+
172+
def create3DNetgenMesh(ngMesh, coordinates, plex, geoInfo):
173+
"""
174+
Method used to generate 3D NetgenMeshes
175+
176+
:arg ngMesh: the netgen mesh to be populated
177+
:arg coordinates: vertices coordinates
178+
:arg plex: PETSc DMPlex
179+
:arg geoInfo: geometric information assosciated with the Netgen mesh
180+
181+
"""
182+
ngMesh.AddPoints(coordinates)
183+
cStart, cEnd = plex.getHeightStratum(0)
184+
vStart, _ = plex.getHeightStratum(3)
185+
# Outside of jitted loop we put all calls to plex
186+
sIndicies = [plex.getCone(i) for i in range(cStart,cEnd)]
187+
f1Indicies = np.array([plex.getCone(s[0]) for s in sIndicies])
188+
f2Indicies = np.array([plex.getCone(s[1]) for s in sIndicies])
189+
fIndicies = np.hstack([f1Indicies,f2Indicies])
190+
191+
cells_indicies = np.vstack([np.hstack([plex.getCone(sIndex[k])-vStart
192+
for k in range(len(sIndex))]) for sIndex in fIndicies])
193+
ngMesh.Add(ngm.FaceDescriptor(bc=1))
194+
ngMesh.Add(ngm.FaceDescriptor(bc=plex.getLabelSize(FACE_SETS_LABEL)+1))
195+
cells = build3DCells(coordinates, sIndicies,
196+
cells_indicies, (cStart, cEnd))
197+
if cells.ndim == 2:
198+
ngMesh.AddElements(dim=3, index=plex.getLabelSize(FACE_SETS_LABEL)+1,
199+
data=cells, base=0)
200+
for bcLabel in range(1,plex.getLabelSize(FACE_SETS_LABEL)+1):
201+
faces = []
202+
if plex.getStratumSize("Face Sets",bcLabel) == 0:
203+
continue
204+
bcIndices = plex.getStratumIS("Face Sets",bcLabel).indices
205+
for j in bcIndices:
206+
sIndex = plex.getCone(j)
207+
if len(sIndex)==3:
208+
S = list(itertools.chain.from_iterable([
209+
list(plex.getCone(sIndex[k])-vStart) for k in range(len(sIndex))]))
210+
face = list(dict.fromkeys(S))
211+
A = np.array([coordinates[face[1]]-coordinates[face[0]],
212+
coordinates[face[2]]-coordinates[face[1]],
213+
coordinates[face[0]]-coordinates[face[2]]])
214+
eig = np.linalg.eig(A)[0].real
215+
if eig[1]*eig[2] > 0:
216+
faces = faces + [face]
217+
else:
218+
faces = faces + [[face[0],face[2],face[1]]]
219+
ngMesh.Add(ngm.FaceDescriptor(bc=bcLabel, surfnr=bcLabel))
220+
ngMesh.AddElements(dim=2, index=bcLabel,
221+
data=np.asarray(faces,dtype=np.int32), base=0,
222+
project_geometry = geoInfo)
223+
224+
24225
class MeshMapping:
25226
'''
26227
This class creates a mapping between Netgen/NGSolve meshes and PETSc DMPlex
@@ -31,9 +232,10 @@ class MeshMapping:
31232
32233
'''
33234

34-
def __init__(self, mesh=None, comm=None, name="Default"):
235+
def __init__(self, mesh=None, comm=None, geo=None, name="Default"):
35236
self.name = name
36237
self.comm = comm if comm is not None else PETSc.COMM_WORLD
238+
self.geo = geo
37239
if isinstance(mesh,(ngs.comp.Mesh,ngm.Mesh)):
38240
self.createPETScDMPlex(mesh)
39241
elif isinstance(mesh,PETSc.DMPlex):
@@ -58,103 +260,16 @@ def createNGSMesh(self, plex):
58260
self.petscPlex = plex
59261
ngMesh = ngm.Mesh(dim=plex.getCoordinateDim())
60262
self.ngMesh = ngMesh
263+
self.geoInfo = False
264+
if self.geo:
265+
self.ngMesh.SetGeometry(self.geo)
266+
self.geoInfo = True
61267

62268
coordinates = coordinates.reshape(-1,plex.getDimension())
63269
if plex.getDimension() == 2:
64-
self.ngMesh.AddPoints(coordinates)
65-
cStart,cEnd = plex.getHeightStratum(0)
66-
vStart, _ = plex.getHeightStratum(2)
67-
cells = []
68-
for i in range(cStart,cEnd):
69-
sIndex = plex.getCone(i)
70-
if len(sIndex)==3:
71-
S = list(itertools.chain.from_iterable([
72-
list(plex.getCone(sIndex[k])-vStart) for k in range(len(sIndex))]))
73-
cell = list(dict.fromkeys(S))
74-
A = np.array([coordinates[cell[1]]-coordinates[cell[0]],
75-
coordinates[cell[2]]-coordinates[cell[1]]])
76-
if np.linalg.det(A) > 0:
77-
cells = cells + [cell]
78-
else:
79-
cells = cells + [[cell[0],cell[2],cell[1]]]
80-
else:
81-
S = [plex.getCone(sIndex[k])-vStart for k in list(range(len(sIndex)))]
82-
k = 0
83-
cell = []
84-
while k < len(sIndex):
85-
if S[k][0] == S[k-1][1]:
86-
cell = cell+[S[k][0]]
87-
else:
88-
if k == 0 and S[k][0] not in list(S[k+1]):
89-
S[-1] = np.array([S[-1][1],S[-1][0]])
90-
else:
91-
S[k] = np.array([S[k][1],S[k][0]])
92-
cell = cell+[S[k][0]]
93-
k = k+1
94-
cells = cells + [cell]
95-
self.ngMesh.Add(ngm.FaceDescriptor(bc=1))
96-
cells = np.asarray(cells, dtype=np.int32)
97-
if cells.ndim == 2:
98-
self.ngMesh.AddElements(dim=2, index=1, data=cells, base=0)
99-
for bcLabel in range(1,plex.getLabelSize(FACE_SETS_LABEL)+1):
100-
if plex.getStratumSize("Face Sets",bcLabel) == 0:
101-
continue
102-
bcIndices = plex.getStratumIS("Face Sets",bcLabel).indices
103-
for j in bcIndices:
104-
bcIndex = plex.getCone(j)-vStart
105-
if len(bcIndex) == 2:
106-
edge = ngm.Element1D([v+1 for v in bcIndex],index=bcLabel)
107-
self.ngMesh.Add(edge)
270+
create2DNetgenMesh(self.ngMesh, coordinates, plex, self.geoInfo)
108271
elif plex.getDimension() == 3:
109-
self.ngMesh.AddPoints(coordinates)
110-
cStart, cEnd = plex.getHeightStratum(0)
111-
vStart, _ = plex.getHeightStratum(3)
112-
cells = []
113-
for i in range(cStart,cEnd):
114-
fIndex = plex.getCone(i)
115-
f1 = plex.getCone(fIndex[0])
116-
f2 = plex.getCone(fIndex[1])
117-
sIndex = [f1[0],f1[1],f1[2],f2[0],f2[1],f2[2]]
118-
if len(fIndex) == 4:
119-
S = list(itertools.chain.from_iterable([
120-
list(plex.getCone(sIndex[k])-vStart) for k in list(range(len(sIndex)))]))
121-
cell = list(dict.fromkeys(S))
122-
A = np.array([coordinates[cell[1]]-coordinates[cell[0]],
123-
coordinates[cell[2]]-coordinates[cell[1]],
124-
coordinates[cell[3]]-coordinates[cell[2]]])
125-
if np.linalg.det(A) > 0:
126-
cells = cells + [cell]
127-
else:
128-
cells = cells + [[cell[0],cell[1],cell[3],cell[2]]]
129-
#fd = self.ngmesh.Add(ngm.FaceDescriptor(bc=plex.getLabelSize(FACE_SETS_LABEL)+1))
130-
self.ngMesh.Add(ngm.FaceDescriptor(bc=plex.getLabelSize(FACE_SETS_LABEL)+1))
131-
cells = np.asarray(cells, dtype=np.int32)
132-
if cells.ndim == 2:
133-
self.ngMesh.AddElements(dim=3, index=plex.getLabelSize(FACE_SETS_LABEL)+1,
134-
data=cells, base=0)
135-
for bcLabel in range(1,plex.getLabelSize(FACE_SETS_LABEL)+1):
136-
faces = []
137-
if plex.getStratumSize("Face Sets",bcLabel) == 0:
138-
continue
139-
bcIndices = plex.getStratumIS("Face Sets",bcLabel).indices
140-
for j in bcIndices:
141-
sIndex = plex.getCone(j)
142-
if len(sIndex)==3:
143-
S = list(itertools.chain.from_iterable([
144-
list(plex.getCone(sIndex[k])-vStart) for k in range(len(sIndex))]))
145-
face = list(dict.fromkeys(S))
146-
A = np.array([coordinates[face[1]]-coordinates[face[0]],
147-
coordinates[face[2]]-coordinates[face[1]],
148-
coordinates[face[0]]-coordinates[face[2]]])
149-
eig = np.linalg.eig(A)[0].real
150-
if eig[1]*eig[2] > 0:
151-
faces = faces + [face]
152-
else:
153-
faces = faces + [[face[0],face[2],face[1]]]
154-
#fd = self.ngmesh.Add(ngm.FaceDescriptor(bc=bcLabel))
155-
self.ngMesh.Add(ngm.FaceDescriptor(bc=bcLabel))
156-
self.ngMesh.AddElements(dim=2, index=bcLabel,
157-
data=np.asarray(faces,dtype=np.int32), base=0)
272+
create3DNetgenMesh(self.ngMesh, coordinates, plex, self.geoInfo)
158273
else:
159274
raise NotImplementedError("No implementation for dimension greater than 3.")
160275

@@ -170,6 +285,7 @@ def createPETScDMPlex(self, mesh):
170285
else:
171286
self.ngMesh = mesh
172287
comm = self.comm
288+
self.geo = self.ngMesh.GetGeometry()
173289
if self.ngMesh.dim == 3:
174290
if comm.rank == 0:
175291
V = self.ngMesh.Coordinates()
@@ -180,16 +296,7 @@ def createPETScDMPlex(self, mesh):
180296
surfMesh, dim = True, 2
181297
T = self.ngMesh.Elements2D().NumPy()["nodes"]
182298

183-
if Version(np.__version__) >= Version("2.2"):
184-
T = np.trim_zeros(T, "b", axis=1).astype(np.int32) - 1
185-
else:
186-
T = (
187-
np.array(
188-
[list(np.trim_zeros(a, "b")) for a in list(T)],
189-
dtype=np.int32,
190-
)
191-
- 1
192-
)
299+
T = trim_util(T)
193300

194301
plex = PETSc.DMPlex().createFromCellList(dim, T, V, comm=comm)
195302
plex.setName(self.name)

0 commit comments

Comments
 (0)