CDMATH/tests/swig/BoySurface/VTKnewReader.py

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