1#!/usr/local/bin/python3.8
2import h5py
3import numpy as np
4import os, sys
5
6class Xdmf:
7  def __init__(self, filename):
8    self.filename = filename
9    self.cellMap  = {1 : {1 : 'Polyvertex', 2 : 'Polyline'}, 2 : {3 : 'Triangle', 4 : 'Quadrilateral'}, 3 : {4 : 'Tetrahedron', 6: 'Wedge', 8 : 'Hexahedron'}}
10
11    # py2/py3 compatibility, see https://github.com/h5py/h5py/issues/379
12    if sys.version_info[0] < 3:
13      self.typeMap = {'scalar' : 'Scalar', 'vector' : 'Vector', 'tensor' : 'Tensor6', 'matrix' : 'Matrix'}
14      self.typeExt = {2 : {'vector' : ['x', 'y'], 'tensor' : ['xx', 'yy', 'xy']}, 3 : {'vector' : ['x', 'y', 'z'], 'tensor' : ['xx', 'yy', 'zz', 'xy', 'yz', 'xz']}}
15    else:
16      self.typeMap = {b'scalar' : 'Scalar', b'vector' : 'Vector', b'tensor' : 'Tensor6', b'matrix' : 'Matrix'}
17      self.typeExt = {2 : {b'vector' : ['x', 'y'], b'tensor' : ['xx', 'yy', 'xy']}, 3 : {b'vector' : ['x', 'y', 'z'], b'tensor' : ['xx', 'yy', 'zz', 'xy', 'yz', 'xz']}}
18    return
19
20  def writeHeader(self, fp, hdfFilename):
21    fp.write('''\
22<?xml version="1.0" ?>
23<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" [
24<!ENTITY HeavyData "%s">
25]>
26''' % os.path.basename(hdfFilename))
27    fp.write('\n<Xdmf>\n  <Domain Name="domain">\n')
28    return
29
30  def writeCells(self, fp, topologyPath, numCells, numCorners, cellsName = "cells"):
31    fp.write('''\
32    <DataItem Name="%s"
33              ItemType="Uniform"
34              Format="HDF"
35              NumberType="Float" Precision="8"
36              Dimensions="%d %d">
37      &HeavyData;:/%s/cells
38    </DataItem>
39''' % (cellsName, numCells, numCorners, topologyPath))
40    return
41
42  def writeVertices(self, fp, geometryPath, numVertices, spaceDim):
43    fp.write('''\
44    <DataItem Name="vertices"
45              Format="HDF"
46              Dimensions="%d %d">
47      &HeavyData;:/%s/vertices
48    </DataItem>
49    <!-- ============================================================ -->
50''' % (numVertices, spaceDim, geometryPath))
51    return
52
53  def writeLocations(self, fp, numParticles, spaceDim):
54    fp.write('''\
55    <DataItem Name="particle_coordinates"
56              Format="HDF"
57              Dimensions="%d %d">
58      &HeavyData;:/particles/coordinates
59    </DataItem>
60''' % (numParticles, spaceDim))
61    return
62
63  def writeTimeGridHeader(self, fp, time):
64    fp.write('''\
65    <Grid Name="TimeSeries" GridType="Collection" CollectionType="Temporal">
66      <Time TimeType="List">
67        <DataItem Format="XML" NumberType="Float" Dimensions="%d">
68          ''' % (len(time)))
69    fp.write(' '.join([str(float(t)) for t in time]))
70    fp.write('''
71        </DataItem>
72      </Time>
73''')
74    return
75
76  #http://www.xdmf.org/index.php/XDMF_Model_and_Format#Topology
77  def writeHybridSpaceGridHeader(self, fp):
78    fp.write('      <Grid Name="domain" GridType="Collection">\n')
79    return
80
81  def writeSpaceGridHeader(self, fp, numCells, numCorners, cellDim, spaceDim, cellsName = "cells"):
82    fp.write('''\
83      <Grid Name="domain" GridType="Uniform">
84        <Topology
85           TopologyType="%s"
86           NumberOfElements="%d">
87          <DataItem Reference="XML">
88            /Xdmf/Domain/DataItem[@Name="%s"]
89          </DataItem>
90        </Topology>
91        <Geometry GeometryType="%s">
92          <DataItem Reference="XML">
93            /Xdmf/Domain/DataItem[@Name="vertices"]
94          </DataItem>
95        </Geometry>
96''' % (self.cellMap[cellDim][numCorners], numCells, cellsName, "XYZ" if spaceDim > 2 else "XY"))
97    return
98
99  def writeFieldSingle(self, fp, numSteps, timestep, spaceDim, name, f, domain):
100    if len(f[1].shape) > 2:
101      dof = f[1].shape[1]
102      bs  = f[1].shape[2]
103    elif len(f[1].shape) > 1:
104      if numSteps > 1:
105        dof = f[1].shape[1]
106        bs  = 1
107      else:
108        dof = f[1].shape[0]
109        bs  = f[1].shape[1]
110    else:
111      dof = f[1].shape[0]
112      bs  = 1
113    fp.write('''\
114        <Attribute
115           Name="%s"
116           Type="%s"
117           Center="%s">
118          <DataItem ItemType="HyperSlab"
119        	    Dimensions="1 %d %d"
120        	    Type="HyperSlab">
121            <DataItem
122               Dimensions="3 3"
123               Format="XML">
124              %d 0 0
125              1 1 1
126              1 %d %d
127            </DataItem>
128            <DataItem
129               DataType="Float" Precision="8"
130               Dimensions="%d %d %d"
131               Format="HDF">
132              &HeavyData;:%s
133            </DataItem>
134          </DataItem>
135        </Attribute>
136''' % (f[0], self.typeMap[f[1].attrs['vector_field_type']], domain, dof, bs, timestep, dof, bs, numSteps, dof, bs, name))
137    return
138
139  def writeFieldComponents(self, fp, numSteps, timestep, spaceDim, name, f, domain):
140    vtype = f[1].attrs['vector_field_type']
141    if len(f[1].shape) > 2:
142      dof    = f[1].shape[1]
143      bs     = f[1].shape[2]
144      cdims  = '1 %d 1' % dof
145      dims   = '%d %d %d' % (numSteps, dof, bs)
146      stride = '1 1 1'
147      size   = '1 %d 1' % dof
148    else:
149      dof    = f[1].shape[0]
150      bs     = f[1].shape[1]
151      cdims  = '%d 1' % dof
152      dims   = '%d %d' % (dof, bs)
153      stride = '1 1'
154      size   = '%d 1' % dof
155    for c in range(bs):
156      ext = self.typeExt[spaceDim][vtype][c]
157      if len(f[1].shape) > 2: start  = '%d 0 %d' % (timestep, c)
158      else:                   start  = '0 %d' % c
159      fp.write('''\
160        <Attribute
161           Name="%s"
162           Type="Scalar"
163           Center="%s">
164          <DataItem ItemType="HyperSlab"
165        	    Dimensions="%s"
166        	    Type="HyperSlab">
167            <DataItem
168               Dimensions="3 %d"
169               Format="XML">
170              %s
171              %s
172              %s
173            </DataItem>
174            <DataItem
175               DataType="Float" Precision="8"
176               Dimensions="%s"
177               Format="HDF">
178              &HeavyData;:%s
179            </DataItem>
180          </DataItem>
181        </Attribute>
182''' % (f[0]+'_'+ext, domain, cdims, len(f[1].shape), start, stride, size, dims, name))
183    return
184
185  def writeField(self, fp, numSteps, timestep, cellDim, spaceDim, name, f, domain):
186    ctypes = ['tensor', 'matrix']
187    if spaceDim == 2 or cellDim != spaceDim: ctypes.append('vector')
188    if f[1].attrs['vector_field_type'] in ctypes:
189      self.writeFieldComponents(fp, numSteps, timestep, spaceDim, name, f, domain)
190    else:
191      self.writeFieldSingle(fp, numSteps, timestep, spaceDim, name, f, domain)
192    return
193
194  def writeSpaceGridFooter(self, fp):
195    fp.write('      </Grid>\n')
196    return
197
198  def writeParticleGridHeader(self, fp, numParticles, spaceDim):
199    fp.write('''\
200      <Grid Name="particle_domain" GridType="Uniform">
201        <Topology TopologyType="Polyvertex" NodesPerElement="%d" />
202        <Geometry GeometryType="%s">
203          <DataItem Reference="XML">/Xdmf/Domain/DataItem[@Name="particle_coordinates"]</DataItem>
204        </Geometry>
205''' % (numParticles, "XYZ" if spaceDim > 2 else "XY"))
206    return
207
208  def writeParticleField(self, fp, fieldname, numParticles, numComp):
209    fp.write('''\
210    <Attribute Name="particles/%s">
211      <DataItem Name="%s"
212                Format="HDF"
213                Dimensions="%d %d">
214                &HeavyData;:/particle_fields/%s
215      </DataItem>
216    </Attribute>
217''' % (fieldname, fieldname, numParticles, numComp, fieldname))
218    return
219
220  def writeTimeGridFooter(self, fp):
221    fp.write('    </Grid>\n')
222    return
223
224  def writeFooter(self, fp):
225    fp.write('  </Domain>\n</Xdmf>\n')
226    return
227
228  def write(self, hdfFilename, topologyPath, numCells, numCorners, cellDim, htopologyPath, numHCells, numHCorners, geometryPath, numVertices, spaceDim, time, vfields, cfields, numParticles, pfields):
229    useTime = not (len(time) < 2 and time[0] == -1)
230    with open(self.filename, 'w') as fp:
231      self.writeHeader(fp, hdfFilename)
232      # Field information
233      self.writeCells(fp, topologyPath, numCells, numCorners)
234      if numHCells:
235        self.writeCells(fp, htopologyPath, numHCells, numHCorners, "hcells")
236      self.writeVertices(fp, geometryPath, numVertices, spaceDim)
237      if useTime: self.writeTimeGridHeader(fp, time)
238      for t in range(len(time)):
239        if numHCells:
240          self.writeHybridSpaceGridHeader(fp)
241          self.writeSpaceGridHeader(fp, numHCells, numHCorners, cellDim, spaceDim, "hcells")
242          self.writeSpaceGridFooter(fp)
243        self.writeSpaceGridHeader(fp, numCells, numCorners, cellDim, spaceDim)
244        for vf in vfields: self.writeField(fp, len(time), t, cellDim, spaceDim, '/vertex_fields/'+vf[0], vf, 'Node')
245        for cf in cfields: self.writeField(fp, len(time), t, cellDim, spaceDim, '/cell_fields/'+cf[0], cf, 'Cell')
246        self.writeSpaceGridFooter(fp)
247        if numHCells:
248          self.writeSpaceGridFooter(fp)
249      if useTime: self.writeTimeGridFooter(fp)
250      if numParticles:
251        self.writeLocations(fp, numParticles, spaceDim)
252        if useTime: self.writeTimeGridHeader(fp, time)
253        for t in range(len(time)):
254          self.writeParticleGridHeader(fp, numParticles, spaceDim)
255          for pf in pfields:
256            self.writeParticleField(fp, pf[0], numParticles, int(pf[1].attrs['Nc']))
257          self.writeSpaceGridFooter(fp)
258        if useTime: self.writeTimeGridFooter(fp)
259      self.writeFooter(fp)
260    return
261
262def generateXdmf(hdfFilename, xdmfFilename = None):
263  if xdmfFilename is None:
264    xdmfFilename = os.path.splitext(hdfFilename)[0] + '.xmf'
265  # Read mesh
266  h5          = h5py.File(hdfFilename, 'r')
267  if 'viz' in h5 and 'geometry' in h5['viz']:
268    geomPath  = 'viz/geometry'
269    geom      = h5['viz']['geometry']
270  else:
271    geomPath  = 'geometry'
272    geom      = h5['geometry']
273  if 'viz' in h5 and 'topology' in h5['viz']:
274    topoPath  = 'viz/topology'
275    topo      = h5['viz']['topology']
276  else:
277    topoPath  = 'topology'
278    topo      = h5['topology']
279  if 'viz' in h5 and 'hybrid_topology' in h5['viz']:
280    htopoPath = 'viz/hybrid_topology'
281    htopo     = h5['viz']['hybrid_topology']
282  else:
283    htopoPath = None
284    htopo     = None
285  vertices    = geom['vertices']
286  numVertices = vertices.shape[0]
287  spaceDim    = vertices.shape[1]
288  cells       = topo['cells']
289  numCells    = cells.shape[0]
290  numCorners  = cells.shape[1]
291  cellDim     = topo['cells'].attrs['cell_dim']
292  if htopo:
293    hcells      = htopo['cells']
294    numHCells   = hcells.shape[0]
295    numHCorners = hcells.shape[1]
296  else:
297    numHCells   = 0
298    numHCorners = 0
299  if 'time' in h5:
300    time      = np.array(h5['time']).flatten()
301  else:
302    time      = [-1]
303  vfields     = []
304  cfields     = []
305  pfields     = []
306  pfields     = []
307  if 'vertex_fields' in h5: vfields = h5['vertex_fields'].items()
308  if 'cell_fields' in h5: cfields = h5['cell_fields'].items()
309  numParticles = 0
310  if 'particles' in h5:
311    numParticles = h5['particles']['coordinates'].shape[0]
312  if 'particle_fields' in h5: pfields = h5['particle_fields'].items()
313
314  # Write Xdmf
315  Xdmf(xdmfFilename).write(hdfFilename, topoPath, numCells, numCorners, cellDim, htopoPath, numHCells, numHCorners, geomPath, numVertices, spaceDim, time, vfields, cfields, numParticles, pfields)
316  h5.close()
317  return
318
319if __name__ == '__main__':
320  for f in sys.argv[1:]:
321    generateXdmf(f)
322