diff options
| -rwxr-xr-x | pyfrm2xdmf | 223 |
1 files changed, 151 insertions, 72 deletions
@@ -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() |