src/MEDLoader/Swig/VTKReader.py

   1 #  -*- coding: iso-8859-1 -*-
   2 # Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20 # Author Anthony GEAY (CEA/DEN/DM2S/STMF)
  21
  22 from MEDLoader import *
  23
  24 class PVDReader:
  25     @classmethod
  26     def New(cls,fileName):
  27         """ Static constructor. """
  28         return PVDReader(fileName)
  29         pass
  30
  31     def __init__(self,fileName):
  32         self._fileName=fileName
  33         pass
  34
  35     def loadTopInfo(self):
  36         fd=open(self._fileName,"r")
  37         return self.__parseXML(fd)
  38
  39     def __parseXML(self,fd):
  40         import xml.sax
  41         class PVD_SAX_Reader(xml.sax.ContentHandler):
  42             def __init__(self):
  43                 self._tsteps=[]
  44                 pass
  45             def startElement(self,name,attrs):
  46                 if name=="VTKFile":
  47                     if attrs["type"]!="Collection":
  48                         raise Exception("Mismatch between reader (PVD) type and file content !")
  49                     return
  50                 if name=="DataSet":
  51                     self._tsteps.append((float(attrs["timestep"]),str(attrs["file"])))
  52                     return
  53                 pass
  54             pass
  55         rd=PVD_SAX_Reader()
  56         parser=xml.sax.make_parser()
  57         parser.setContentHandler(rd)
  58         parser.parse(fd)
  59         return rd
  60     pass
  61
  62 class PVTUReader:
  63     @classmethod
  64     def New(cls,fileName):
  65         """ Static constructor. """
  66         return PVTUReader(fileName)
  67         pass
  68
  69     def __init__(self,fileName):
  70         self._fileName=fileName
  71         pass
  72
  73     def loadParaInfo(self):
  74         fd=open(self._fileName,"r")
  75         return self.__parseXML(fd)
  76
  77     def __parseXML(self,fd):
  78         import xml.sax
  79         class PVTU_SAX_Reader(xml.sax.ContentHandler):
  80             def __init__(self):
  81                 self._data_array={2:self.DAPointData,3:self.DACellData}
  82                 self._node_fields=[]
  83                 self._cell_fields=[]
  84                 self._pfiles=[]
  85                 self._tmp=None
  86                 pass
  87             def DAPointData(self,attrs):
  88                 self._node_fields.append((str(attrs["Name"]),int(attrs["NumberOfComponents"])))
  89                 pass
  90             def DACellData(self,attrs):
  91                 self._cell_fields.append((str(attrs["Name"]),int(attrs["NumberOfComponents"])))
  92                 pass
  93             def startElement(self,name,attrs):
  94                 if name=="VTKFile":
  95                     if attrs["type"]!="PUnstructuredGrid":
  96                         raise Exception("Mismatch between reader (PVTU) type and file content !")
  97                     return
  98                 if name=="Piece":
  99                     self._pfiles.append(str(attrs["Source"]))
 100                     return
 101                 if name=="PPointData":
 102                     self._tmp=2
 103                     return
 104                 if name=="PCellData":
 105                     self._tmp=3
 106                     return
 107                 if name=="PDataArray":
 108                     if self._tmp in self._data_array.keys():
 109                         self._data_array[self._tmp](attrs)
 110                         pass
 111                     return
 112                 pass
 113             pass
 114         rd=PVTU_SAX_Reader()
 115         parser=xml.sax.make_parser()
 116         parser.setContentHandler(rd)
 117         parser.parse(fd)
 118         return rd
 119     pass
 120
 121 class VTURawReader:
 122     """ Converting a VTU file in raw mode into the MED format.
 123     """
 124     VTKTypes_2_MC=[-1,0,-1,1,33,3,-1,5,-1,4,14,-1,NORM_HEXA8,16,15,-1,22,-1,-1,-1,-1,2,6,8,20,30,25,23,9,27,-1,-1,-1,-1,7,-1,-1,-1,-1,-1,-1,-1,31]
 125
 126     class NormalException(Exception):
 127         pass
 128
 129     class NotRawVTUException(Exception):
 130         pass
 131
 132     def loadInMEDFileDS(self):
 133         import numpy as np
 134         fd=open(self._fileName,"r")
 135         ref,rd=self.__parseXML(fd)
 136         #
 137         ret=MEDFileData()
 138         ms=MEDFileMeshes() ; ret.setMeshes(ms)
 139         fs=MEDFileFields() ; ret.setFields(fs)
 140         #
 141         types=np.memmap(fd,dtype=rd._type_types,mode='r',offset=ref+rd._off_types,shape=(rd._nb_cells,))
 142         types=self.__swapIfNecessary(rd._bo,types)
 143         # mesh dimension detection
 144         types2=types.copy() ; types2.sort() ; types2=np.unique(types2)
 145         meshDim=MEDCouplingMesh.GetDimensionOfGeometricType(self.VTKTypes_2_MC[types2[0]])
 146         for typ in types2[1:]:
 147             md=MEDCouplingMesh.GetDimensionOfGeometricType(self.VTKTypes_2_MC[typ])
 148             if md!=meshDim:
 149                 raise Exception("MultiLevel umeshes not managed yet !")
 150             pass
 151         m=MEDCouplingUMesh("mesh",meshDim)
 152         # coordinates
 153         coo=np.memmap(fd,dtype=rd._type_coords,mode='r',offset=ref+rd._off_coords,shape=(rd._nb_nodes*rd._space_dim,))
 154         coo=self.__swapIfNecessary(rd._bo,coo) ; coo=DataArrayDouble(np.array(coo,dtype='float64')) ; coo.rearrange(rd._space_dim)
 155         m.setCoords(coo)
 156         # connectivity
 157         offsets=np.memmap(fd,dtype=rd._type_off,mode='r',offset=ref+rd._off_off,shape=(rd._nb_cells,))
 158         offsets=self.__swapIfNecessary(rd._bo,offsets) ; connLgth=offsets[-1] ; offsets2=DataArrayInt(rd._nb_cells+1) ; offsets2.setIJ(0,0,0)
 159         offsets2[1:]=DataArrayInt(offsets)
 160         offsets3=offsets2.deltaShiftIndex() ; offsets2=offsets3.deepCpy() ; offsets3+=1 ; offsets3.computeOffsets2()
 161         offsets=offsets3
 162         tmp1=DataArrayInt(len(offsets2),2) ; tmp1[:,0]=1 ; tmp1[:,1]=offsets2 ; tmp1.rearrange(1) ; tmp1.computeOffsets2()
 163         tmp1=DataArrayInt.Range(1,2*len(offsets2),2).buildExplicitArrByRanges(tmp1)
 164         conn=np.memmap(fd,dtype=rd._type_conn,mode='r',offset=ref+rd._off_conn,shape=(connLgth,))
 165         conn=self.__swapIfNecessary(rd._bo,conn)
 166         types=np.array(types,dtype='int32') ; types=DataArrayInt(types) ; types.transformWithIndArr(self.VTKTypes_2_MC)
 167         conn2=DataArrayInt(offsets.back())
 168         conn2[offsets[0:-1]]=types
 169         conn2[tmp1]=DataArrayInt(conn)
 170         m.setConnectivity(conn2,offsets,True)
 171         m.checkCoherency() ; mm=MEDFileUMesh() ; mm.setMeshAtLevel(0,m) ; ms.pushMesh(mm)
 172         # Fields on nodes and on cells
 173         for spatialDisc,nbEnt,fields in [(ON_NODES,rd._nb_nodes,rd._node_fields),(ON_CELLS,rd._nb_cells,rd._cell_fields)]:
 174             for name,typ,nbCompo,off in fields:
 175                 ff=MEDFileFieldMultiTS()
 176                 f=MEDCouplingFieldDouble(spatialDisc,ONE_TIME)
 177                 f.setName(name) ; f.setMesh(m)
 178                 vals=np.memmap(fd,dtype=typ,mode='r',offset=ref+off,shape=(nbEnt*nbCompo))
 179                 arr=DataArrayDouble(np.array(vals,dtype='float64')) ; arr.rearrange(nbCompo)
 180                 f.setArray(arr) ; f.checkCoherency()
 181                 f.setTime(self._time[0],self._time[1],0)
 182                 ff.appendFieldNoProfileSBT(f)
 183                 fs.pushField(ff)
 184                 pass
 185             pass
 186         return ret
 187
 188     def __parseXML(self,fd):
 189         import xml.sax
 190         class VTU_SAX_Reader(xml.sax.ContentHandler):
 191             def __init__(self):
 192                 self._data_array={0:self.DAPoints,1:self.DACells,2:self.DAPointData,3:self.DACellData}
 193                 self._node_fields=[]
 194                 self._cell_fields=[]
 195                 pass
 196             def DAPoints(self,attrs):
 197                 self._space_dim=int(attrs["NumberOfComponents"])
 198                 self._type_coords=str(attrs["type"]).lower()
 199                 self._off_coords=int(attrs["offset"])
 200                 pass
 201             def DACells(self,attrs):
 202                 if attrs["Name"]=="connectivity":
 203                     self._type_conn=str(attrs["type"]).lower()
 204                     self._off_conn=int(attrs["offset"])
 205                     pass
 206                 if attrs["Name"]=="offsets":
 207                     self._type_off=str(attrs["type"]).lower()
 208                     self._off_off=int(attrs["offset"])
 209                     pass
 210                 if attrs["Name"]=="types":
 211                     self._type_types=str(attrs["type"]).lower()
 212                     self._off_types=int(attrs["offset"])
 213                     pass
 214                 pass
 215             def DAPointData(self,attrs):
 216                 self._node_fields.append((str(attrs["Name"]),str(attrs["type"]).lower(),int(attrs["NumberOfComponents"]),int(attrs["offset"])))
 217                 pass
 218             def DACellData(self,attrs):
 219                 self._cell_fields.append((str(attrs["Name"]),str(attrs["type"]).lower(),int(attrs["NumberOfComponents"]),int(attrs["offset"])))
 220                 pass
 221             def startElement(self,name,attrs):
 222                 if name=="VTKFile":
 223                     if attrs["type"]!="UnstructuredGrid":
 224                         raise Exception("Mismatch between reader (VTU) type and file content !")
 225                     self._bo=bool(["LittleEndian","BigEndian"].index(attrs["byte_order"]))
 226                     pass
 227                 if name=="Piece":
 228                     self._nb_cells=int(attrs["NumberOfCells"])
 229                     self._nb_nodes=int(attrs["NumberOfPoints"])
 230                     return
 231                 if name=="Points":
 232                     self._tmp=0
 233                     return
 234                 if name=="Cells":
 235                     self._tmp=1
 236                     return
 237                 if name=="PointData":
 238                     self._tmp=2
 239                     return
 240                 if name=="CellData":
 241                     self._tmp=3
 242                     return
 243                 if name=="DataArray":
 244                     self._data_array[self._tmp](attrs)
 245                     return
 246                 if name=="AppendedData":
 247                     if str(attrs["encoding"])=="raw":
 248                         raise VTURawReader.NormalException("")
 249                     else:
 250                         raise VTURawReader.NotRawVTUException("The file is not a raw VTU ! Change reader !")
 251                 pass
 252             pass
 253         rd=VTU_SAX_Reader()
 254         parser=xml.sax.make_parser()
 255         parser.setContentHandler(rd)
 256         isOK=False
 257         try:
 258             parser.parse(fd)
 259         except self.NormalException as e:
 260             isOK=True
 261             fd.seek(0)
 262             for i in xrange(31): fd.readline()
 263             ref=fd.tell()+5
 264             pass
 265         if not isOK:
 266             raise Exception("Error in VTURawReader : not a raw format ?")
 267         return ref,rd
 268
 269     @classmethod
 270     def New(cls,fileName,tim=(0.,0)):
 271         """ Static constructor. """
 272         return VTURawReader(fileName,tim)
 273         pass
 274
 275     def __init__(self,fileName,tim=(0.,0)):
 276         msg="The time specified in constructor as 2nd arg should be a tuple containing 2 values 1 float and 1 int !"
 277         if type(tim)!=tuple:
 278             raise Exception(msg)
 279         if len(tim)!=2:
 280             raise Exception(msg)
 281         if type(tim[0])!=float or type(tim[1])!=int:
 282             raise Exception(msg)
 283         self._fileName=fileName
 284         self._time=tim
 285         pass
 286
 287     def __swapIfNecessary(self,b,arr):
 288         if b:
 289             ret=arr.copy()
 290             ret.byteswap(True)
 291             return ret
 292         else:
 293             return arr
 294         pass
 295     pass