src/MEDLoader/Swig/CaseReader.py

   1 #  -*- coding: iso-8859-1 -*-
   2 # Copyright (C) 2007-2015  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/LGLS)
  21
  22 # http://www-vis.lbl.gov/NERSC/Software/ensight/doc/OnlineHelp/UM-C11.pdf
  23 import numpy as np
  24 from MEDLoader import *
  25 from CaseIO import CaseIO
  26 import sys,re
  27
  28 class CaseReader(CaseIO):
  29     """ Converting a file in the Case format (Ensight) to the MED format.
  30     A new file with the same base name and the .med extension is created.
  31     """
  32
  33     @classmethod
  34     def New(cls,fileName):
  35         """ Static constructor. """
  36         return CaseReader(fileName)
  37         pass
  38
  39     def __init__(self,fileName):
  40         """ Constructor """
  41         self._fileName=fileName
  42         self._dirName=os.path.dirname(self._fileName)
  43         pass
  44
  45     def __traduceMesh(self,name,typ,coords,cells):
  46         """ Convert a CASE mesh into a MEDCouplingUMesh. """
  47         nbCoords=len(coords)
  48         coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
  49         coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
  50         ct=self.dictMCTyp2[typ]
  51         m=MEDCouplingUMesh(name,MEDCouplingUMesh.GetDimensionOfGeometricType(ct))
  52         m.setCoords(coo)
  53         nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(ct)
  54         cI=DataArrayInt(len(cells)+1) ; cI.iota() ; cI*=nbNodesPerCell+1
  55         #
  56         cells2=cells.reshape(len(cells),nbNodesPerCell)
  57         if cells2.dtype=='int32':
  58             c2=DataArrayInt(cells2)
  59         else:
  60             c2=DataArrayInt(np.array(cells2,dtype="int32"))
  61             pass
  62         c=DataArrayInt(len(cells),nbNodesPerCell+1) ; c[:,0]=ct ; c[:,1:]=c2-1 ; c.rearrange(1)
  63         m.setConnectivity(c,cI,True)
  64         m.checkCoherency2()
  65         return m
  66
  67     def __traduceMeshForPolyhed(self,name,coords,arr0,arr1,arr2):
  68         nbCoords=len(coords)
  69         coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
  70         coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
  71         m=MEDCouplingUMesh(name,3)
  72         m.setCoords(coo)
  73         #
  74         arr2=arr2[:]-1
  75         arr0mc0=DataArrayInt(arr0) ; arr0mc0.computeOffsets2()
  76         arr0mc1=DataArrayInt(arr0).deepCpy()
  77         arr0mc2=DataArrayInt(len(arr0),2) ; arr0mc2[:,0]=DataArrayInt(arr0)-1 ; arr0mc2[:,1]=1 ; arr0mc2.rearrange(1) ; arr0mc2.computeOffsets2()
  78         arr0mc3=DataArrayInt.Range(0,2*len(arr0),2).buildExplicitArrByRanges(arr0mc2)
  79         arr1mc0=DataArrayInt(arr1) ; arr1mc0.computeOffsets2()
  80         arr1mc1=arr1mc0[arr0mc0] ; arr1mc1[1:]+=arr0mc0[1:]
  81         arr1mc2=DataArrayInt(arr1).deepCpy() ; arr1mc2+=1 ; arr1mc2.computeOffsets2()
  82         arr2mc0=(arr1mc2[1:])[arr0mc3]
  83         #
  84         c=DataArrayInt(arr1.size+arr2.size)
  85         c[arr1mc1[:-1]]=NORM_POLYHED
  86         c[arr2mc0]=-1
  87         a=arr2mc0.buildUnion(arr1mc1[:-1]).buildComplement(len(c))
  88         c[a]=DataArrayInt(arr2)
  89         #
  90         m.setConnectivity(c,arr1mc1,True)
  91         m.checkCoherency2()
  92         return m
  93
  94     def __traduceMeshForPolygon(self,name,coords,arr0,arr1):
  95         nbCoords=len(coords)
  96         coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
  97         coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
  98         m=MEDCouplingUMesh(name,2)
  99         m.setCoords(coo)
 100         #
 101         arr0_0=DataArrayInt(arr0+1) ; arr0_0.computeOffsets2()
 102         arr0_1=DataArrayInt(len(arr0),2) ; arr0_1[:,1]=DataArrayInt(arr0) ; arr0_1[:,0]=1 ; arr0_1.rearrange(1) ; arr0_1.computeOffsets2()
 103         arr0_2=DataArrayInt.Range(1,2*len(arr0),2).buildExplicitArrByRanges(arr0_1)
 104         c=DataArrayInt(len(arr0)+len(arr1)) ; c[:]=0 ; c[arr0_0[:-1]]=NORM_POLYGON
 105         c[arr0_2]=DataArrayInt(arr1-1)
 106         #
 107         m.setConnectivity(c,arr0_0,True)
 108         m.checkCoherency2()
 109         return m
 110
 111     def __convertGeo2MED(self,geoFileName):
 112         """ Convert all the geometry (all the meshes) contained in teh CASE file into MEDCouplingUMesh'es. """
 113         fd=open(os.path.join(self._dirName,geoFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0) ; title=fd.read(80)
 114         title=title.strip().lower()
 115         if "binary" not in title:
 116             raise Exception("Error only binary geo files are supported for the moment !")
 117             pass
 118         zeType=True
 119         if "fortran" in title:
 120             mcmeshes=self.__convertGeo2MEDFortran(fd,end) ; zeType=False
 121         else:
 122             mcmeshes=self.__convertGeo2MEDC(fd,end)
 123         #
 124         ms=MEDFileMeshes()
 125         ms.resize(len(mcmeshes))
 126         for i,m in enumerate(mcmeshes):
 127             mlm=MEDFileUMesh()
 128             mlm.setMeshAtLevel(0,m)
 129             ms.setMeshAtPos(i,mlm)
 130             pass
 131         return mcmeshes,ms,zeType
 132
 133     def __convertGeo2MEDFortran(self,fd,end):
 134         mcmeshes=[]
 135         fd.read(80) # comment 1
 136         fd.read(80) # comment 2
 137         fd.read(80) # node id
 138         fd.read(80) # element id
 139         pos=fd.tell()
 140         elt=fd.read(80) ; elt=elt.strip() ; pos=fd.tell()
 141         mcmeshes2=[]
 142         typ="part"
 143         nbOfTurn=0
 144         while abs(pos-end)>8 and "part" in typ:
 145             if "part" not in elt:
 146                 raise Exception("Error on reading mesh fortran #1 !")
 147             fd.seek(fd.tell()+4)# skip #
 148             tmp=fd.read(80) ; meshName=tmp.split("P")[-1]
 149             tmp=fd.read(80)
 150             if "coordinates" not in tmp:
 151                 raise Exception("Error on reading mesh fortran #2 !")
 152             pos=fd.tell() # 644
 153             if nbOfTurn==0:
 154                 pos+=76 # what else ?
 155             else:
 156                 pos+=40
 157                 pass
 158             nbNodes=np.memmap(fd,dtype='>i4',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
 159             pos+=12 # what else ?
 160             a=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos),shape=(nbNodes))
 161             b=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos+nbNodes*4+2*4),shape=(nbNodes))
 162             c=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos+nbNodes*2*4+4*4),shape=(nbNodes))
 163             coo=np.zeros(dtype=">f4",shape=(nbNodes*3))
 164             coo[:nbNodes]=a ; coo[nbNodes:2*nbNodes]=b ; coo[2*nbNodes:]=c
 165             coo=coo.reshape(nbNodes,3)
 166             pos+=nbNodes*3*4 ; fd.seek(pos)#np.array(0,dtype='float%i'%(typeOfCoo)).nbytes
 167             typ=fd.read(80).strip() ; pos=fd.tell()
 168             zeK=""
 169             for k in self.dictMCTyp2.keys():
 170                 if k in typ:
 171                     zeK=k
 172                     break
 173                     pass
 174                 pass
 175             pos+=8*4 # yeh man !
 176             nbCellsOfType=np.memmap(fd,dtype='>i4',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
 177             pos+=4 # for the number of cells
 178             pos+=2*4 # because it's great !
 179             nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(self.dictMCTyp2[zeK])
 180             nodalConn=np.memmap(fd,dtype='>i4',mode='r',offset=pos,shape=(nbCellsOfType,nbNodesPerCell))
 181             meshName=meshName.strip()
 182             mcmeshes2.append(self.__traduceMesh(meshName,zeK,coo,nodalConn))
 183             pos+=nbNodesPerCell*nbCellsOfType*4
 184             if abs(pos-end)>8:
 185                 fd.seek(pos) ;elt=fd.read(80) ; typ=elt[:] ; pos+=80
 186                 pass
 187             nbOfTurn+=1
 188             pass
 189         #coo=mcmeshes2[0].getCoords() ; name=mcmeshes2[0].getName()
 190         #for itmesh in mcmeshes2: itmesh.setCoords(coo)
 191         #m=MEDCouplingUMesh.MergeUMeshesOnSameCoords(mcmeshes2) ; m.setName(name)
 192         #mcmeshes.append(m)
 193         return mcmeshes2
 194
 195     def __convertGeo2MEDC(self,fd,end):
 196         fd.readline()
 197         name=fd.readline().strip() ; fd.readline() ; fd.readline()
 198         pos=fd.tell()
 199         mcmeshes=[]
 200         elt=fd.read(80) ; elt=elt.strip() ; pos+=80
 201         while pos!=end:
 202             if "part" not in elt:
 203                 raise Exception("Error on reading mesh #1 !")
 204             fd.seek(fd.tell()+4)
 205             meshName=fd.read(80).strip()
 206             if fd.read(len("coordinates"))!="coordinates":
 207                 raise Exception("Error on reading mesh #2 !")
 208             pos=fd.tell()
 209             typeOfCoo=np.memmap(fd,dtype='byte',mode='r',offset=int(pos),shape=(1)).tolist()[0]
 210             pos+=1+17*4
 211             nbNodes=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
 212             pos+=4
 213             coo=np.memmap(fd,dtype='float32',mode='r',offset=int(pos),shape=(nbNodes,3))
 214             pos+=nbNodes*3*4 ; fd.seek(pos)#np.array(0,dtype='float%i'%(typeOfCoo)).nbytes
 215             typ=fd.read(80).strip() ; pos=fd.tell()
 216             mcmeshes2=[]
 217             while pos!=end and typ!="part":
 218                 mctyp=self.dictMCTyp2[typ]
 219                 nbCellsOfType=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
 220                 pos+=4
 221                 if mctyp!=NORM_POLYHED and mctyp!=NORM_POLYGON:
 222                     nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(mctyp)
 223                     cells=np.memmap(fd,dtype='int32',mode='r',offset=pos,shape=(nbCellsOfType,nbNodesPerCell))
 224                     pos+=nbCellsOfType*nbNodesPerCell*4
 225                     fd.seek(pos)
 226                     mcmeshes2.append(self.__traduceMesh(meshName,typ,coo,cells))
 227                 elif mctyp==NORM_POLYHED:
 228                     nbOfFacesPerCell=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(nbCellsOfType,))
 229                     pos+=nbCellsOfType*4
 230                     szOfNbOfNodesPerFacePerCellArr=int(nbOfFacesPerCell.sum())
 231                     arr1=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNbOfNodesPerFacePerCellArr,))#arr1 -> nbOfNodesPerFacePerCellArr
 232                     pos+=szOfNbOfNodesPerFacePerCellArr*4
 233                     szOfNodesPerFacePerCellArr=arr1.sum()
 234                     arr2=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNodesPerFacePerCellArr,))#arr2 -> nodesPerFacePerCellArr
 235                     pos+=szOfNodesPerFacePerCellArr*4 ; fd.seek(pos)
 236                     mcmeshes2.append(self.__traduceMeshForPolyhed(meshName,coo,nbOfFacesPerCell,arr1,arr2))
 237                     pass
 238                 else:
 239                     nbOfNodesPerCell=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(nbCellsOfType,))
 240                     pos+=nbCellsOfType*4
 241                     szOfNbOfNodesPerCellArr=int(nbOfNodesPerCell.sum())
 242                     arr1=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNbOfNodesPerCellArr,))
 243                     pos+=szOfNbOfNodesPerCellArr*4  ; fd.seek(pos)
 244                     mcmeshes2.append(self.__traduceMeshForPolygon(meshName,coo,nbOfNodesPerCell,arr1))
 245                 if pos!=end:
 246                     elt=fd.read(80) ; elt=elt.strip() ; typ=elt[:] ; pos+=80
 247                     pass
 248                 pass
 249             coo=mcmeshes2[0].getCoords() ; name=mcmeshes2[0].getName()
 250             for itmesh in mcmeshes2: itmesh.setCoords(coo)
 251             m=MEDCouplingUMesh.MergeUMeshesOnSameCoords(mcmeshes2) ; m.setName(name)
 252             mcmeshes.append(m)
 253             pass
 254         return mcmeshes
 255
 256
 257     def __convertField(self,mlfields, mcmeshes, fileName, fieldName, discr, nbCompo, locId, it):
 258         """ Convert the fields. """
 259         stars=re.search("[\*]+",fileName).group()
 260         st="%0"+str(len(stars))+"i"
 261         trueFileName=fileName.replace(stars,st%(it))
 262         fd=open(os.path.join(self._dirName,trueFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0)
 263         name=fd.readline().strip().split(" ")[0]
 264         if name!=fieldName:
 265             raise Exception("ConvertField : mismatch")
 266         pos=fd.tell()
 267         st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
 268         while pos!=end:
 269             if st!="part":
 270                 raise Exception("ConvertField : mismatch #2")
 271             fdisc=MEDCouplingFieldDiscretization.New(self.discSpatial2[discr])
 272             meshId=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1)).tolist()[0]-1
 273             nbOfValues=fdisc.getNumberOfTuples(mcmeshes[meshId])
 274             vals2=DataArrayDouble(nbOfValues,nbCompo)
 275             fd.seek(pos+4)
 276             st=fd.read(80).strip() ; pos=fd.tell()
 277             offset=0
 278             while pos!=end and st!="part":
 279                 if st!="coordinates":
 280                     nbOfValsOfTyp=mcmeshes[meshId].getNumberOfCellsWithType(self.dictMCTyp2[st])
 281                 else:
 282                     nbOfValsOfTyp=nbOfValues
 283                     pass
 284                 vals=np.memmap(fd,dtype='float32',mode='r',offset=int(pos),shape=(nbOfValsOfTyp,nbCompo))#np.memmap(fd,dtype='int32',mode='r',offset=159,shape=(1))
 285                 vals2[offset:offset+nbOfValsOfTyp]=DataArrayDouble(np.array(vals,dtype='float64')).fromNoInterlace()
 286                 pos+=nbOfValsOfTyp*nbCompo*4 ; fd.seek(pos)
 287                 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
 288                 offset+=nbOfValsOfTyp
 289                 pass
 290             f=MEDCouplingFieldDouble(self.discSpatial2[discr],ONE_TIME) ; f.setName("%s_%s"%(fieldName,mcmeshes[meshId].getName()))
 291             f.setMesh(mcmeshes[meshId]) ; f.setArray(vals2) ; f.setTime(float(it),it,-1)
 292             f.checkCoherency()
 293             mlfields[locId+meshId].appendFieldNoProfileSBT(f)
 294             pass
 295
 296     def __convertFieldFortran(self,mlfields, mcmeshes, fileName, fieldName, discr, nbCompo, locId, it):
 297         """ Convert the fields. """
 298         if re.search("[\*]+",fileName):
 299             stars=re.search("[\*]+",fileName).group()
 300             st="%0"+str(len(stars))+"i"
 301             trueFileName=fileName.replace(stars,st%(it))
 302             pass
 303         else:
 304             trueFileName=fileName
 305             pass
 306         fd=open(os.path.join(self._dirName,trueFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0)
 307         name=fd.read(80)
 308         if fieldName not in name:
 309             raise Exception("ConvertField : mismatch")
 310         pos=fd.tell()
 311         st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
 312         if "part" not in st:
 313             raise Exception("ConvertField : mismatch #2")
 314         st=fd.read(80).strip() ; pos=fd.tell()
 315         pos+=12 # I love it
 316         offset=0
 317         nbTurn=0
 318         while pos!=end and "part" not in st:
 319             fdisc=MEDCouplingFieldDiscretization.New(self.discSpatial2[discr])
 320             nbOfValues=fdisc.getNumberOfTuples(mcmeshes[nbTurn])
 321             vals2=DataArrayDouble(nbOfValues,nbCompo)
 322             pos+=24 # I love it again !
 323             nbOfValsOfTyp=np.memmap(fd,dtype='>i4',mode='r',offset=pos,shape=(1)).tolist()[0]/4
 324             pos+=4
 325             vals=np.zeros(dtype=">f4",shape=(nbOfValsOfTyp*nbCompo))
 326             for iii in xrange(nbCompo):
 327                 valsTmp=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos),shape=(nbOfValsOfTyp))
 328                 vals[iii*nbOfValsOfTyp:(iii+1)*nbOfValsOfTyp]=valsTmp
 329                 pos+=nbOfValsOfTyp*4
 330                 pos+=2*4 ## hey hey, that is the ultimate class !
 331                 vals2.setInfoOnComponent(iii,chr(ord('X')+iii))
 332                 pass
 333             if pos>end:
 334                 pos=end
 335                 pass
 336             vals=vals.reshape(nbOfValsOfTyp,nbCompo)
 337             vals2[offset:offset+nbOfValsOfTyp]=DataArrayDouble(np.array(vals,dtype='float64')).fromNoInterlace()
 338             if pos!=end:
 339                 fd.seek(pos)
 340                 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
 341                 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
 342                 pass
 343             f=MEDCouplingFieldDouble(self.discSpatial2[discr],ONE_TIME) ; f.setName("%s_%s"%(fieldName,mcmeshes[nbTurn].getName()))
 344             f.setMesh(mcmeshes[nbTurn]) ; f.setArray(vals2) ; f.setTime(float(it),it,-1)
 345             f.checkCoherency()
 346             mlfields[locId+nbTurn].appendFieldNoProfileSBT(f)
 347             nbTurn+=1
 348             pass
 349         pass
 350
 351     def loadInMEDFileDS(self):
 352         """ Load a CASE file into a MEDFileData object. """
 353         f=file(self._fileName)
 354         lines=f.readlines()
 355         ind=lines.index("GEOMETRY\n")
 356         if ind==-1:
 357             raise Exception("Error with file %s"%(fname))
 358         geoName=re.match("model:([\W]*)([\w\.]+)",lines[ind+1]).group(2)
 359         m1,m2,typeOfFile=self.__convertGeo2MED(geoName)
 360         fieldsInfo=[] ; nbOfTimeSteps=0
 361         if "VARIABLE\n" in lines:
 362             ind=lines.index("VARIABLE\n")
 363             end=len(lines)-1
 364             if "TIME\n" in lines:
 365                 end=lines.index("TIME\n")
 366                 pass
 367             for i in xrange(ind+1,end):
 368                 m=re.match("^([\w]+)[\s]+\per[\s]+([\w]+)[\s]*\:[\s]*([\w]+)[\s]+([\S]+)$",lines[i])
 369                 if m:
 370                     if m.groups()[0]=="constant":
 371                         continue
 372                     spatialDisc=m.groups()[1] ; fieldName=m.groups()[2] ; nbOfCompo=self.dictCompo2[m.groups()[0]] ; fieldFileName=m.groups()[3]
 373                     fieldsInfo.append((fieldName,spatialDisc,nbOfCompo,fieldFileName))
 374                     pass
 375                 pass
 376
 377             expr=re.compile("number[\s]+of[\s]+steps[\s]*\:[\s]*([\d]+)")
 378             tmp=filter(expr.search,lines)
 379             if len(tmp)!=0:
 380                 nbOfTimeSteps=int(expr.search(filter(expr.search,lines)[0]).group(1))
 381                 expr=re.compile("filename[\s]+start[\s]+number[\s]*\:[\s]*([\d]+)")
 382                 startIt=int(expr.search(filter(expr.search,lines)[0]).group(1))
 383                 expr=re.compile("filename[\s]+increment[\s]*\:[\s]*([\d]+)")
 384                 incrIt=int(expr.search(filter(expr.search,lines)[0]).group(1))
 385             else:
 386                 nbOfTimeSteps=1
 387                 startIt=0
 388                 incrIt=1
 389                 pass
 390             curIt=startIt
 391             pass
 392         mlfields=MEDFileFields()
 393         mlfields.resize(len(fieldsInfo)*len(m1))
 394         i=0
 395         for field in fieldsInfo:
 396             for m in m1:
 397                 mlfields.setFieldAtPos(i,MEDFileFieldMultiTS())
 398                 i+=1
 399                 pass
 400             pass
 401         for ts in xrange(nbOfTimeSteps):
 402             i=0
 403             for field in fieldsInfo:
 404                 if typeOfFile:
 405                     self.__convertField(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt);
 406                 else:
 407                     self.__convertFieldFortran(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt)
 408                     pass
 409                 i+=len(m1)
 410                 pass
 411             curIt+=incrIt
 412             pass
 413         ret=MEDFileData()
 414         ret.setMeshes(m2)
 415         del mlfields[filter(lambda x: len(mlfields[x])==0,range(len(mlfields)))]
 416         ret.setFields(mlfields)
 417         return ret
 418
 419     pass