Salome HOME
d1dde8d4257d8bd3cc80fe51e83cb09f8962e46d
[tools/medcoupling.git] / src / MEDLoader / Swig / CaseReader.py
1 #  -*- coding: iso-8859-1 -*-
2 # Copyright (C) 2007-2016  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.checkConsistency()
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.computeOffsetsFull()
76         arr0mc1=DataArrayInt(arr0).deepCopy()
77         arr0mc2=DataArrayInt(len(arr0),2) ; arr0mc2[:,0]=DataArrayInt(arr0)-1 ; arr0mc2[:,1]=1 ; arr0mc2.rearrange(1) ; arr0mc2.computeOffsetsFull()
78         arr0mc3=DataArrayInt.Range(0,2*len(arr0),2).buildExplicitArrByRanges(arr0mc2)
79         arr1mc0=DataArrayInt(arr1) ; arr1mc0.computeOffsetsFull()
80         arr1mc1=arr1mc0[arr0mc0] ; arr1mc1[1:]+=arr0mc0[1:] 
81         arr1mc2=DataArrayInt(arr1).deepCopy() ; arr1mc2+=1 ; arr1mc2.computeOffsetsFull()
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.checkConsistency()
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.computeOffsetsFull()
102         arr0_1=DataArrayInt(len(arr0),2) ; arr0_1[:,1]=DataArrayInt(arr0) ; arr0_1[:,0]=1 ; arr0_1.rearrange(1) ; arr0_1.computeOffsetsFull()
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.checkConsistency()
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                 if typ[0]=='\0': pos+=1; continue
219                 mctyp=self.dictMCTyp2[typ]
220                 nbCellsOfType=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
221                 pos+=4
222                 if mctyp!=NORM_POLYHED and mctyp!=NORM_POLYGON:
223                     nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(mctyp)
224                     cells=np.memmap(fd,dtype='int32',mode='r',offset=pos,shape=(nbCellsOfType,nbNodesPerCell))
225                     pos+=nbCellsOfType*nbNodesPerCell*4
226                     fd.seek(pos)
227                     mcmeshes2.append(self.__traduceMesh(meshName,typ,coo,cells))
228                 elif mctyp==NORM_POLYHED:
229                     nbOfFacesPerCell=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(nbCellsOfType,))
230                     pos+=nbCellsOfType*4
231                     szOfNbOfNodesPerFacePerCellArr=int(nbOfFacesPerCell.sum())
232                     arr1=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNbOfNodesPerFacePerCellArr,))#arr1 -> nbOfNodesPerFacePerCellArr
233                     pos+=szOfNbOfNodesPerFacePerCellArr*4
234                     szOfNodesPerFacePerCellArr=arr1.sum()
235                     arr2=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNodesPerFacePerCellArr,))#arr2 -> nodesPerFacePerCellArr
236                     pos+=szOfNodesPerFacePerCellArr*4 ; fd.seek(pos)
237                     mcmeshes2.append(self.__traduceMeshForPolyhed(meshName,coo,nbOfFacesPerCell,arr1,arr2))
238                     pass
239                 else:
240                     nbOfNodesPerCell=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(nbCellsOfType,))
241                     pos+=nbCellsOfType*4
242                     szOfNbOfNodesPerCellArr=int(nbOfNodesPerCell.sum())
243                     arr1=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(szOfNbOfNodesPerCellArr,))
244                     pos+=szOfNbOfNodesPerCellArr*4  ; fd.seek(pos)
245                     mcmeshes2.append(self.__traduceMeshForPolygon(meshName,coo,nbOfNodesPerCell,arr1))
246                 if pos!=end:
247                     elt=fd.read(80) ; elt=elt.strip() ; typ=elt[:] ; pos+=80
248                     pass
249                 pass
250             if mcmeshes2:
251                 coo=mcmeshes2[0].getCoords() ; name=mcmeshes2[0].getName()
252                 for itmesh in mcmeshes2: itmesh.setCoords(coo)
253                 m=MEDCouplingUMesh.MergeUMeshesOnSameCoords(mcmeshes2) ; m.setName(name)
254                 mcmeshes.append(m)
255             pass
256         return mcmeshes
257         
258     
259     def __convertField(self,mlfields, mcmeshes, fileName, fieldName, discr, nbCompo, locId, it):
260         """ Convert the fields. """
261         stars=re.search("[\*]+",fileName).group()
262         st="%0"+str(len(stars))+"i"
263         trueFileName=fileName.replace(stars,st%(it))
264         fd=open(os.path.join(self._dirName,trueFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0)
265         name=fd.readline().strip().split(" ")[0]
266         if name!=fieldName:
267             raise Exception("ConvertField : mismatch")
268         pos=fd.tell()
269         st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
270         while pos!=end:
271             if st!="part":
272                 raise Exception("ConvertField : mismatch #2")
273             fdisc=MEDCouplingFieldDiscretization.New(self.discSpatial2[discr])
274             meshId=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1)).tolist()[0]-1
275             nbOfValues=fdisc.getNumberOfTuples(mcmeshes[meshId])
276             vals2=DataArrayDouble(nbOfValues,nbCompo)
277             fd.seek(pos+4)
278             st=fd.read(80).strip() ; pos=fd.tell()
279             offset=0
280             while pos!=end and st!="part":
281                 if st!="coordinates":
282                     nbOfValsOfTyp=mcmeshes[meshId].getNumberOfCellsWithType(self.dictMCTyp2[st])
283                 else:
284                     nbOfValsOfTyp=nbOfValues
285                     pass
286                 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))
287                 vals2[offset:offset+nbOfValsOfTyp]=DataArrayDouble(np.array(vals,dtype='float64')).fromNoInterlace()
288                 pos+=nbOfValsOfTyp*nbCompo*4 ; fd.seek(pos)
289                 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
290                 offset+=nbOfValsOfTyp
291                 pass
292             f=MEDCouplingFieldDouble(self.discSpatial2[discr],ONE_TIME) ; f.setName("%s_%s"%(fieldName,mcmeshes[meshId].getName()))
293             f.setMesh(mcmeshes[meshId]) ; f.setArray(vals2) ; f.setTime(float(it),it,-1)
294             f.checkConsistencyLight()
295             mlfields[locId+meshId].appendFieldNoProfileSBT(f)
296             pass
297
298     def __convertFieldFortran(self,mlfields, mcmeshes, fileName, fieldName, discr, nbCompo, locId, it):
299         """ Convert the fields. """
300         if re.search("[\*]+",fileName):
301             stars=re.search("[\*]+",fileName).group()
302             st="%0"+str(len(stars))+"i"
303             trueFileName=fileName.replace(stars,st%(it))
304             pass
305         else:
306             trueFileName=fileName
307             pass
308         fd=open(os.path.join(self._dirName,trueFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0)
309         name=fd.read(80)
310         if fieldName not in name:
311             raise Exception("ConvertField : mismatch")
312         pos=fd.tell()
313         st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
314         if "part" not in st:
315             raise Exception("ConvertField : mismatch #2")
316         st=fd.read(80).strip() ; pos=fd.tell()
317         pos+=12 # I love it
318         offset=0
319         nbTurn=0
320         while pos!=end and "part" not in st:
321             fdisc=MEDCouplingFieldDiscretization.New(self.discSpatial2[discr])
322             nbOfValues=fdisc.getNumberOfTuples(mcmeshes[nbTurn])
323             vals2=DataArrayDouble(nbOfValues,nbCompo)
324             pos+=24 # I love it again !
325             nbOfValsOfTyp=np.memmap(fd,dtype='>i4',mode='r',offset=pos,shape=(1)).tolist()[0]/4
326             pos+=4
327             vals=np.zeros(dtype=">f4",shape=(nbOfValsOfTyp*nbCompo))
328             for iii in xrange(nbCompo):
329                 valsTmp=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos),shape=(nbOfValsOfTyp))
330                 vals[iii*nbOfValsOfTyp:(iii+1)*nbOfValsOfTyp]=valsTmp
331                 pos+=nbOfValsOfTyp*4
332                 pos+=2*4 ## hey hey, that is the ultimate class !
333                 vals2.setInfoOnComponent(iii,chr(ord('X')+iii))
334                 pass
335             if pos>end:
336                 pos=end
337                 pass
338             vals=vals.reshape(nbOfValsOfTyp,nbCompo)
339             vals2[offset:offset+nbOfValsOfTyp]=DataArrayDouble(np.array(vals,dtype='float64')).fromNoInterlace()
340             if pos!=end:
341                 fd.seek(pos)
342                 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
343                 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
344                 pass
345             f=MEDCouplingFieldDouble(self.discSpatial2[discr],ONE_TIME) ; f.setName("%s_%s"%(fieldName,mcmeshes[nbTurn].getName()))
346             f.setMesh(mcmeshes[nbTurn]) ; f.setArray(vals2) ; f.setTime(float(it),it,-1)
347             f.checkConsistencyLight()
348             mlfields[locId+nbTurn].appendFieldNoProfileSBT(f)
349             nbTurn+=1
350             pass
351         pass
352     
353     def loadInMEDFileDS(self):
354         """ Load a CASE file into a MEDFileData object. """
355         f=file(self._fileName)
356         lines=f.readlines()
357         ind=lines.index("GEOMETRY\n")
358         if ind==-1:
359             raise Exception("Error with file %s"%(fname))
360         geoName=re.match("model:([\W]*)([\w\.]+)",lines[ind+1]).group(2)
361         m1,m2,typeOfFile=self.__convertGeo2MED(geoName)
362         fieldsInfo=[] ; nbOfTimeSteps=0
363         if "VARIABLE\n" in lines:
364             ind=lines.index("VARIABLE\n")
365             end=len(lines)-1
366             if "TIME\n" in lines:
367                 end=lines.index("TIME\n")
368                 pass
369             for i in xrange(ind+1,end):
370                 m=re.match("^([\w]+)[\s]+\per[\s]+([\w]+)[\s]*\:[\s]*([\w]+)[\s]+([\S]+)$",lines[i])
371                 if m:
372                     if m.groups()[0]=="constant":
373                         continue
374                     spatialDisc=m.groups()[1] ; fieldName=m.groups()[2] ; nbOfCompo=self.dictCompo2[m.groups()[0]] ; fieldFileName=m.groups()[3]
375                     fieldsInfo.append((fieldName,spatialDisc,nbOfCompo,fieldFileName))
376                     pass
377                 pass
378             
379             expr=re.compile("number[\s]+of[\s]+steps[\s]*\:[\s]*([\d]+)")
380             tmp=filter(expr.search,lines)
381             if len(tmp)!=0:
382                 nbOfTimeSteps=int(expr.search(filter(expr.search,lines)[0]).group(1))
383                 expr=re.compile("filename[\s]+start[\s]+number[\s]*\:[\s]*([\d]+)")
384                 startIt=int(expr.search(filter(expr.search,lines)[0]).group(1))
385                 expr=re.compile("filename[\s]+increment[\s]*\:[\s]*([\d]+)")
386                 incrIt=int(expr.search(filter(expr.search,lines)[0]).group(1))
387             else:
388                 nbOfTimeSteps=1
389                 startIt=0
390                 incrIt=1
391                 pass
392             curIt=startIt
393             pass
394         mlfields=MEDFileFields()
395         mlfields.resize(len(fieldsInfo)*len(m1))
396         i=0
397         for field in fieldsInfo:
398             for m in m1:
399                 mlfields.setFieldAtPos(i,MEDFileFieldMultiTS())
400                 i+=1
401                 pass
402             pass
403         for ts in xrange(nbOfTimeSteps):
404             i=0
405             for field in fieldsInfo:
406                 if typeOfFile:
407                     self.__convertField(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt);
408                 else:
409                     self.__convertFieldFortran(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt)
410                     pass
411                 i+=len(m1)
412                 pass
413             curIt+=incrIt
414             pass
415         ret=MEDFileData()
416         ret.setMeshes(m2)
417         del mlfields[filter(lambda x: len(mlfields[x])==0,range(len(mlfields)))]
418         ret.setFields(mlfields)
419         return ret
420
421     pass