src/MEDLoader/Swig/CaseReader.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/LGLS)
  21
  22 import numpy as np
  23 from MEDLoader import *
  24 from CaseIO import CaseIO
  25 import sys,re
  26
  27 class CaseReader(CaseIO):
  28     """ Converting a file in the Case format (Ensight) to the MED format.
  29     A new file with the same base name and the .med extension is created.
  30     """
  31
  32     @classmethod
  33     def New(cls,fileName):
  34         """ Static constructor. """
  35         return CaseReader(fileName)
  36         pass
  37
  38     def __init__(self,fileName):
  39         """ Constructor """
  40         self._fileName=fileName
  41         self._dirName=os.path.dirname(self._fileName)
  42         pass
  43
  44     def __traduceMesh(self,name,typ,coords,cells):
  45         """ Convert a CASE mesh into a MEDCouplingUMesh. """
  46         nbCoords=len(coords)
  47         coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
  48         coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
  49         ct=self.dictMCTyp2[typ]
  50         m=MEDCouplingUMesh(name,MEDCouplingUMesh.GetDimensionOfGeometricType(ct))
  51         m.setCoords(coo)
  52         nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(ct)
  53         cI=DataArrayInt(len(cells)+1) ; cI.iota() ; cI*=nbNodesPerCell+1
  54         #
  55         cells2=cells.reshape(len(cells),nbNodesPerCell)
  56         c2=DataArrayInt(cells2)
  57         c=DataArrayInt(len(cells),nbNodesPerCell+1) ; c[:,0]=ct ; c[:,1:]=c2-1 ; c.rearrange(1)
  58         m.setConnectivity(c,cI,True)
  59         m.checkCoherency2()
  60         return m
  61
  62     def __traduceMeshForPolyhed(self,name,coords,arr0,arr1,arr2):
  63         nbCoords=len(coords)
  64         coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
  65         coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
  66         m=MEDCouplingUMesh(name,3)
  67         m.setCoords(coo)
  68         #
  69         arr2=arr2[:]-1
  70         arr0mc0=DataArrayInt(arr0) ; arr0mc0.computeOffsets2()
  71         arr0mc1=DataArrayInt(arr0).deepCpy()
  72         arr0mc2=DataArrayInt(len(arr0),2) ; arr0mc2[:,0]=DataArrayInt(arr0)-1 ; arr0mc2[:,1]=1 ; arr0mc2.rearrange(1) ; arr0mc2.computeOffsets2()
  73         arr0mc3=DataArrayInt.Range(0,2*len(arr0),2).buildExplicitArrByRanges(arr0mc2)
  74         arr1mc0=DataArrayInt(arr1) ; arr1mc0.computeOffsets2()
  75         arr1mc1=arr1mc0[arr0mc0] ; arr1mc1[1:]+=arr0mc0[1:]
  76         arr1mc2=DataArrayInt(arr1).deepCpy() ; arr1mc2+=1 ; arr1mc2.computeOffsets2()
  77         arr2mc0=(arr1mc2[1:])[arr0mc3]
  78         #
  79         c=DataArrayInt(arr1.size+arr2.size)
  80         c[arr1mc1[:-1]]=NORM_POLYHED
  81         c[arr2mc0]=-1
  82         a=arr2mc0.buildUnion(arr1mc1[:-1]).buildComplement(len(c))
  83         c[a]=DataArrayInt(arr2)
  84         #
  85         m.setConnectivity(c,arr1mc1,True)
  86         m.checkCoherency2()
  87         return m
  88
  89     def __traduceMeshForPolygon(self,name,coords,arr0,arr1):
  90         nbCoords=len(coords)
  91         coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
  92         coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
  93         m=MEDCouplingUMesh(name,2)
  94         m.setCoords(coo)
  95         #
  96         arr0_0=DataArrayInt(arr0+1) ; arr0_0.computeOffsets2()
  97         arr0_1=DataArrayInt(len(arr0),2) ; arr0_1[:,1]=DataArrayInt(arr0) ; arr0_1[:,0]=1 ; arr0_1.rearrange(1) ; arr0_1.computeOffsets2()
  98         arr0_2=DataArrayInt.Range(1,2*len(arr0),2).buildExplicitArrByRanges(arr0_1)
  99         c=DataArrayInt(len(arr0)+len(arr1)) ; c[:]=0 ; c[arr0_0[:-1]]=NORM_POLYGON
 100         c[arr0_2]=DataArrayInt(arr1-1)
 101         #
 102         m.setConnectivity(c,arr0_0,True)
 103         m.checkCoherency2()
 104         return m
 105
 106     def __convertGeo2MED(self,geoFileName):
 107         """ Convert all the geometry (all the meshes) contained in teh CASE file into MEDCouplingUMesh'es. """
 108         fd=open(os.path.join(self._dirName,geoFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0) ; fd.readline() ; fd.readline()
 109         name=fd.readline().strip() ; fd.readline() ; fd.readline()
 110         pos=fd.tell()
 111         mcmeshes=[]
 112         elt=fd.read(80) ; elt=elt.strip() ; pos+=80
 113         while pos!=end:
 114             if elt!="part":
 115                 raise Exception("Error on reading mesh #1 !")
 116             fd.seek(fd.tell()+4)
 117             meshName=fd.read(80).strip()
 118             if fd.read(len("coordinates"))!="coordinates":
 119                 raise Exception("Error on reading mesh #2 !")
 120             pos=fd.tell()
 121             typeOfCoo=np.memmap(fd,dtype='byte',mode='r',offset=int(pos),shape=(1)).tolist()[0]
 122             pos+=1+17*4
 123             nbNodes=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
 124             pos+=4
 125             coo=np.memmap(fd,dtype='float32',mode='r',offset=int(pos),shape=(nbNodes,3))
 126             pos+=nbNodes*3*4 ; fd.seek(pos)#np.array(0,dtype='float%i'%(typeOfCoo)).nbytes
 127             typ=fd.read(80).strip() ; pos=fd.tell()
 128             mcmeshes2=[]
 129             while pos!=end and typ!="part":
 130                 mctyp=self.dictMCTyp2[typ]
 131                 nbCellsOfType=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
 132                 pos+=4
 133                 if mctyp!=NORM_POLYHED and mctyp!=NORM_POLYGON:
 134                     nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(mctyp)
 135                     cells=np.memmap(fd,dtype='int32',mode='r',offset=pos,shape=(nbCellsOfType,nbNodesPerCell))
 136                     pos+=nbCellsOfType*nbNodesPerCell*4
 137                     fd.seek(pos)
 138                     mcmeshes2.append(self.__traduceMesh(meshName,typ,coo,cells))
 139                 elif mctyp==NORM_POLYHED:
 140                     nbOfFacesPerCell=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(nbCellsOfType,))
 141                     pos+=nbCellsOfType*4
 142                     szOfNbOfNodesPerFacePerCellArr=int(nbOfFacesPerCell.sum())
 143                     arr1=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNbOfNodesPerFacePerCellArr,))#arr1 -> nbOfNodesPerFacePerCellArr
 144                     pos+=szOfNbOfNodesPerFacePerCellArr*4
 145                     szOfNodesPerFacePerCellArr=arr1.sum()
 146                     arr2=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNodesPerFacePerCellArr,))#arr2 -> nodesPerFacePerCellArr
 147                     pos+=szOfNodesPerFacePerCellArr*4 ; fd.seek(pos)
 148                     mcmeshes2.append(self.__traduceMeshForPolyhed(meshName,coo,nbOfFacesPerCell,arr1,arr2))
 149                     pass
 150                 else:
 151                     nbOfNodesPerCell=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(nbCellsOfType,))
 152                     pos+=nbCellsOfType*4
 153                     szOfNbOfNodesPerCellArr=int(nbOfNodesPerCell.sum())
 154                     arr1=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNbOfNodesPerCellArr,))
 155                     pos+=szOfNbOfNodesPerCellArr*4  ; fd.seek(pos)
 156                     mcmeshes2.append(self.__traduceMeshForPolygon(meshName,coo,nbOfNodesPerCell,arr1))
 157                 if pos!=end:
 158                     elt=fd.read(80) ; elt=elt.strip() ; typ=elt[:] ; pos+=80
 159                     pass
 160                 pass
 161             coo=mcmeshes2[0].getCoords() ; name=mcmeshes2[0].getName()
 162             for itmesh in mcmeshes2: itmesh.setCoords(coo)
 163             m=MEDCouplingUMesh.MergeUMeshesOnSameCoords(mcmeshes2) ; m.setName(name)
 164             mcmeshes.append(m)
 165             pass
 166         #
 167         ms=MEDFileMeshes()
 168         ms.resize(len(mcmeshes))
 169         for i,m in enumerate(mcmeshes):
 170             mlm=MEDFileUMesh()
 171             mlm.setMeshAtLevel(0,m)
 172             ms.setMeshAtPos(i,mlm)
 173             pass
 174         return mcmeshes,ms
 175
 176     def __convertField(self,mlfields, mcmeshes, fileName, fieldName, discr, nbCompo, locId, it):
 177         """ Convert the fields. """
 178         stars=re.search("[\*]+",fileName).group()
 179         st="%0"+str(len(stars))+"i"
 180         trueFileName=fileName.replace(stars,st%(it))
 181         fd=open(os.path.join(self._dirName,trueFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0)
 182         name=fd.readline().strip().split(" ")[0]
 183         if name!=fieldName:
 184             raise Exception("ConvertField : mismatch")
 185         pos=fd.tell()
 186         st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
 187         while pos!=end:
 188             if st!="part":
 189                 raise Exception("ConvertField : mismatch #2")
 190             fdisc=MEDCouplingFieldDiscretization.New(self.discSpatial2[discr])
 191             meshId=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1)).tolist()[0]-1
 192             nbOfValues=fdisc.getNumberOfTuples(mcmeshes[meshId])
 193             vals2=DataArrayDouble(nbOfValues,nbCompo)
 194             fd.seek(pos+4)
 195             st=fd.read(80).strip() ; pos=fd.tell()
 196             offset=0
 197             while pos!=end and st!="part":
 198                 if st!="coordinates":
 199                     nbOfValsOfTyp=mcmeshes[meshId].getNumberOfCellsWithType(self.dictMCTyp2[st])
 200                 else:
 201                     nbOfValsOfTyp=nbOfValues
 202                     pass
 203                 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))
 204                 vals2[offset:offset+nbOfValsOfTyp]=DataArrayDouble(np.array(vals,dtype='float64')).fromNoInterlace()
 205                 pos+=nbOfValsOfTyp*nbCompo*4 ; fd.seek(pos)
 206                 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
 207                 offset+=nbOfValsOfTyp
 208                 pass
 209             f=MEDCouplingFieldDouble(self.discSpatial2[discr],ONE_TIME) ; f.setName("%s_%s"%(fieldName,mcmeshes[meshId].getName()))
 210             f.setMesh(mcmeshes[meshId]) ; f.setArray(vals2) ; f.setTime(float(it),it,-1)
 211             f.checkCoherency()
 212             mlfields[locId+meshId].appendFieldNoProfileSBT(f)
 213             pass
 214         pass
 215
 216     def loadInMEDFileDS(self):
 217         """ Load a CASE file into a MEDFileData object. """
 218         f=file(self._fileName)
 219         lines=f.readlines()
 220         ind=lines.index("GEOMETRY\n")
 221         if ind==-1:
 222             raise Exception("Error with file %s"%(fname))
 223         geoName=re.match("model:([\W]*)([\w\.]+)",lines[ind+1]).group(2)
 224         m1,m2=self.__convertGeo2MED(geoName)
 225         ind=lines.index("VARIABLE\n")
 226         fieldsInfo=[]
 227         for i in xrange(ind+1,lines.index("TIME\n")):
 228             m=re.match("^([\w]+)[\s]+\per[\s]+([\w]+)[\s]*\:[\s]*([\w]+)[\s]+([\S]+)$",lines[i])
 229             if m:
 230                 spatialDisc=m.groups()[1] ; fieldName=m.groups()[2] ; nbOfCompo=self.dictCompo2[m.groups()[0]] ; fieldFileName=m.groups()[3]
 231                 fieldsInfo.append((fieldName,spatialDisc,nbOfCompo,fieldFileName))
 232                 pass
 233             pass
 234
 235         expr=re.compile("number[\s]+of[\s]+steps[\s]*\:[\s]*([\d]+)")
 236         nbOfTimeSteps=int(expr.search(filter(expr.search,lines)[0]).group(1))
 237
 238         expr=re.compile("filename[\s]+start[\s]+number[\s]*\:[\s]*([\d]+)")
 239         startIt=int(expr.search(filter(expr.search,lines)[0]).group(1))
 240
 241         expr=re.compile("filename[\s]+increment[\s]*\:[\s]*([\d]+)")
 242         incrIt=int(expr.search(filter(expr.search,lines)[0]).group(1))
 243
 244         curIt=startIt
 245         mlfields=MEDFileFields()
 246         mlfields.resize(len(fieldsInfo)*len(m1))
 247         i=0
 248         for field in fieldsInfo:
 249             for m in m1:
 250                 mlfields.setFieldAtPos(i,MEDFileFieldMultiTS())
 251                 i+=1
 252                 pass
 253             pass
 254         for ts in xrange(nbOfTimeSteps):
 255             i=0
 256             for field in fieldsInfo:
 257                 self.__convertField(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt)
 258                 i+=len(m1)
 259                 pass
 260             curIt+=incrIt
 261             pass
 262         ret=MEDFileData()
 263         ret.setMeshes(m2)
 264         del mlfields[filter(lambda x: len(mlfields[x])==0,range(len(mlfields)))]
 265         ret.setFields(mlfields)
 266         return ret
 267
 268     pass