Skip to content

Commit a0222b7

Browse files
alexbenedictoalexbenedicto
authored
fix: Invalid calculations in split mesh (#180)
* Correct invalid count of cells and false documentation * Add 2 test cases to check element counts created
1 parent 31b4c44 commit a0222b7

File tree

4 files changed

+281
-39
lines changed

4 files changed

+281
-39
lines changed

‎geos-processing/src/geos/processing/generic_processing_tools/SplitMesh.py‎

Lines changed: 27 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -27,14 +27,14 @@
2727
VTK_POLYHEDRON,
2828
VTK_POLYGON,
2929
)
30-
3130
from vtkmodules.util.numpy_support import ( numpy_to_vtk, vtk_to_numpy )
32-
3331
from geos.processing.pre_processing.CellTypeCounterEnhanced import CellTypeCounterEnhanced
3432
from geos.mesh.model.CellTypeCounts import CellTypeCounts
3533

3634
__doc__ = """
37-
SplitMesh module is a vtk filter that split cells of a mesh composed of Tetrahedra, pyramids, and hexahedra.
35+
SplitMesh module is a vtk filter that splits cells of a mesh composed of tetrahedra, pyramids, hexahedra, triangles, and quads.
36+
37+
Warning: Current implementation only supports meshes composed of either polygons or polyhedra, not both together.
3838
3939
Filter input and output types are vtkUnstructuredGrid.
4040
@@ -133,21 +133,24 @@ def RequestData(
133133
assert self.inData is not None, "Input mesh is undefined."
134134
assert output is not None, "Output mesh is undefined."
135135

136+
# Count the number of cells before splitting. Then we will be able to know how many new cells and points
137+
# to allocate because each cell type is splitted in a known number of new cells and points.
136138
nbCells: int = self.inData.GetNumberOfCells()
137139
counts: CellTypeCounts = self._getCellCounts()
138-
nbTet: int = counts.getTypeCount( VTK_TETRA )
139-
nbPyr: int = counts.getTypeCount( VTK_PYRAMID )
140-
nbHex: int = counts.getTypeCount( VTK_HEXAHEDRON )
141-
nbTriangles: int = counts.getTypeCount( VTK_TRIANGLE )
142-
nbQuad: int = counts.getTypeCount( VTK_QUAD )
140+
nbTet: int = counts.getTypeCount( VTK_TETRA ) # will divide into 8 tets
141+
nbPyr: int = counts.getTypeCount( VTK_PYRAMID ) # will divide into 6 pyramids and 4 tets so 10 new cells
142+
nbHex: int = counts.getTypeCount( VTK_HEXAHEDRON ) # will divide into 8 hexes
143+
nbTriangles: int = counts.getTypeCount( VTK_TRIANGLE ) # will divide into 4 triangles
144+
nbQuad: int = counts.getTypeCount( VTK_QUAD ) # will divide into 4 quads
143145
nbPolygon: int = counts.getTypeCount( VTK_POLYGON )
144146
nbPolyhedra: int = counts.getTypeCount( VTK_POLYHEDRON )
145147
assert counts.getTypeCount( VTK_WEDGE ) == 0, "Input mesh contains wedges that are not currently supported."
148+
# Current implementation only supports meshes composed of either polygons or polyhedra
146149
assert nbPolyhedra * nbPolygon == 0, ( "Input mesh is composed of both polygons and polyhedra,"
147150
" but it must contains only one of the two." )
148151
nbNewPoints: int = 0
149152
nbNewPoints = nbHex * 19 + nbTet * 6 + nbPyr * 9 if nbPolyhedra > 0 else nbTriangles * 3 + nbQuad * 5
150-
nbNewCells: int = nbHex * 8 + nbTet * 8 + nbPyr * 10 * nbTriangles * 4 + nbQuad * 4
153+
nbNewCells: int = nbHex * 8 + nbTet * 8 + nbPyr * 10 + nbTriangles * 4 + nbQuad * 4
151154

152155
self.points = vtkPoints()
153156
self.points.DeepCopy( self.inData.GetPoints() )
@@ -159,19 +162,20 @@ def RequestData(
159162
self.originalId.SetName( "OriginalID" )
160163
self.originalId.Allocate( nbNewCells )
161164
self.cellTypes = []
165+
# Define cell type to splitting method mapping
166+
splitMethods = {
167+
VTK_HEXAHEDRON: self._splitHexahedron,
168+
VTK_TETRA: self._splitTetrahedron,
169+
VTK_PYRAMID: self._splitPyramid,
170+
VTK_TRIANGLE: self._splitTriangle,
171+
VTK_QUAD: self._splitQuad,
172+
}
162173
for c in range( nbCells ):
163174
cell: vtkCell = self.inData.GetCell( c )
164175
cellType: int = cell.GetCellType()
165-
if cellType == VTK_HEXAHEDRON:
166-
self._splitHexahedron( cell, c )
167-
elif cellType == VTK_TETRA:
168-
self._splitTetrahedron( cell, c )
169-
elif cellType == VTK_PYRAMID:
170-
self._splitPyramid( cell, c )
171-
elif cellType == VTK_TRIANGLE:
172-
self._splitTriangle( cell, c )
173-
elif cellType == VTK_QUAD:
174-
self._splitQuad( cell, c )
176+
splitMethod = splitMethods.get( cellType )
177+
if splitMethod is not None:
178+
splitMethod( cell, c )
175179
else:
176180
raise TypeError( f"Cell type {vtkCellTypes.GetClassNameFromTypeId(cellType)} is not supported." )
177181
# add points and cells
@@ -261,7 +265,7 @@ def _splitPyramid( self: Self, cell: vtkCell, index: int ) -> None:
261265
r"""Split a pyramid.
262266
263267
Let's suppose an input pyramid composed of nodes (0, 1, 2, 3, 4),
264-
the cell is splitted in 8 pyramids using edge centers.
268+
the cell is split into 6 pyramids and 4 tetrahedra using edge centers.
265269
266270
4
267271
,/|\
@@ -310,7 +314,8 @@ def _splitPyramid( self: Self, cell: vtkCell, index: int ) -> None:
310314
self.cells.InsertNextCell( 4, [ pt12, pt7, pt6, pt13 ] )
311315
for _ in range( 10 ):
312316
self.originalId.InsertNextValue( index )
313-
self.cellTypes.extend( [ VTK_PYRAMID ] * 8 )
317+
self.cellTypes.extend( [ VTK_PYRAMID ] * 6 )
318+
self.cellTypes.extend( [ VTK_TETRA ] * 4 )
314319

315320
def _splitHexahedron( self: Self, cell: vtkCell, index: int ) -> None:
316321
r"""Split a hexahedron.
@@ -378,7 +383,7 @@ def _splitTriangle( self: Self, cell: vtkCell, index: int ) -> None:
378383
r"""Split a triangle.
379384
380385
Let's suppose an input triangle composed of nodes (0, 1, 2),
381-
the cell is splitted in 3 triangles using edge centers.
386+
the cell is split into 4 triangles using edge centers.
382387
383388
2
384389
|\

‎geos-processing/tests/data/hexs3_tets36_pyrs18.vtu‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
2+
<UnstructuredGrid>
3+
<Piece NumberOfPoints="38" NumberOfCells="57">
4+
<PointData>
5+
<DataArray type="Float32" Name="DistanceToCenter" format="binary" RangeMin="0" RangeMax="0.8660253882408142">
6+
mAAAAAAAAADXs10/17NdP9ezXT/Xs10/17NdP9ezXT/Xs10/17NdP9ezXT/Xs10/17NdP9ezXT/Xs10/17NdP9ezXT/Xs10/17NdP9ezXT/Xs10/17NdP9ezXT/Xs10/AAAAANezXT/Xs10/AAAAANezXT/Xs10/17NdP9ezXT/Xs10/17NdPwAAAADXs10/17NdPwAAAAAAAAAAAAAAAA==
7+
</DataArray>
8+
<DataArray type="Float32" Name="Polynomial" format="binary" RangeMin="1" RangeMax="4">
9+
mAAAAAAAAAAAAIA/AAAAQAAAAEAAAEBAAAAAQAAAQEAAAEBAAACAQAAAAEAAAABAAABAQAAAQEAAAEBAAACAQAAAAEAAAABAAABAQAAAQEAAAEBAAACAQAAAAEAAAEBAAACAPwAAAEAAAEBAAACAPwAAgD8AAABAAABAQAAAAEAAAEBAAACAQAAAgD8AAABAAABAQAAAgD8AAIA/AACAPw==
10+
</DataArray>
11+
</PointData>
12+
<CellData>
13+
</CellData>
14+
<Points>
15+
<DataArray type="Float32" Name="Points" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="4.358898943540674">
16+
yAEAAAAAAAAAAAAAAAAAAAAAAAAAAIA/AAAAAAAAAAAAAAAAAACAPwAAAAAAAIA/AACAPwAAAAAAAAAAAAAAAAAAgD8AAIA/AAAAAAAAgD8AAAAAAACAPwAAgD8AAIA/AACAPwAAgD8AAABAAACAPwAAAAAAAIA/AAAAQAAAAAAAAABAAAAAQAAAAAAAAABAAACAPwAAgD8AAIA/AAAAQAAAgD8AAABAAAAAQAAAgD8AAEBAAAAAQAAAAAAAAABAAABAQAAAAAAAAEBAAABAQAAAAAAAAEBAAAAAQAAAgD8AAABAAABAQAAAgD8AAEBAAABAQAAAgD8AAABAAAAAAAAAAAAAAABAAAAAAAAAgD8AAMA/AAAAPwAAAD8AAEBAAACAPwAAAAAAAEBAAACAPwAAgD8AACBAAADAPwAAAD8AAAAAAAAAQAAAAAAAAAAAAABAQAAAAAAAAIA/AABAQAAAAAAAAAAAAAAAQAAAgD8AAAAAAABAQAAAgD8AAIA/AABAQAAAgD8AAAA/AAAgQAAAAD8AAEBAAAAAAAAAAAAAAEBAAAAAAAAAgD8AACBAAAAAPwAAAD8AAAA/AADAPwAAAD8AAMA/AAAgQAAAAD8=
17+
<InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
18+
<Value index="0">
19+
0
20+
</Value>
21+
<Value index="1">
22+
4.3588989435
23+
</Value>
24+
</InformationKey>
25+
<InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
26+
<Value index="0">
27+
0
28+
</Value>
29+
<Value index="1">
30+
4.3588989435
31+
</Value>
32+
</InformationKey>
33+
</DataArray>
34+
</Points>
35+
<Cells>
36+
<DataArray type="Int64" Name="connectivity" format="binary" RangeMin="0" RangeMax="37">
37+
EAgAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAwAAAAAAAAACAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAHAAAAAAAAAAYAAAAAAAAAAwAAAAAAAAAIAAAAAAAAAAoAAAAAAAAACQAAAAAAAAAHAAAAAAAAAAsAAAAAAAAADQAAAAAAAAAMAAAAAAAAAAoAAAAAAAAADgAAAAAAAAAQAAAAAAAAAA8AAAAAAAAADQAAAAAAAAARAAAAAAAAABMAAAAAAAAAEgAAAAAAAAABAAAAAAAAABQAAAAAAAAACAAAAAAAAAAWAAAAAAAAAAEAAAAAAAAACAAAAAAAAAADAAAAAAAAABYAAAAAAAAACwAAAAAAAAAVAAAAAAAAAAUAAAAAAAAAFgAAAAAAAAALAAAAAAAAAAUAAAAAAAAABwAAAAAAAAAWAAAAAAAAABQAAAAAAAAAFQAAAAAAAAALAAAAAAAAABYAAAAAAAAAFAAAAAAAAAALAAAAAAAAAAgAAAAAAAAAFgAAAAAAAAABAAAAAAAAAAUAAAAAAAAAFQAAAAAAAAAWAAAAAAAAAAEAAAAAAAAAFQAAAAAAAAAUAAAAAAAAABYAAAAAAAAAAQAAAAAAAAADAAAAAAAAAAcAAAAAAAAAFgAAAAAAAAABAAAAAAAAAAcAAAAAAAAABQAAAAAAAAAWAAAAAAAAAAsAAAAAAAAABwAAAAAAAAADAAAAAAAAABYAAAAAAAAACwAAAAAAAAADAAAAAAAAAAgAAAAAAAAAFgAAAAAAAAAIAAAAAAAAABcAAAAAAAAADgAAAAAAAAAZAAAAAAAAAAgAAAAAAAAADgAAAAAAAAAKAAAAAAAAABkAAAAAAAAAEQAAAAAAAAAYAAAAAAAAAAsAAAAAAAAAGQAAAAAAAAARAAAAAAAAAAsAAAAAAAAADQAAAAAAAAAZAAAAAAAAABcAAAAAAAAAGAAAAAAAAAARAAAAAAAAABkAAAAAAAAAFwAAAAAAAAARAAAAAAAAAA4AAAAAAAAAGQAAAAAAAAAIAAAAAAAAAAsAAAAAAAAAGAAAAAAAAAAZAAAAAAAAAAgAAAAAAAAAGAAAAAAAAAAXAAAAAAAAABkAAAAAAAAACAAAAAAAAAAKAAAAAAAAAA0AAAAAAAAAGQAAAAAAAAAIAAAAAAAAAA0AAAAAAAAACwAAAAAAAAAZAAAAAAAAABEAAAAAAAAADQAAAAAAAAAKAAAAAAAAABkAAAAAAAAAEQAAAAAAAAAKAAAAAAAAAA4AAAAAAAAAGQAAAAAAAAAaAAAAAAAAAAkAAAAAAAAAHAAAAAAAAAAgAAAAAAAAABoAAAAAAAAAHAAAAAAAAAAbAAAAAAAAACAAAAAAAAAAHwAAAAAAAAAMAAAAAAAAAB0AAAAAAAAAIAAAAAAAAAAfAAAAAAAAAB0AAAAAAAAAHgAAAAAAAAAgAAAAAAAAAAkAAAAAAAAADAAAAAAAAAAfAAAAAAAAACAAAAAAAAAACQAAAAAAAAAfAAAAAAAAABwAAAAAAAAAIAAAAAAAAAAaAAAAAAAAAB0AAAAAAAAADAAAAAAAAAAgAAAAAAAAABoAAAAAAAAADAAAAAAAAAAJAAAAAAAAACAAAAAAAAAAGgAAAAAAAAAbAAAAAAAAAB4AAAAAAAAAIAAAAAAAAAAaAAAAAAAAAB4AAAAAAAAAHQAAAAAAAAAgAAAAAAAAAB8AAAAAAAAAHgAAAAAAAAAbAAAAAAAAACAAAAAAAAAAHwAAAAAAAAAbAAAAAAAAABwAAAAAAAAAIAAAAAAAAAAUAAAAAAAAACEAAAAAAAAAFwAAAAAAAAAIAAAAAAAAACMAAAAAAAAAGAAAAAAAAAAiAAAAAAAAABUAAAAAAAAACwAAAAAAAAAjAAAAAAAAACEAAAAAAAAAIgAAAAAAAAAYAAAAAAAAABcAAAAAAAAAIwAAAAAAAAAUAAAAAAAAABUAAAAAAAAAIgAAAAAAAAAhAAAAAAAAACMAAAAAAAAAFAAAAAAAAAAIAAAAAAAAAAsAAAAAAAAAFQAAAAAAAAAjAAAAAAAAABgAAAAAAAAACwAAAAAAAAAIAAAAAAAAABcAAAAAAAAAIwAAAAAAAAACAAAAAAAAAAMAAAAAAAAACQAAAAAAAAAaAAAAAAAAACQAAAAAAAAADAAAAAAAAAAHAAAAAAAAAAYAAAAAAAAAHQAAAAAAAAAkAAAAAAAAAAMAAAAAAAAABwAAAAAAAAAMAAAAAAAAAAkAAAAAAAAAJAAAAAAAAAACAAAAAAAAAAYAAAAAAAAABwAAAAAAAAADAAAAAAAAACQAAAAAAAAAAgAAAAAAAAAaAAAAAAAAAB0AAAAAAAAABgAAAAAAAAAkAAAAAAAAAAwAAAAAAAAAHQAAAAAAAAAaAAAAAAAAAAkAAAAAAAAAJAAAAAAAAAAJAAAAAAAAAAoAAAAAAAAADwAAAAAAAAAcAAAAAAAAACUAAAAAAAAAEgAAAAAAAAANAAAAAAAAAAwAAAAAAAAAHwAAAAAAAAAlAAAAAAAAAAoAAAAAAAAADQAAAAAAAAASAAAAAAAAAA8AAAAAAAAAJQAAAAAAAAAJAAAAAAAAAAwAAAAAAAAADQAAAAAAAAAKAAAAAAAAACUAAAAAAAAACQAAAAAAAAAcAAAAAAAAAB8AAAAAAAAADAAAAAAAAAAlAAAAAAAAABIAAAAAAAAAHwAAAAAAAAAcAAAAAAAAAA8AAAAAAAAAJQAAAAAAAAA=
38+
</DataArray>
39+
<DataArray type="Int64" Name="offsets" format="binary" RangeMin="8" RangeMax="258">
40+
yAEAAAAAAAAIAAAAAAAAABAAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAoAAAAAAAAACwAAAAAAAAAMAAAAAAAAAA0AAAAAAAAADgAAAAAAAAAPAAAAAAAAABAAAAAAAAAAEQAAAAAAAAASAAAAAAAAABMAAAAAAAAAFAAAAAAAAAAVAAAAAAAAABYAAAAAAAAAFwAAAAAAAAAYAAAAAAAAABkAAAAAAAAAGgAAAAAAAAAbAAAAAAAAABwAAAAAAAAAHQAAAAAAAAAeAAAAAAAAAB8AAAAAAAAAIAAAAAAAAAAhAAAAAAAAACIAAAAAAAAAIwAAAAAAAAAkAAAAAAAAACUAAAAAAAAAJgAAAAAAAAAnAAAAAAAAACgAAAAAAAAAKQAAAAAAAAAqAAAAAAAAACtAAAAAAAAALIAAAAAAAAAtwAAAAAAAAC8AAAAAAAAAMEAAAAAAAAAxgAAAAAAAADLAAAAAAAAANAAAAAAAAAA1QAAAAAAAADaAAAAAAAAAN8AAAAAAAAA5AAAAAAAAADpAAAAAAAAAO4AAAAAAAAA8wAAAAAAAAD4AAAAAAAAAP0AAAAAAAAAAgEAAAAAAAA=
41+
</DataArray>
42+
<DataArray type="UInt8" Name="types" format="binary" RangeMin="10" RangeMax="14">
43+
OQAAAAAAAAAMDAwKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoODg4ODg4ODg4ODg4ODg4ODg4=
44+
</DataArray>
45+
</Cells>
46+
</Piece>
47+
</UnstructuredGrid>
48+
</VTKFile>

‎geos-processing/tests/data/quads2_tris4.vtu‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
2+
<UnstructuredGrid>
3+
<Piece NumberOfPoints="9" NumberOfCells="6">
4+
<PointData>
5+
<DataArray type="Float32" Name="DistanceToCenter" format="binary" RangeMin="0.7071067690849304" RangeMax="0.7071067690849304">
6+
JAAAAAAAAADzBDU/8wQ1P/MENT/zBDU/8wQ1P/MENT/zBDU/8wQ1P/MENT8=
7+
</DataArray>
8+
<DataArray type="Float32" Name="Polynomial" format="binary" RangeMin="1" RangeMax="3">
9+
JAAAAAAAAAAAAIA/AAAAQAAAAEAAAEBAAAAAQAAAAEAAAEBAAAAAQAAAAEA=
10+
</DataArray>
11+
</PointData>
12+
<CellData>
13+
</CellData>
14+
<Points>
15+
<DataArray type="Float32" Name="Points" NumberOfComponents="3" format="binary" RangeMin="0" RangeMax="2.8284271247461903">
16+
bAAAAAAAAAAAAAAAAAAAAAAAAAAAAIA/AAAAAAAAAAAAAAAAAACAPwAAAAAAAIA/AACAPwAAAAAAAABAAACAPwAAAAAAAIA/AAAAQAAAAAAAAABAAAAAQAAAAAAAAABAAAAAAAAAAAAAAAAAAAAAQAAAAAA=
17+
<InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
18+
<Value index="0">
19+
0
20+
</Value>
21+
<Value index="1">
22+
2.8284271247
23+
</Value>
24+
</InformationKey>
25+
<InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
26+
<Value index="0">
27+
0
28+
</Value>
29+
<Value index="1">
30+
2.8284271247
31+
</Value>
32+
</InformationKey>
33+
</DataArray>
34+
</Points>
35+
<Cells>
36+
<DataArray type="Int64" Name="connectivity" format="binary" RangeMin="0" RangeMax="8">
37+
oAAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAwAAAAAAAAACAAAAAAAAAAMAAAAAAAAABAAAAAAAAAAGAAAAAAAAAAUAAAAAAAAAAQAAAAAAAAAHAAAAAAAAAAMAAAAAAAAABwAAAAAAAAAEAAAAAAAAAAMAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAgAAAAAAAAAAwAAAAAAAAAFAAAAAAAAAAgAAAAAAAAA
38+
</DataArray>
39+
<DataArray type="Int64" Name="offsets" format="binary" RangeMin="4" RangeMax="20">
40+
MAAAAAAAAAAEAAAAAAAAAAgAAAAAAAAACwAAAAAAAAAOAAAAAAAAABEAAAAAAAAAFAAAAAAAAAA=
41+
</DataArray>
42+
<DataArray type="UInt8" Name="types" format="binary" RangeMin="5" RangeMax="9">
43+
BgAAAAAAAAAJCQUFBQU=
44+
</DataArray>
45+
</Cells>
46+
</Piece>
47+
</UnstructuredGrid>
48+
</VTKFile>

0 commit comments

Comments
 (0)