1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# run within virtual environment (uses Python 3 syntax)
import argparse
import os
import re
from array import array
from collections import Counter
from subprocess import check_output
def meshFile(param):
global base
base, ext = os.path.splitext(param)
if ext.lower() != '.pyfrm':
raise argparse.ArgumentTypeError('Mesh file must have a .pyfrm extension')
return param
def writeHeader(xdmfFile):
"write XDMF header"
xdmfFile.write('<?xml version="1.0" ?>\n')
xdmfFile.write('<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>\n')
xdmfFile.write('<Xdmf xmlns:xi="http://www.w3.org/2003/XInclude" Version="2.2">\n')
xdmfFile.write(' <Domain>\n')
return
def writeTopology(xdmfFile, nCells, connFile):
"write Topology element"
xdmfFile.write(' <Topology TopologyType="Quadrilateral" NumberOfElements="{}">\n'.format(nCells))
xdmfFile.write(' <xi:include href="{}"/>\n'.format(connFile))
xdmfFile.write(' </Topology>\n')
return
def writeGeometry(xdmfFile, nDims, nCells, nVerts, pyfrm, dataset):
"write Geometry element"
if nDims == 2:
xdmfFile.write(' <Geometry GeometryType="X_Y">\n') # co-ordinates in separate arrays
else:
xdmfFile.write(' <Geometry GeometryType="X_Y_Z">\n') # co-ordinates in separate arrays
writeHyperSlab(xdmfFile, nDims, nCells, nVerts, pyfrm, dataset)
xdmfFile.write(' </Geometry>\n')
return
def writeHyperSlab(xdmfFile, nDims, nCells, nVerts, pyfrm, dataset):
"write HyperSlab element"
for coord in range(nDims):
xdmfFile.write(' <DataItem ItemType="HyperSlab"\n')
xdmfFile.write(' Dimensions="{} 1 1"\n'.format(nCells*nVerts))
xdmfFile.write(' Type="HyperSlab">\n')
xdmfFile.write(' <DataItem\n') # start, stride and count of hyperslab region
xdmfFile.write(' Dimensions="3 3"\n')
xdmfFile.write(' Format="XML">\n')
xdmfFile.write(' 0 0 {}\n'.format(coord)) # select co-ordinate
xdmfFile.write(' 1 1 1\n') # select every vertex in every cell
xdmfFile.write(' {:<3} {} 1\n'.format(nVerts, nCells)) # loop over cells (first) and vertices (second)
xdmfFile.write(' </DataItem>\n')
xdmfFile.write(' <DataItem\n')
xdmfFile.write(' Name="Points" \n')
xdmfFile.write(' Dimensions="{} {} {}"\n'.format(nVerts, nCells, nDims))
xdmfFile.write(' Format="HDF">\n')
xdmfFile.write(' {}:/{}\n'.format(pyfrm, dataset))
xdmfFile.write(' </DataItem>\n')
xdmfFile.write(' </DataItem>\n')
return
def writeAttribute(xdmfFile, tag):
"write Attribute element"
xdmfFile.write(' <Attribute Name="Partition" Center="Grid">\n')
xdmfFile.write(' <DataItem\n')
xdmfFile.write(' Dimensions="1"\n')
xdmfFile.write(' Format="XML">\n')
xdmfFile.write(' {}\n'.format(tag)) # tag with partition number
xdmfFile.write(' </DataItem>\n')
xdmfFile.write(' </Attribute>\n')
return
def writeConnectivities(connFile, nCells, nVerts, orderDict):
"write connectivities to xml file"
cf = open(connFile, 'w')
cf.write('<DataItem DataType="Int"\n')
cf.write(' Dimensions="{} {}"\n'.format(nCells, nVerts))
cf.write(' Format="XML">\n')
for i in range (0, nCells):
cf.write(' ')
for j in range (1, nVerts+1):
cf.write(' ' + repr(orderDict[j]*nCells+i).ljust(1))
cf.write('\n')
cf.write('</DataItem>\n')
cf.close()
print('connectivities written to ' + connFile)
return
def writeFooter(xdmfFile):
"write XDMF footer"
xdmfFile.write(' </Domain>\n')
xdmfFile.write('</Xdmf>\n')
return
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
# cell types
cellTypes = ['quad', 'tri']
numCellTypes = len(cellTypes)
nverts = {cellTypes[0]: 4, cellTypes[1]: 3}
ndims = {cellTypes[0]: 2, cellTypes[1]: 2}
# 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
# 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()
|
