pyfr2xdmf: allow first and second order hexahedral cell types

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()