1 # -*- coding: iso-8859-1 -*-
2 # Copyright (C) 2007-2020 CEA/DEN, EDF R&D
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.
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.
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
18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 # Author : Anthony GEAY (CEA/DEN/DM2S/STMF/LGLS)
22 # http://www-vis.lbl.gov/NERSC/Software/ensight/doc/OnlineHelp/UM-C11.pdf
24 from MEDLoader import *
25 from CaseIO import CaseIO
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.
34 def New(cls,fileName):
35 """ Static constructor. """
36 return CaseReader(fileName)
39 def __init__(self,fileName):
42 self._fileName=fileName
43 self._dirName=os.path.dirname(self._fileName)
46 def __traduceMesh(self,name,typ,coords,cells):
47 """ Convert a CASE mesh into a MEDCouplingUMesh. """
48 name = name.decode("ascii")
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))
55 nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(ct)
56 cI=DataArrayInt(len(cells)+1) ; cI.iota() ; cI*=nbNodesPerCell+1
58 cells2=cells.reshape(len(cells),nbNodesPerCell)
59 if cells2.dtype=='int32':
60 c2=DataArrayInt(cells2)
62 c2=DataArrayInt(np.array(cells2,dtype="int32"))
64 c=DataArrayInt(len(cells),nbNodesPerCell+1) ; c[:,0]=ct ; c[:,1:]=c2-1 ; c.rearrange(1)
65 m.setConnectivity(c,cI,True)
69 def __traduceMeshForPolyhed(self,name,coords,arr0,arr1,arr2):
71 coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
72 coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
73 m=MEDCouplingUMesh(name,3)
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]
86 c=DataArrayInt(arr1.size+arr2.size)
87 c[arr1mc1[:-1]]=NORM_POLYHED
89 a=arr2mc0.buildUnion(arr1mc1[:-1]).buildComplement(len(c))
90 c[a]=DataArrayInt(arr2)
92 m.setConnectivity(c,arr1mc1,True)
96 def __traduceMeshForPolygon(self,name,coords,arr0,arr1):
98 coo=np.array(coords,dtype="float64") ; coo=coo.reshape(nbCoords,3)
99 coo=DataArrayDouble(coo) ; coo=coo.fromNoInterlace()
100 m=MEDCouplingUMesh(name,2)
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)
109 m.setConnectivity(c,arr0_0,True)
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 !")
121 if b"fortran" in title:
122 mcmeshes=self.__convertGeo2MEDFortran(fd,end) ; zeType=False
124 mcmeshes=self.__convertGeo2MEDC(fd,end)
127 ms.resize(len(mcmeshes))
128 for i,m in enumerate(mcmeshes):
130 mlm.setMeshAtLevel(0,m)
131 ms.setMeshAtPos(i,mlm)
133 return mcmeshes,ms,zeType
135 def __convertGeo2MEDFortran(self,fd,end):
137 fd.read(80) # comment 1
138 fd.read(80) # comment 2
139 fd.read(80) # node id
140 fd.read(80) # element id
142 elt=fd.read(80) ; elt=elt.strip() ; pos=fd.tell()
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]
152 if b"coordinates" not in tmp:
153 raise Exception("Error on reading mesh fortran #2 !")
156 pos+=76 # what else ?
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()
171 for k in self.dictMCTyp2:
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
187 fd.seek(pos) ;elt=fd.read(80) ; typ=elt[:] ; pos+=80
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)
197 def __convertGeo2MEDC(self,fd,end):
199 #name=fd.readline().strip() ; fd.readline() ; fd.readline()
201 descrip=fd.read(80).strip() ; fd.read(80) ; fd.read(80)
204 elt=fd.read(80) ; elt=elt.strip() ; pos+=80
206 if b"part" not in elt:
207 raise Exception("Error on reading mesh #1 !")
209 meshName=fd.read(80).strip()
210 if fd.read(len("coordinates"))!=b"coordinates":
211 raise Exception("Error on reading mesh #2 !")
213 typeOfCoo=np.memmap(fd,dtype='byte',mode='r',offset=int(pos),shape=(1)).tolist()[0]
215 nbNodes=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(1,)).tolist()[0]
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()
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]
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
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,))
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))
244 nbOfNodesPerCell=np.memmap(fd,dtype='int32',mode='r',offset=int(pos),shape=(nbCellsOfType,))
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))
251 elt=fd.read(80) ; elt=elt.strip() ; typ=elt[:] ; pos+=80
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)
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")
271 raise Exception("ConvertField : mismatch")
273 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
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 ):
281 nbOfValues=fdisc.getNumberOfTuples(mcmeshes[meshId])
282 vals2=DataArrayDouble(nbOfValues,nbCompo)
284 st=fd.read(80).strip() ; pos=fd.tell()
286 while pos!=end and st!=b"part":
287 if st!=b"coordinates":
288 nbOfValsOfTyp=mcmeshes[meshId].getNumberOfCellsWithType(self.dictMCTyp2[st])
290 nbOfValsOfTyp=nbOfValues
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
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)
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))
312 trueFileName=fileName
314 fd=open(os.path.join(self._dirName,trueFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0)
316 if fieldName not in name:
317 raise Exception("ConvertField : mismatch")
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()
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
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
338 pos+=2*4 ## hey hey, that is the ultimate class !
339 vals2.setInfoOnComponent(iii,chr(ord('X')+iii))
344 vals=vals.reshape(nbOfValsOfTyp,nbCompo)
345 vals2[offset:offset+nbOfValsOfTyp]=DataArrayDouble(np.array(vals,dtype='float64')).fromNoInterlace()
348 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
349 st=fd.read(80) ; st=st.strip() ; pos=fd.tell()
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)
359 def loadInMEDFileDS(self):
360 """ Load a CASE file into a MEDFileData object. """
361 f=open(self._fileName)
363 ind=lines.index("GEOMETRY\n")
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")
372 if "TIME\n" in lines:
373 end=lines.index("TIME\n")
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])
378 if m.groups()[0]=="constant":
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))
386 expr=re.compile("number[\s]+of[\s]+steps[\s]*\:[\s]*([\d]+)")
387 tmp = [line for line in lines if expr.search(line)]
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))
401 mlfields=MEDFileFields()
402 mlfields.resize(len(fieldsInfo)*len(m1))
404 for field in fieldsInfo:
406 mlfields.setFieldAtPos(i,MEDFileFieldMultiTS())
410 for ts in range(nbOfTimeSteps):
412 for field in fieldsInfo:
414 self.__convertField(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt);
416 self.__convertFieldFortran(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt)
424 del mlfields[[x for x in range(len(mlfields)) if len(mlfields[x]) == 0]]
425 ret.setFields(mlfields)