src/MEDLoader/Swig/VTKReader.py

   1 #  -*- coding: iso-8859-1 -*-
   2 # Copyright (C) 2007-2020  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, or (at your option) any later version.
   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:
 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         def __init__(self,lineNb):
 128             Exception.__init__(self)
 129             self._line_nb=lineNb
 130         def getLineNb(self):
 131             return self._line_nb
 132         pass
 133
 134     class NotRawVTUException(Exception):
 135         pass
 136
 137     def loadInMEDFileDS(self):
 138         import numpy as np
 139         fd=open(self._fileName,"r")
 140         ref,rd=self.__parseXML(fd)
 141         #
 142         ret=MEDFileData()
 143         ms=MEDFileMeshes() ; ret.setMeshes(ms)
 144         fs=MEDFileFields() ; ret.setFields(fs)
 145         #
 146         types=np.memmap(fd,dtype=rd._type_types,mode='r',offset=ref+rd._off_types,shape=(rd._nb_cells,))
 147         types=self.__swapIfNecessary(rd._bo,types)
 148         # mesh dimension detection
 149         types2=types.copy() ; types2.sort() ; types2=np.unique(types2)
 150         meshDim=MEDCouplingMesh.GetDimensionOfGeometricType(self.VTKTypes_2_MC[types2[0]])
 151         for typ in types2[1:]:
 152             md=MEDCouplingMesh.GetDimensionOfGeometricType(self.VTKTypes_2_MC[typ])
 153             if md!=meshDim:
 154                 raise Exception("MultiLevel umeshes not managed yet !")
 155             pass
 156         m=MEDCouplingUMesh("mesh",meshDim)
 157         # coordinates
 158         coo=np.memmap(fd,dtype=rd._type_coords,mode='r',offset=ref+rd._off_coords,shape=(rd._nb_nodes*rd._space_dim,))
 159         coo=self.__swapIfNecessary(rd._bo,coo) ; coo=DataArrayDouble(np.array(coo,dtype='float64')) ; coo.rearrange(rd._space_dim)
 160         m.setCoords(coo)
 161         # connectivity
 162         offsets=np.memmap(fd,dtype=rd._type_off,mode='r',offset=ref+rd._off_off,shape=(rd._nb_cells,))
 163         offsets=self.__swapIfNecessary(rd._bo,offsets) ; connLgth=offsets[-1] ; offsets2=DataArrayInt(rd._nb_cells+1) ; offsets2.setIJ(0,0,0)
 164         offsets2[1:]=DataArrayInt(offsets)
 165         offsets3=offsets2.deltaShiftIndex() ; offsets2=offsets3.deepCopy() ; offsets3+=1 ; offsets3.computeOffsetsFull()
 166         offsets=offsets3
 167         tmp1=DataArrayInt(len(offsets2),2) ; tmp1[:,0]=1 ; tmp1[:,1]=offsets2 ; tmp1.rearrange(1) ; tmp1.computeOffsetsFull()
 168         tmp1=DataArrayInt.Range(1,2*len(offsets2),2).buildExplicitArrByRanges(tmp1)
 169         conn=np.memmap(fd,dtype=rd._type_conn,mode='r',offset=ref+rd._off_conn,shape=(connLgth,))
 170         conn=self.__swapIfNecessary(rd._bo,conn)
 171         types=np.array(types,dtype='int32') ; types=DataArrayInt(types) ; types.transformWithIndArr(self.VTKTypes_2_MC)
 172         conn2=DataArrayInt(offsets.back())
 173         conn2[offsets[0:-1]]=types
 174         conn2[tmp1]=DataArrayInt(conn)
 175         m.setConnectivity(conn2,offsets,True)
 176         m.checkConsistencyLight() ; mm=MEDFileUMesh() ; mm.setMeshAtLevel(0,m) ; ms.pushMesh(mm)
 177         # Fields on nodes and on cells
 178         for spatialDisc,nbEnt,fields in [(ON_NODES,rd._nb_nodes,rd._node_fields),(ON_CELLS,rd._nb_cells,rd._cell_fields)]:
 179             for name,typ,nbCompo,off in fields:
 180                 ff=MEDFileFieldMultiTS()
 181                 f=MEDCouplingFieldDouble(spatialDisc,ONE_TIME)
 182                 f.setName(name) ; f.setMesh(m)
 183                 vals=np.memmap(fd,dtype=typ,mode='r',offset=ref+off,shape=(nbEnt*nbCompo))
 184                 vals=self.__swapIfNecessary(rd._bo,vals)
 185                 arr=DataArrayDouble(np.array(vals,dtype='float64')) ; arr.rearrange(nbCompo)
 186                 f.setArray(arr) ; f.checkConsistencyLight()
 187                 f.setTime(self._time[0],self._time[1],0)
 188                 ff.appendFieldNoProfileSBT(f)
 189                 fs.pushField(ff)
 190                 pass
 191             pass
 192         return ret
 193
 194     def __parseXML(self,fd):
 195         import xml.sax
 196         class VTU_SAX_Reader(xml.sax.ContentHandler):
 197             def __init__(self):
 198                 self._loc=None
 199                 self._data_array={0:self.DAPoints,1:self.DACells,2:self.DAPointData,3:self.DACellData}
 200                 self._node_fields=[]
 201                 self._cell_fields=[]
 202                 pass
 203             def setLocator(self,loc):
 204                 self._loc=loc
 205             def DAPoints(self,attrs):
 206                 self._space_dim=int(attrs["NumberOfComponents"])
 207                 self._type_coords=str(attrs["type"]).lower()
 208                 self._off_coords=int(attrs["offset"])
 209                 pass
 210             def DACells(self,attrs):
 211                 if attrs["Name"]=="connectivity":
 212                     self._type_conn=str(attrs["type"]).lower()
 213                     self._off_conn=int(attrs["offset"])
 214                     pass
 215                 if attrs["Name"]=="offsets":
 216                     self._type_off=str(attrs["type"]).lower()
 217                     self._off_off=int(attrs["offset"])
 218                     pass
 219                 if attrs["Name"]=="types":
 220                     self._type_types=str(attrs["type"]).lower()
 221                     self._off_types=int(attrs["offset"])
 222                     pass
 223                 pass
 224             def DAPointData(self,attrs):
 225                 self._node_fields.append((str(attrs["Name"]),str(attrs["type"]).lower(),int(attrs["NumberOfComponents"]),int(attrs["offset"])))
 226                 pass
 227             def DACellData(self,attrs):
 228                 self._cell_fields.append((str(attrs["Name"]),str(attrs["type"]).lower(),int(attrs["NumberOfComponents"]),int(attrs["offset"])))
 229                 pass
 230             def startElement(self,name,attrs):
 231                 if name=="VTKFile":
 232                     if attrs["type"]!="UnstructuredGrid":
 233                         raise Exception("Mismatch between reader (VTU) type and file content !")
 234                     self._bo=bool(["LittleEndian","BigEndian"].index(attrs["byte_order"]))
 235                     pass
 236                 if name=="Piece":
 237                     self._nb_cells=int(attrs["NumberOfCells"])
 238                     self._nb_nodes=int(attrs["NumberOfPoints"])
 239                     return
 240                 if name=="Points":
 241                     self._tmp=0
 242                     return
 243                 if name=="Cells":
 244                     self._tmp=1
 245                     return
 246                 if name=="PointData":
 247                     self._tmp=2
 248                     return
 249                 if name=="CellData":
 250                     self._tmp=3
 251                     return
 252                 if name=="DataArray":
 253                     self._data_array[self._tmp](attrs)
 254                     return
 255                 if name=="AppendedData":
 256                     if str(attrs["encoding"])=="raw":
 257                         raise VTURawReader.NormalException(self._loc.getLineNumber())
 258                     else:
 259                         raise VTURawReader.NotRawVTUException("The file is not a raw VTU ! Change reader !")
 260                 pass
 261             pass
 262         rd=VTU_SAX_Reader()
 263         parser=xml.sax.make_parser()
 264         parser.setContentHandler(rd)
 265         locator=xml.sax.expatreader.ExpatLocator(parser)
 266         rd.setLocator(locator)
 267         isOK=False
 268         try:
 269             parser.parse(fd)
 270         except self.NormalException as e:
 271             isOK=True
 272             fd.seek(0)
 273             for i in range(e.getLineNb()): fd.readline()
 274             ref=fd.tell()+5
 275             pass
 276         if not isOK:
 277             raise Exception("Error in VTURawReader : not a raw format ?")
 278         return ref,rd
 279
 280     @classmethod
 281     def New(cls,fileName,tim=(0.,0)):
 282         """ Static constructor. """
 283         return VTURawReader(fileName,tim)
 284         pass
 285
 286     def __init__(self,fileName,tim=(0.,0)):
 287         msg="The time specified in constructor as 2nd arg should be a tuple containing 2 values 1 float and 1 int !"
 288         if not isinstance(tim, tuple):
 289             raise Exception(msg)
 290         if len(tim)!=2:
 291             raise Exception(msg)
 292         if not isinstance(tim[0], float) or not isinstance(tim[1], int):
 293             raise Exception(msg)
 294         self._fileName=fileName
 295         self._time=tim
 296         pass
 297
 298     def __swapIfNecessary(self,b,arr):
 299         if b:
 300             ret=arr.copy()
 301             ret.byteswap(True)
 302             return ret
 303         else:
 304             return arr
 305         pass
 306     pass