Skip to content

Commit 6f507af

Browse files
RomainBavilleRomainBaville
authored
feat: Add the possibility to map cells for same type meshes but with differents cell dimension (#191)
* Deals with mesh of same type but with cells with different dimensions for the element map functions
1 parent b5ce4f4 commit 6f507af

File tree

4 files changed

+198
-77
lines changed

4 files changed

+198
-77
lines changed

‎geos-mesh/src/geos/mesh/utils/arrayHelpers.py‎

Lines changed: 28 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,8 @@
1111
from vtkmodules.util.numpy_support import vtk_to_numpy
1212
from vtkmodules.vtkCommonCore import vtkDataArray, vtkPoints
1313
from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkFieldData, vtkMultiBlockDataSet, vtkDataSet,
14-
vtkCompositeDataSet, vtkDataObject, vtkPointData, vtkCellData, vtkPolyData )
14+
vtkCompositeDataSet, vtkDataObject, vtkPointData, vtkCellData, vtkPolyData,
15+
vtkCell )
1516
from vtkmodules.vtkFiltersCore import vtkCellCenters
1617
from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten
1718

@@ -275,42 +276,49 @@ def UpdateDictElementMappingFromDataSetToDataSet(
275276
for idElementTo in range( nbElementsTo ):
276277
# Test if the element of the final mesh is already mapped.
277278
if -1 in elementMap[ flatIdDataSetTo ][ idElementTo ]:
278-
coordElementTo: list[ tuple[ float, ...] ] = []
279+
typeElemTo: int
280+
coordElementTo: set[ tuple[ float, ...] ] = set()
279281
if points:
280-
coordElementTo.append( dataSetTo.GetPoint( idElementTo ) )
282+
typeElemTo = 0
283+
coordElementTo.add( dataSetTo.GetPoint( idElementTo ) )
281284
else:
285+
cellTo: vtkCell = dataSetTo.GetCell( idElementTo )
286+
typeElemTo = cellTo.GetCellType()
282287
# Get the coordinates of each points of the cell.
283-
nbPointsTo: int = dataSetTo.GetCell( idElementTo ).GetNumberOfPoints()
284-
cellPointsTo: vtkPoints = dataSetTo.GetCell( idElementTo ).GetPoints()
288+
nbPointsTo: int = cellTo.GetNumberOfPoints()
289+
cellPointsTo: vtkPoints = cellTo.GetPoints()
285290
for idPointTo in range( nbPointsTo ):
286-
coordElementTo.append( cellPointsTo.GetPoint( idPointTo ) )
291+
coordElementTo.add( cellPointsTo.GetPoint( idPointTo ) )
287292

288293
idElementFrom: int = 0
289294
ElementFromFund: bool = False
290295
while idElementFrom < nbElementsFrom and not ElementFromFund:
291296
# Test if the element of the source mesh is already mapped.
292297
if idElementFrom not in idElementsFromFund:
293-
coordElementFrom: list[ tuple[ float, ...] ] = []
298+
typeElemFrom: int
299+
coordElementFrom: set[ tuple[ float, ...] ] = set()
294300
if points:
295-
coordElementFrom.append( dataSetFrom.GetPoint( idElementFrom ) )
301+
typeElemFrom = 0
302+
coordElementFrom.add( dataSetFrom.GetPoint( idElementFrom ) )
296303
else:
297-
# Get the coordinates of each points of the cell.
298-
nbPointsFrom: int = dataSetFrom.GetCell( idElementFrom ).GetNumberOfPoints()
299-
cellPointsFrom: vtkPoints = dataSetFrom.GetCell( idElementFrom ).GetPoints()
304+
cellFrom: vtkCell = dataSetFrom.GetCell( idElementFrom )
305+
typeElemFrom = cellFrom.GetCellType()
306+
# Get the coordinates of each points of the face.
307+
nbPointsFrom: int = cellFrom.GetNumberOfPoints()
308+
cellPointsFrom: vtkPoints = cellFrom.GetPoints()
300309
for idPointFrom in range( nbPointsFrom ):
301-
coordElementFrom.append( cellPointsFrom.GetPoint( idPointFrom ) )
310+
coordElementFrom.add( cellPointsFrom.GetPoint( idPointFrom ) )
302311

303312
pointShared: bool = True
304-
if dataSetTo.GetClassName() == dataSetFrom.GetClassName():
313+
if typeElemTo == typeElemFrom:
305314
if coordElementTo != coordElementFrom:
306315
pointShared = False
307-
elif isinstance( dataSetTo, vtkPolyData ):
308-
for coordPointsTo in coordElementTo:
309-
if coordPointsTo not in coordElementFrom:
316+
else:
317+
if nbPointsTo < nbPointsFrom:
318+
if not coordElementTo.issubset( coordElementFrom ):
310319
pointShared = False
311-
elif isinstance( dataSetFrom, vtkPolyData ):
312-
for coordPointsFrom in coordElementFrom:
313-
if coordPointsFrom not in coordElementTo:
320+
else:
321+
if not coordElementTo.issuperset( coordElementFrom ):
314322
pointShared = False
315323

316324
if pointShared:
@@ -319,6 +327,7 @@ def UpdateDictElementMappingFromDataSetToDataSet(
319327
idElementsFromFund.append( idElementFrom )
320328

321329
idElementFrom += 1
330+
return
322331

323332

324333
def hasArray( mesh: vtkUnstructuredGrid, arrayNames: list[ str ] ) -> bool:

‎geos-mesh/tests/conftest.py‎

Lines changed: 98 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -182,6 +182,8 @@ def _get_dataset( datasetType: str ) -> Union[ vtkMultiBlockDataSet, vtkPolyData
182182
vtkFilename = "data/meshGeosExtractBlockTmp.vtm"
183183
elif datasetType == "well":
184184
vtkFilename = "data/well.vtu"
185+
elif datasetType == "emptyWell":
186+
vtkFilename = "data/well_empty.vtu"
185187
datapath: str = os.path.join( os.path.dirname( os.path.realpath( __file__ ) ), vtkFilename )
186188
reader.SetFileName( datapath )
187189
reader.Update()
@@ -210,62 +212,110 @@ def _get_elementMap( meshFromName: str, meshToName: str, points: bool ) -> Dict[
210212
Returns:
211213
elementMap (Dict[int, npt.NDArray[np.int64]]): The element mapping dictionary.
212214
"""
215+
sharedCells2D3DId: npt.NDArray[ np.int64 ] = np.array(
216+
[ [ 0, 0 ], [ 1, 1 ], [ 2, 2 ], [ 3, 3 ], [ 4, 4 ], [ 5, 5 ], [ 6, 6 ], [ 7, 7 ], [ 8, 8 ], [ 9, 9 ],
217+
[ 10, 10 ], [ 11, 11 ], [ 12, 12 ], [ 13, 13 ], [ 14, 14 ], [ 15, 15 ], [ 16, 16 ], [ 17, 17 ],
218+
[ 18, 18 ], [ 19, 19 ], [ 20, 20 ], [ 21, 21 ], [ 22, 22 ], [ 23, 23 ], [ 24, 48 ], [ 25, 50 ],
219+
[ 26, 51 ], [ 27, 54 ], [ 28, 56 ], [ 29, 57 ], [ 30, 58 ], [ 31, 59 ], [ 32, 60 ], [ 33, 61 ],
220+
[ 34, 62 ], [ 35, 63 ], [ 36, 64 ], [ 37, 65 ], [ 38, 66 ], [ 39, 67 ], [ 40, 68 ], [ 41, 69 ],
221+
[ 42, 70 ], [ 43, 71 ], [ 44, 72 ], [ 45, 73 ], [ 46, 74 ], [ 47, 75 ], [ 48, 76 ], [ 49, 77 ],
222+
[ 50, 78 ], [ 51, 79 ], [ 52, 580 ], [ 53, 581 ], [ 54, 582 ], [ 55, 583 ], [ 56, 584 ], [ 57, 585 ],
223+
[ 58, 586 ], [ 59, 587 ], [ 60, 588 ], [ 61, 589 ], [ 62, 590 ], [ 63, 591 ], [ 64, 592 ], [ 65, 593 ],
224+
[ 66, 594 ], [ 67, 595 ], [ 68, 596 ], [ 69, 597 ], [ 70, 598 ], [ 71, 599 ], [ 72, 600 ], [ 73, 601 ],
225+
[ 74, 602 ], [ 75, 603 ], [ 76, 628 ], [ 77, 630 ], [ 78, 631 ], [ 79, 634 ], [ 80, 636 ], [ 81, 637 ],
226+
[ 82, 638 ], [ 83, 639 ], [ 84, 640 ], [ 85, 641 ], [ 86, 642 ], [ 87, 643 ], [ 88, 644 ], [ 89, 645 ],
227+
[ 90, 646 ], [ 91, 647 ], [ 92, 648 ], [ 93, 649 ], [ 94, 650 ], [ 95, 651 ], [ 96, 652 ], [ 97, 653 ],
228+
[ 98, 654 ], [ 99, 655 ], [ 100, 656 ], [ 101, 657 ], [ 102, 658 ], [ 103, 659 ], [ 104, 1160 ],
229+
[ 105, 1161 ], [ 106, 1162 ], [ 107, 1163 ], [ 108, 1164 ], [ 109, 1165 ], [ 110, 1166 ], [ 111, 1167 ],
230+
[ 112, 1168 ], [ 113, 1169 ], [ 114, 1170 ], [ 115, 1171 ], [ 116, 1172 ], [ 117, 1173 ], [ 118, 1174 ],
231+
[ 119, 1175 ], [ 120, 1176 ], [ 121, 1177 ], [ 122, 1178 ], [ 123, 1179 ], [ 124, 1180 ], [ 125, 1181 ],
232+
[ 126, 1182 ], [ 127, 1183 ], [ 128, 1208 ], [ 129, 1210 ], [ 130, 1211 ], [ 131, 1214 ], [ 132, 1216 ],
233+
[ 133, 1217 ], [ 134, 1218 ], [ 135, 1219 ], [ 136, 1220 ], [ 137, 1221 ], [ 138, 1222 ], [ 139, 1223 ],
234+
[ 140, 1224 ], [ 141, 1225 ], [ 142, 1226 ], [ 143, 1227 ], [ 144, 1228 ], [ 145, 1229 ], [ 146, 1230 ],
235+
[ 147, 1231 ], [ 148, 1232 ], [ 149, 1233 ], [ 150, 1234 ], [ 151, 1235 ], [ 152, 1236 ], [ 153, 1237 ],
236+
[ 154, 1238 ], [ 155, 1239 ] ],
237+
dtype=np.int64,
238+
)
239+
sharedPoints1D2DId: npt.NDArray[ np.int64 ] = np.array( [ [ 0, 26 ] ], dtype=np.int64 )
240+
sharedPoints1D3DId: npt.NDArray[ np.int64 ] = np.array( [ [ 0, 475 ] ], dtype=np.int64 )
213241
elementMap: Dict[ int, npt.NDArray[ np.int64 ] ] = {}
214-
nbElements: Tuple[ int, int ] = ( 4092, 212 ) if points else ( 1740, 156 )
215-
if meshFromName == "multiblock":
216-
if meshToName == "emptymultiblock":
217-
elementMap[ 1 ] = np.array( [ [ 1, element ] for element in range( nbElements[ 0 ] ) ] )
218-
elementMap[ 3 ] = np.array( [ [ 3, element ] for element in range( nbElements[ 1 ] ) ] )
219-
elif meshToName == "emptyFracture":
220-
elementMap[ 0 ] = np.array( [ [ 3, element ] for element in range( nbElements[ 1 ] ) ] )
221-
elif meshFromName == "dataset":
222-
if meshToName == "emptydataset":
223-
elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 0 ] ) ] )
224-
elif meshToName == "emptyFracture":
242+
nbElements: Tuple[ int, int, int ] = ( 4092, 212, 11 ) if points else ( 1740, 156, 10 )
243+
if meshFromName == "well":
244+
if meshToName == "emptyWell":
245+
elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 2 ] ) ] )
246+
elif meshToName == "emptyFracture" or meshToName == "emptypolydata":
225247
elementMap[ 0 ] = np.full( ( nbElements[ 1 ], 2 ), -1, np.int64 )
226-
elif meshToName == "emptypolydata":
227-
elementMap[ 0 ] = np.array( [ [ 0, 0 ], [ 0, 1 ], [ 0, 2 ], [ 0, 3 ], [ 0, 4 ], [ 0, 5 ], [ 0, 6 ],
228-
[ 0, 7 ], [ 0, 8 ], [ 0, 9 ], [ 0, 10 ], [ 0, 11 ], [ 0, 12 ], [ 0, 13 ],
229-
[ 0, 14 ], [ 0, 15 ], [ 0, 16 ], [ 0, 17 ], [ 0, 18 ], [ 0,
230-
19 ], [ 0, 20 ],
231-
[ 0, 21 ], [ 0, 22 ], [ 0, 23 ], [ 0, 48 ], [ 0, 50 ], [ 0,
232-
51 ], [ 0, 54 ],
233-
[ 0, 56 ], [ 0, 57 ], [ 0, 58 ], [ 0, 59 ], [ 0, 60 ], [ 0,
234-
61 ], [ 0, 62 ],
235-
[ 0, 63 ], [ 0, 64 ], [ 0, 65 ], [ 0, 66 ], [ 0, 67 ], [ 0,
236-
68 ], [ 0, 69 ],
237-
[ 0, 70 ], [ 0, 71 ], [ 0, 72 ], [ 0, 73 ], [ 0, 74 ],
238-
[ 0, 75 ], [ 0, 76 ], [ 0, 77 ], [ 0, 78 ], [ 0, 79 ], [ 0, 580 ],
239-
[ 0, 581 ], [ 0, 582 ], [ 0, 583 ], [ 0, 584 ], [ 0, 585 ], [ 0, 586 ],
240-
[ 0, 587 ], [ 0, 588 ], [ 0, 589 ], [ 0, 590 ], [ 0, 591 ], [ 0, 592 ],
241-
[ 0, 593 ], [ 0, 594 ], [ 0, 595 ], [ 0, 596 ], [ 0, 597 ], [ 0, 598 ],
242-
[ 0, 599 ], [ 0, 600 ], [ 0, 601 ], [ 0, 602 ], [ 0, 603 ], [ 0, 628 ],
243-
[ 0, 630 ], [ 0, 631 ], [ 0, 634 ], [ 0, 636 ], [ 0, 637 ], [ 0, 638 ],
244-
[ 0, 639 ], [ 0, 640 ], [ 0, 641 ], [ 0, 642 ], [ 0, 643 ], [ 0, 644 ],
245-
[ 0, 645 ], [ 0, 646 ], [ 0, 647 ], [ 0, 648 ], [ 0, 649 ], [ 0, 650 ],
246-
[ 0, 651 ], [ 0, 652 ], [ 0, 653 ], [ 0, 654 ], [ 0, 655 ], [ 0, 656 ],
247-
[ 0, 657 ], [ 0, 658 ], [ 0, 659 ], [ 0, 1160 ], [ 0, 1161 ], [ 0, 1162 ],
248-
[ 0, 1163 ], [ 0, 1164 ], [ 0, 1165 ], [ 0, 1166 ], [ 0, 1167 ],
249-
[ 0, 1168 ], [ 0, 1169 ], [ 0, 1170 ], [ 0, 1171 ], [ 0, 1172 ],
250-
[ 0, 1173 ], [ 0, 1174 ], [ 0, 1175 ], [ 0, 1176 ], [ 0, 1177 ],
251-
[ 0, 1178 ], [ 0, 1179 ], [ 0, 1180 ], [ 0, 1181 ], [ 0, 1182 ],
252-
[ 0, 1183 ], [ 0, 1208 ], [ 0, 1210 ], [ 0, 1211 ], [ 0, 1214 ],
253-
[ 0, 1216 ], [ 0, 1217 ], [ 0, 1218 ], [ 0, 1219 ], [ 0, 1220 ],
254-
[ 0, 1221 ], [ 0, 1222 ], [ 0, 1223 ], [ 0, 1224 ], [ 0, 1225 ],
255-
[ 0, 1226 ], [ 0, 1227 ], [ 0, 1228 ], [ 0, 1229 ], [ 0, 1230 ],
256-
[ 0, 1231 ], [ 0, 1232 ], [ 0, 1233 ], [ 0, 1234 ], [ 0, 1235 ],
257-
[ 0, 1236 ], [ 0, 1237 ], [ 0, 1238 ], [ 0, 1239 ] ] )
248+
for id2DElem in range( nbElements[ 1 ] ):
249+
for sharedPoint1D2DId in sharedPoints1D2DId:
250+
if id2DElem == sharedPoint1D2DId[ 1 ]:
251+
elementMap[ 0 ][ id2DElem ] = [ 0, sharedPoint1D2DId[ 0 ] ]
252+
elif meshToName == "emptydataset":
253+
elementMap[ 0 ] = np.full( ( nbElements[ 0 ], 2 ), -1, np.int64 )
254+
for id3DElem in range( nbElements[ 0 ] ):
255+
for sharedPoint1D3DId in sharedPoints1D3DId:
256+
if id3DElem == sharedPoint1D3DId[ 1 ]:
257+
elementMap[ 0 ][ id3DElem ] = [ 0, sharedPoint1D3DId[ 0 ] ]
258258
elif meshToName == "emptymultiblock":
259-
elementMap[ 1 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 0 ] ) ] )
259+
elementMap[ 1 ] = np.full( ( nbElements[ 0 ], 2 ), -1, np.int64 )
260+
for id3DElem in range( nbElements[ 0 ] ):
261+
for sharedPoint1D3DId in sharedPoints1D3DId:
262+
if id3DElem == sharedPoint1D3DId[ 1 ]:
263+
elementMap[ 1 ][ id3DElem ] = [ 0, sharedPoint1D3DId[ 0 ] ]
260264
elementMap[ 3 ] = np.full( ( nbElements[ 1 ], 2 ), -1, np.int64 )
261-
elif meshFromName == "fracture":
262-
if meshToName == "emptyFracture":
265+
for id2DElem in range( nbElements[ 1 ] ):
266+
for sharedPoint1D2DId in sharedPoints1D2DId:
267+
if id2DElem == sharedPoint1D2DId[ 1 ]:
268+
elementMap[ 3 ][ id2DElem ] = [ 0, sharedPoint1D2DId[ 0 ] ]
269+
elif meshFromName == "fracture" or meshFromName == "polydata":
270+
if meshToName == "emptyFracture" or meshToName == "emptypolydata":
263271
elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 1 ] ) ] )
272+
elif meshToName == "emptyWell":
273+
elementMap[ 0 ] = np.full( ( nbElements[ 2 ], 2 ), -1, np.int64 )
274+
for id1DElem in range( nbElements[ 2 ] ):
275+
for sharedPoint1D2DId in sharedPoints1D2DId:
276+
if id1DElem == sharedPoint1D2DId[ 0 ]:
277+
elementMap[ 0 ][ id1DElem ] = [ 0, sharedPoint1D2DId[ 1 ] ]
278+
elif meshToName == "emptydataset":
279+
elementMap[ 0 ] = np.full( ( nbElements[ 0 ], 2 ), -1, np.int64 )
280+
for id3DElem in range( nbElements[ 0 ] ):
281+
for sharedCell2D3DId in sharedCells2D3DId:
282+
if id3DElem == sharedCell2D3DId[ 1 ]:
283+
elementMap[ 0 ][ id3DElem ] = [ 0, sharedCell2D3DId[ 0 ] ]
264284
elif meshToName == "emptymultiblock":
265285
elementMap[ 1 ] = np.full( ( nbElements[ 0 ], 2 ), -1, np.int64 )
286+
for id3DElem in range( nbElements[ 0 ] ):
287+
for sharedCell2D3DId in sharedCells2D3DId:
288+
if id3DElem == sharedCell2D3DId[ 1 ]:
289+
elementMap[ 1 ][ id3DElem ] = [ 0, sharedCell2D3DId[ 0 ] ]
266290
elementMap[ 3 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 1 ] ) ] )
267-
elif meshFromName == "polydata" and meshToName == "emptypolydata":
268-
elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 1 ] ) ] )
291+
elif meshFromName == "dataset":
292+
if meshToName == "emptydataset":
293+
elementMap[ 0 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 0 ] ) ] )
294+
elif meshToName == "emptyWell":
295+
elementMap[ 0 ] = np.full( ( nbElements[ 2 ], 2 ), -1, np.int64 )
296+
for id1DElem in range( nbElements[ 2 ] ):
297+
for sharedPoint1D3DId in sharedPoints1D3DId:
298+
if id1DElem == sharedPoint1D3DId[ 0 ]:
299+
elementMap[ 0 ][ id1DElem ] = [ 0, sharedPoint1D3DId[ 1 ] ]
300+
elif meshToName == "emptypolydata" or meshToName == "emptyFracture":
301+
elementMap[ 0 ] = np.array( [ [ 0, element ] for element in sharedCells2D3DId[ :, 1 ] ] )
302+
elif meshToName == "emptymultiblock":
303+
elementMap[ 1 ] = np.array( [ [ 0, element ] for element in range( nbElements[ 0 ] ) ] )
304+
elementMap[ 3 ] = np.array( [ [ 0, element ] for element in sharedCells2D3DId[ :, 1 ] ] )
305+
elif meshFromName == "multiblock":
306+
if meshToName == "emptymultiblock":
307+
elementMap[ 1 ] = np.array( [ [ 1, element ] for element in range( nbElements[ 0 ] ) ] )
308+
elementMap[ 3 ] = np.array( [ [ 1, element ] for element in sharedCells2D3DId[ :, 1 ] ] )
309+
elif meshToName == "emptyWell":
310+
elementMap[ 0 ] = np.full( ( nbElements[ 2 ], 2 ), -1, np.int64 )
311+
for id1DElem in range( nbElements[ 2 ] ):
312+
for sharedPoint1D3DId in sharedPoints1D3DId:
313+
if id1DElem == sharedPoint1D3DId[ 0 ]:
314+
elementMap[ 0 ][ id1DElem ] = [ 1, sharedPoint1D3DId[ 1 ] ]
315+
elif meshToName == "emptyFracture" or meshToName == "emptypolydata":
316+
elementMap[ 0 ] = np.array( [ [ 1, element ] for element in sharedCells2D3DId[ :, 1 ] ] )
317+
elif meshToName == "emptydataset":
318+
elementMap[ 0 ] = np.array( [ [ 1, element ] for element in range( nbElements[ 0 ] ) ] )
269319

270320
return elementMap
271321

0 commit comments

Comments
 (0)