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