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