From b13b451f1f8359087f2dce7e7ec444c617b2d1bc Mon Sep 17 00:00:00 2001 From: Paul Garlick Date: Thu, 24 Mar 2016 19:37:07 +0000 Subject: pyfr2xdmf: allow first and second order hexahedral cell types --- pyfrm2xdmf | 223 +++++++++++++++++++++++++++++++++++++++++-------------------- 1 file 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(' \n') return -def writeTopology(xdmfFile, nCells, connFile): +def writeTopology(xdmfFile, tType, nCells, connFile): "write Topology element" - xdmfFile.write(' \n'.format(nCells)) + xdmfFile.write(' \n'.format(tType, nCells)) xdmfFile.write(' \n'.format(connFile)) xdmfFile.write(' \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(' \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(' \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(' \n') xdmfFile.write(' \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(' \n') @@ -78,7 +80,7 @@ def writeAttribute(xdmfFile, tag): xdmfFile.write(' \n') return -def writeConnectivities(connFile, nCells, nVerts, orderDict): +def writeConnectivities(connFile, nCells, nVerts): "write connectivities to xml file" cf = open(connFile, 'w') cf.write('\n') @@ -102,76 +104,153 @@ def writeFooter(xdmfFile): xdmfFile.write('\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(' \n'.format(part)) - else: - g.write(' \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(' \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(' \n') - # connectivities file - writeConnectivities(xfname, partitions[cellType][part], - nverts[cellTypes[cellType]], order[cellType]) - - writeAttribute(g, part) - g.write(' \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(' \n'.format(part)) + else: + g.write(' \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(' \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(' \n') + # connectivities file + writeConnectivities(xfname, partitions[cellType][part], + len(nodeIDs[cellType])) + + writeAttribute(g, part) + g.write(' \n') + + writeFooter(g) + g.close() -- cgit