aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorPaul Garlick <pgarlick@tourbillion-technology.com>2016-03-24 19:37:07 +0000
committerPaul Garlick <pgarlick@tourbillion-technology.com>2016-03-24 19:37:07 +0000
commitb13b451f1f8359087f2dce7e7ec444c617b2d1bc (patch)
tree64f27961106fba5da615a2eecfce610d25c65579
parentc98547931d611bc5d5dd754ef29594b863f28eb9 (diff)
downloadpyfrUtils-b13b451f1f8359087f2dce7e7ec444c617b2d1bc.tar.gz
pyfr2xdmf: allow first and second order hexahedral cell types
-rwxr-xr-xpyfrm2xdmf223
1 files changed, 151 insertions, 72 deletions
diff --git a/pyfrm2xdmf b/pyfrm2xdmf
index 7b88716..172b76c 100755
--- a/pyfrm2xdmf
+++ b/pyfrm2xdmf
@@ -24,15 +24,16 @@ def writeHeader(xdmfFile):
xdmfFile.write(' <Domain>\n')
return
-def writeTopology(xdmfFile, nCells, connFile):
+def writeTopology(xdmfFile, tType, nCells, connFile):
"write Topology element"
- xdmfFile.write(' <Topology TopologyType="Quadrilateral" NumberOfElements="{}">\n'.format(nCells))
+ xdmfFile.write(' <Topology TopologyType="{}" NumberOfElements="{}">\n'.format(tType, nCells))
xdmfFile.write(' <xi:include href="{}"/>\n'.format(connFile))
xdmfFile.write(' </Topology>\n')
return
-def writeGeometry(xdmfFile, nDims, nCells, nVerts, pyfrm, dataset):
+def writeGeometry(xdmfFile, nDims, nCells, nodeArray, pyfrm, dataset):
"write Geometry element"
+ nVerts = len(nodeArray)
if nDims == 2:
xdmfFile.write(' <Geometry GeometryType="X_Y">\n') # co-ordinates in separate arrays
else:
@@ -40,14 +41,15 @@ def writeGeometry(xdmfFile, nDims, nCells, nVerts, pyfrm, dataset):
for coord in range(nDims):
xdmfFile.write(' <DataItem ItemType="Function" Dimensions="{}"\n'.format(nVerts*nCells)) # 1D-array
xdmfFile.write(' Function="JOIN({})">\n'.format(' ; '.join("$" + str(k) for k in range(nVerts))))
- writeHyperSlab(xdmfFile, coord, nDims, nCells, nVerts, pyfrm, dataset)
+ writeHyperSlab(xdmfFile, coord, nDims, nCells, nodeArray, pyfrm, dataset)
xdmfFile.write(' </DataItem>\n')
xdmfFile.write(' </Geometry>\n')
return
-def writeHyperSlab(xdmfFile, coord, nDims, nCells, nVerts, pyfrm, dataset):
+def writeHyperSlab(xdmfFile, coord, nDims, nCells, nodeArray, pyfrm, dataset):
"write HyperSlab element"
- for vert in range(nVerts):
+ nVerts = len(nodeArray)
+ for vert in nodeArray:
xdmfFile.write(' <DataItem ItemType="HyperSlab"\n')
xdmfFile.write(' Dimensions="{} 1 1"\n'.format(nCells))
xdmfFile.write(' Type="HyperSlab">\n')
@@ -78,7 +80,7 @@ def writeAttribute(xdmfFile, tag):
xdmfFile.write(' </Attribute>\n')
return
-def writeConnectivities(connFile, nCells, nVerts, orderDict):
+def writeConnectivities(connFile, nCells, nVerts):
"write connectivities to xml file"
cf = open(connFile, 'w')
cf.write('<DataItem DataType="Int"\n')
@@ -87,8 +89,8 @@ def writeConnectivities(connFile, nCells, nVerts, orderDict):
for i in range (0, nCells):
cf.write(' ')
- for j in range (1, nVerts+1):
- cf.write(' ' + repr(orderDict[j]*nCells+i).ljust(1))
+ for j in range(nVerts):
+ cf.write(' ' + repr(j*nCells+i).ljust(1))
cf.write('\n')
cf.write('</DataItem>\n')
@@ -102,76 +104,153 @@ def writeFooter(xdmfFile):
xdmfFile.write('</Xdmf>\n')
return
+def readH5lsOutput(os_output):
+ "read output from 'h5ls' command"
+ tet4, tet10, tet20, tet35, tet56, tet84 = ({} for i in range(6)) # tet cells
+ pri6, pri18, pri40, pri75, pri126, pri196 = ({} for i in range(6)) # pri cells
+ pyr5, pyr14, pyr30, pyr55, pyr91, pyr140 = ({} for i in range(6)) # pyr cells
+ hex8, hex27, hex64, hex125, hex216 = ({} for i in range(5)) # hex cells
+ tri3, tri6, tri10, tri15, tri21 = ({} for i in range(5)) # tri cells
+ quad4, quad9, quad16, quad25, quad36 = ({} for i in range(5)) # quad cells
+
+ for line in os_output.splitlines():
+ spt = re.search('spt', line.decode()) # restrict to 'spt' arrays
+ if spt:
+ chunk = line.decode().split()
+ partno = int(re.search('\d+', chunk[0]).group())
+ nnodes = int(re.search('(\d+),', line.decode()).group(1))
+ ncells = int(re.search(' (\d+),', line.decode()).group(1))
+ if re.search('quad', line.decode()): # check whether cell is quadrangular
+ if nnodes == 4:
+ quad4[partno] = ncells
+ elif nnodes == 9:
+ quad9[partno] = ncells
+ elif nnodes == 16:
+ quad16[partno] = ncells
+ elif nnodes == 25:
+ quad25[partno] = ncells
+ elif nnodes == 36:
+ quad36[partno] = ncells
+ else:
+ print("unknown cell order")
+ elif re.search('tri', line.decode()): # check whether cell is triangular
+ if nnodes == 3:
+ tri3[partno] = ncells
+ elif nnodes == 6:
+ tri6[partno] = ncells
+ elif nnodes == 10:
+ tri10[partno] = ncells
+ elif nnodes == 15:
+ tri15[partno] = ncells
+ elif nnodes == 21:
+ tri21[partno] = ncells
+ else:
+ print("unknown cell order")
+ elif re.search('hex', line.decode()): # check whether cell is hexahedral
+ if nnodes == 8:
+ hex8[partno] = ncells
+ elif nnodes == 27:
+ hex27[partno] = ncells
+ elif nnodes == 64:
+ hex64[partno] = ncells
+ elif nnodes == 125:
+ hex125[partno] = ncells
+ elif nnodes == 216:
+ hex216[partno] = ncells
+ else:
+ print("unknown cell order")
+ else:
+ print("unknown cell type")
+ break
+
+ return (tet4, tet10, tet20, tet35, tet56, tet84,
+ pri6, pri18, pri40, pri75, pri126, pri196,
+ pyr5, pyr14, pyr30, pyr55, pyr91, pyr140,
+ hex8, hex27, hex64, hex125, hex216,
+ tri3, tri6, tri10, tri15, tri21,
+ quad4, quad9, quad16, quad25, quad36)
+
+# read command line arguments
parser = argparse.ArgumentParser(description="extract connectivities from mesh file")
parser.add_argument("mesh", help="mesh file (.pyfrm)", type=meshFile)
args = parser.parse_args()
# use 'h5ls' command to provide array dimensions
-h5ls_output = check_output(["h5ls", args.mesh])
-
-nquads = {}
-ntris = {}
-for line in h5ls_output.splitlines():
- spt = re.search('spt', line.decode()) # restrict to 'spt' arrays
- if spt:
- chunk = line.decode().split()
- npart = int(re.search('\d+', chunk[0]).group())
- ncells = int(re.search(' (\d+),', line.decode()).group(1))
- if re.search('quad', line.decode()): # check whether cell is quadrilateral
- nquads[npart] = ncells
- elif re.search('tri', line.decode()): # check whether cell is triangular
- ntris[npart] = ncells
- else:
- print("unknown cell type")
- break
+h5lsOutput = check_output(["h5ls", args.mesh])
+
+(ntet4, ntet10, ntet20, ntet35, ntet56, ntet84,
+npri6, npri18, npri40, npri75, npri126, npri196,
+npyr5, npyr14, npyr30, npyr55, npyr91, npyr140,
+nhex8, nhex27, nhex64, nhex125, nhex216,
+ntri3, ntri6, ntri10, ntri15, ntri21,
+nquad4, nquad9, nquad16, nquad25, nquad36) = readH5lsOutput(h5lsOutput)
# cell types
-cellTypes = ['quad', 'tri']
-numCellTypes = len(cellTypes)
-nverts = {cellTypes[0]: 4, cellTypes[1]: 3}
-ndims = {cellTypes[0]: 2, cellTypes[1]: 2}
+firstOrderCellType = ['tet', 'tet', 'tet', 'tet', 'tet', 'tet',
+ 'pri', 'pri', 'pri', 'pri', 'pri', 'pri',
+ 'pyr', 'pyr', 'pyr', 'pyr', 'pyr', 'pyr',
+ 'hex', 'hex', 'hex', 'hex', 'hex',
+ 'tri', 'tri', 'tri', 'tri', 'tri',
+ 'quad', 'quad', 'quad', 'quad', 'quad']
+ndims = {'quad': 2, 'tri': 2, 'hex': 3, 'pyr': 3, 'pri': 3, 'tet':3}
+xdmfTopologyType = {'quad': 'Quadrilateral', 'tri': 'Triangle', 'hex': 'Hexahedron',
+ 'pyr' : 'Pyramid', 'pri': 'Wedge', 'tet': 'Tetrahedron'}
-# XDMF:PyFR vertex numbering
-order = []
-order.append({1:0, 2:1, 3:3, 4:2}) # quad order
-order.append({1:0, 2:1, 3:2}) # tri order
+# node identification (see pyfr/readers/nodemaps.py): reduce high-order cells to first order
+nodeIDs = {}
+nodeIDs[29] = [0, 2, 8, 6] # quad9
+nodeIDs[28] = [0, 1, 3, 2] # quad4
+nodeIDs[23] = [0, 1, 2] # tri3
+nodeIDs[19] = [0, 2, 8, 6, 18, 20, 26, 24] # hex27
+nodeIDs[18] = [0, 1, 3, 2, 4, 5, 7, 6] # hex8
# sort datasets
-quadKeys = list(nquads.keys()) # keys are partition numbers
-triKeys = list(ntris.keys()) # keys are partition numbers
-allKeys = quadKeys + triKeys # concatenate keys
-numTypes = Counter(allKeys) # number of types present in each partition
-partKeys = list(numTypes.keys()) # partition keys
-partitions = [nquads, ntris] # list of dictionaries
-
-# write files
-g = open(os.path.join(base + '.xdmf'), 'w')
-writeHeader(g)
-
-for part in partKeys:
- if numTypes[part] > 1: # check whether partition contain multiple cell types
- g.write(' <Grid Name="Partition{}" GridType="Collection">\n'.format(part))
- else:
- g.write(' <Grid Name="Partition{}" GridType="Uniform">\n'.format(part))
-
- for cellType in range(numCellTypes):
- if part in partitions[cellType]: # check whether these cells exist in this partition
- xfname = os.path.join('con_' + cellTypes[cellType] + '_p' + str(part) + '.xml')
- dname = os.path.join('spt_' + cellTypes[cellType] + '_p' + str(part))
- if numTypes[part] > 1:
- g.write(' <Grid Name="cellType{}" GridType="Uniform">\n'.format(cellType))
- writeTopology(g, partitions[cellType][part], xfname)
- writeGeometry(g, ndims[cellTypes[cellType]],
- partitions[cellType][part],
- nverts[cellTypes[cellType]], args.mesh, dname)
- if numTypes[part] > 1:
- g.write(' </Grid>\n')
- # connectivities file
- writeConnectivities(xfname, partitions[cellType][part],
- nverts[cellTypes[cellType]], order[cellType])
-
- writeAttribute(g, part)
- g.write(' </Grid>\n')
-
-writeFooter(g)
-g.close()
+hex8Keys = list(nhex8.keys()) # keys are partition numbers
+hex27Keys = list(nhex27.keys())
+quad4Keys = list(nquad4.keys())
+tri3Keys = list(ntri3.keys())
+allKeys = quad4Keys + tri3Keys + hex8Keys + hex27Keys # concatenate keys
+numTypes = Counter(allKeys) # number of types present in each partition
+partKeys = list(numTypes.keys()) # partition keys
+partitions = [ntet4, ntet10, ntet20, ntet35, ntet56, ntet84,
+ npri6, npri18, npri40, npri75, npri126, npri196,
+ npyr5, npyr14, npyr30, npyr55, npyr91, npyr140,
+ nhex8, nhex27, nhex64, nhex125, nhex216,
+ ntri3, ntri6, ntri10, ntri15, ntri21,
+ nquad4, nquad9, nquad16, nquad25, nquad36] # list of cell dictionaries
+numCellTypes = len(partitions)
+
+if partKeys:
+ # write files
+ g = open(os.path.join(base + '.xdmf'), 'w')
+ writeHeader(g)
+
+ for part in partKeys:
+ if numTypes[part] > 1: # check whether partition contain multiple cell types
+ g.write(' <Grid Name="Partition{}" GridType="Collection">\n'.format(part))
+ else:
+ g.write(' <Grid Name="Partition{}" GridType="Uniform">\n'.format(part))
+
+ for cellType in range(numCellTypes):
+ if part in partitions[cellType]: # check whether these cells exist in this partition
+ xfname = os.path.join('con_' + firstOrderCellType[cellType] + '_p' + str(part) + '.xml')
+ dname = os.path.join('spt_' + firstOrderCellType[cellType] + '_p' + str(part))
+ if numTypes[part] > 1:
+ g.write(' <Grid Name="cellType{}" GridType="Uniform">\n'.format(cellType))
+ writeTopology(g, xdmfTopologyType[firstOrderCellType[cellType]],
+ partitions[cellType][part], xfname)
+ writeGeometry(g, ndims[firstOrderCellType[cellType]],
+ partitions[cellType][part],
+ nodeIDs[cellType],
+ args.mesh, dname)
+ if numTypes[part] > 1:
+ g.write(' </Grid>\n')
+ # connectivities file
+ writeConnectivities(xfname, partitions[cellType][part],
+ len(nodeIDs[cellType]))
+
+ writeAttribute(g, part)
+ g.write(' </Grid>\n')
+
+ writeFooter(g)
+ g.close()