From bebd0c9d73e0b52bb0d4c4b51c0bcb8fd79e5c6a Mon Sep 17 00:00:00 2001 From: ageay Date: Fri, 31 Jan 2014 15:54:28 +0000 Subject: [PATCH] "Supporting" Fortran case file. --- src/MEDLoader/Swig/CaseReader.py | 197 +++++++++++++++++++++++++++---- 1 file changed, 173 insertions(+), 24 deletions(-) diff --git a/src/MEDLoader/Swig/CaseReader.py b/src/MEDLoader/Swig/CaseReader.py index 87d2fbcce..7dc270e57 100644 --- a/src/MEDLoader/Swig/CaseReader.py +++ b/src/MEDLoader/Swig/CaseReader.py @@ -19,6 +19,7 @@ # # Author Anthony GEAY (CEA/DEN/DM2S/STMF/LGLS) +# http://www-vis.lbl.gov/NERSC/Software/ensight/doc/OnlineHelp/UM-C11.pdf import numpy as np from MEDLoader import * from CaseIO import CaseIO @@ -53,7 +54,11 @@ class CaseReader(CaseIO): cI=DataArrayInt(len(cells)+1) ; cI.iota() ; cI*=nbNodesPerCell+1 # cells2=cells.reshape(len(cells),nbNodesPerCell) - c2=DataArrayInt(cells2) + if cells2.dtype=='int32': + c2=DataArrayInt(cells2) + else: + c2=DataArrayInt(np.array(cells2,dtype="int32")) + pass c=DataArrayInt(len(cells),nbNodesPerCell+1) ; c[:,0]=ct ; c[:,1:]=c2-1 ; c.rearrange(1) m.setConnectivity(c,cI,True) m.checkCoherency2() @@ -105,13 +110,96 @@ class CaseReader(CaseIO): def __convertGeo2MED(self,geoFileName): """ Convert all the geometry (all the meshes) contained in teh CASE file into MEDCouplingUMesh'es. """ - fd=open(os.path.join(self._dirName,geoFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0) ; fd.readline() ; fd.readline() + fd=open(os.path.join(self._dirName,geoFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0) ; title=fd.read(80) + title=title.strip().lower() + if "binary" not in title: + raise Exception("Error only binary geo files are supported for the moment !") + pass + zeType=True + if "fortran" in title: + mcmeshes=self.__convertGeo2MEDFortran(fd,end) ; zeType=False + else: + mcmeshes=self.__convertGeo2MEDC(fd,end) + # + ms=MEDFileMeshes() + ms.resize(len(mcmeshes)) + for i,m in enumerate(mcmeshes): + mlm=MEDFileUMesh() + mlm.setMeshAtLevel(0,m) + ms.setMeshAtPos(i,mlm) + pass + return mcmeshes,ms,zeType + + def __convertGeo2MEDFortran(self,fd,end): + mcmeshes=[] + fd.read(80) # comment 1 + fd.read(80) # comment 2 + fd.read(80) # node id + fd.read(80) # element id + pos=fd.tell() + elt=fd.read(80) ; elt=elt.strip() ; pos=fd.tell() + mcmeshes2=[] + typ="part" + nbOfTurn=0 + while abs(pos-end)>8 and "part" in typ: + if "part" not in elt: + raise Exception("Error on reading mesh fortran #1 !") + fd.seek(fd.tell()+4)# skip # + tmp=fd.read(80) ; meshName=tmp.split("P")[-1] + tmp=fd.read(80) + if "coordinates" not in tmp: + raise Exception("Error on reading mesh fortran #2 !") + pos=fd.tell() # 644 + if nbOfTurn==0: + pos+=76 # what else ? + else: + pos+=40 + pass + nbNodes=np.memmap(fd,dtype='>i4',mode='r',offset=int(pos),shape=(1,)).tolist()[0] + pos+=12 # what else ? + a=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos),shape=(nbNodes)) + b=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos+nbNodes*4+2*4),shape=(nbNodes)) + c=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos+nbNodes*2*4+4*4),shape=(nbNodes)) + coo=np.zeros(dtype=">f4",shape=(nbNodes*3)) + coo[:nbNodes]=a ; coo[nbNodes:2*nbNodes]=b ; coo[2*nbNodes:]=c + coo=coo.reshape(nbNodes,3) + pos+=nbNodes*3*4 ; fd.seek(pos)#np.array(0,dtype='float%i'%(typeOfCoo)).nbytes + typ=fd.read(80).strip() ; pos=fd.tell() + zeK="" + for k in self.dictMCTyp2.keys(): + if k in typ: + zeK=k + break + pass + pass + pos+=8*4 # yeh man ! + nbCellsOfType=np.memmap(fd,dtype='>i4',mode='r',offset=int(pos),shape=(1,)).tolist()[0] + pos+=4 # for the number of cells + pos+=2*4 # because it's great ! + nbNodesPerCell=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(self.dictMCTyp2[zeK]) + nodalConn=np.memmap(fd,dtype='>i4',mode='r',offset=pos,shape=(nbCellsOfType,nbNodesPerCell)) + meshName=meshName.strip() + mcmeshes2.append(self.__traduceMesh(meshName,zeK,coo,nodalConn)) + pos+=nbNodesPerCell*nbCellsOfType*4 + if abs(pos-end)>8: + fd.seek(pos) ;elt=fd.read(80) ; typ=elt[:] ; pos+=80 + pass + nbOfTurn+=1 + pass + #coo=mcmeshes2[0].getCoords() ; name=mcmeshes2[0].getName() + #for itmesh in mcmeshes2: itmesh.setCoords(coo) + #m=MEDCouplingUMesh.MergeUMeshesOnSameCoords(mcmeshes2) ; m.setName(name) + #mcmeshes.append(m) + return mcmeshes2 + + def __convertGeo2MEDC(self,fd,end): + fd.readline() name=fd.readline().strip() ; fd.readline() ; fd.readline() pos=fd.tell() mcmeshes=[] elt=fd.read(80) ; elt=elt.strip() ; pos+=80 while pos!=end: - if elt!="part": + if "part" not in elt: raise Exception("Error on reading mesh #1 !") fd.seek(fd.tell()+4) meshName=fd.read(80).strip() @@ -163,15 +251,8 @@ class CaseReader(CaseIO): m=MEDCouplingUMesh.MergeUMeshesOnSameCoords(mcmeshes2) ; m.setName(name) mcmeshes.append(m) pass - # - ms=MEDFileMeshes() - ms.resize(len(mcmeshes)) - for i,m in enumerate(mcmeshes): - mlm=MEDFileUMesh() - mlm.setMeshAtLevel(0,m) - ms.setMeshAtPos(i,mlm) - pass - return mcmeshes,ms + return mcmeshes + def __convertField(self,mlfields, mcmeshes, fileName, fieldName, discr, nbCompo, locId, it): """ Convert the fields. """ @@ -211,6 +292,60 @@ class CaseReader(CaseIO): f.checkCoherency() mlfields[locId+meshId].appendFieldNoProfileSBT(f) pass + + def __convertFieldFortran(self,mlfields, mcmeshes, fileName, fieldName, discr, nbCompo, locId, it): + """ Convert the fields. """ + if re.search("[\*]+",fileName): + stars=re.search("[\*]+",fileName).group() + st="%0"+str(len(stars))+"i" + trueFileName=fileName.replace(stars,st%(it)) + pass + else: + trueFileName=fileName + pass + fd=open(os.path.join(self._dirName,trueFileName),"r+b") ; fd.seek(0,2) ; end=fd.tell() ; fd.seek(0) + name=fd.read(80) + if fieldName not in name: + raise Exception("ConvertField : mismatch") + pos=fd.tell() + st=fd.read(80) ; st=st.strip() ; pos=fd.tell() + if "part" not in st: + raise Exception("ConvertField : mismatch #2") + st=fd.read(80).strip() ; pos=fd.tell() + pos+=12 # I love it + offset=0 + nbTurn=0 + while pos!=end and "part" not in st: + fdisc=MEDCouplingFieldDiscretization.New(self.discSpatial2[discr]) + nbOfValues=fdisc.getNumberOfTuples(mcmeshes[nbTurn]) + vals2=DataArrayDouble(nbOfValues,nbCompo) + pos+=24 # I love it again ! + nbOfValsOfTyp=np.memmap(fd,dtype='>i4',mode='r',offset=pos,shape=(1)).tolist()[0]/4 + pos+=4 + vals=np.zeros(dtype=">f4",shape=(nbOfValsOfTyp*nbCompo)) + for iii in xrange(nbCompo): + valsTmp=np.memmap(fd,dtype='>f4',mode='r',offset=int(pos),shape=(nbOfValsOfTyp)) + vals[iii*nbOfValsOfTyp:(iii+1)*nbOfValsOfTyp]=valsTmp + pos+=nbOfValsOfTyp*4 + pos+=2*4 ## hey hey, that is the ultimate class ! + vals2.setInfoOnComponent(iii,chr(ord('X')+iii)) + pass + if pos>end: + pos=end + pass + vals=vals.reshape(nbOfValsOfTyp,nbCompo) + vals2[offset:offset+nbOfValsOfTyp]=DataArrayDouble(np.array(vals,dtype='float64')).fromNoInterlace() + if pos!=end: + fd.seek(pos) + st=fd.read(80) ; st=st.strip() ; pos=fd.tell() + st=fd.read(80) ; st=st.strip() ; pos=fd.tell() + pass + f=MEDCouplingFieldDouble(self.discSpatial2[discr],ONE_TIME) ; f.setName("%s_%s"%(fieldName,mcmeshes[nbTurn].getName())) + f.setMesh(mcmeshes[nbTurn]) ; f.setArray(vals2) ; f.setTime(float(it),it,-1) + f.checkCoherency() + mlfields[locId+nbTurn].appendFieldNoProfileSBT(f) + nbTurn+=1 + pass pass def loadInMEDFileDS(self): @@ -221,26 +356,36 @@ class CaseReader(CaseIO): if ind==-1: raise Exception("Error with file %s"%(fname)) geoName=re.match("model:([\W]*)([\w\.]+)",lines[ind+1]).group(2) - m1,m2=self.__convertGeo2MED(geoName) - ind=lines.index("VARIABLE\n") + m1,m2,typeOfFile=self.__convertGeo2MED(geoName) fieldsInfo=[] - for i in xrange(ind+1,lines.index("TIME\n")): + ind=lines.index("VARIABLE\n") + end=len(lines)-1 + if "TIME\n" in lines: + end=lines.index("TIME\n") + pass + for i in xrange(ind+1,end): m=re.match("^([\w]+)[\s]+\per[\s]+([\w]+)[\s]*\:[\s]*([\w]+)[\s]+([\S]+)$",lines[i]) if m: + if m.groups()[0]=="constant": + continue spatialDisc=m.groups()[1] ; fieldName=m.groups()[2] ; nbOfCompo=self.dictCompo2[m.groups()[0]] ; fieldFileName=m.groups()[3] fieldsInfo.append((fieldName,spatialDisc,nbOfCompo,fieldFileName)) pass pass expr=re.compile("number[\s]+of[\s]+steps[\s]*\:[\s]*([\d]+)") - nbOfTimeSteps=int(expr.search(filter(expr.search,lines)[0]).group(1)) - - expr=re.compile("filename[\s]+start[\s]+number[\s]*\:[\s]*([\d]+)") - startIt=int(expr.search(filter(expr.search,lines)[0]).group(1)) - - expr=re.compile("filename[\s]+increment[\s]*\:[\s]*([\d]+)") - incrIt=int(expr.search(filter(expr.search,lines)[0]).group(1)) - + tmp=filter(expr.search,lines) + if len(tmp)!=0: + nbOfTimeSteps=int(expr.search(filter(expr.search,lines)[0]).group(1)) + expr=re.compile("filename[\s]+start[\s]+number[\s]*\:[\s]*([\d]+)") + startIt=int(expr.search(filter(expr.search,lines)[0]).group(1)) + expr=re.compile("filename[\s]+increment[\s]*\:[\s]*([\d]+)") + incrIt=int(expr.search(filter(expr.search,lines)[0]).group(1)) + else: + nbOfTimeSteps=1 + startIt=0 + incrIt=1 + pass curIt=startIt mlfields=MEDFileFields() mlfields.resize(len(fieldsInfo)*len(m1)) @@ -254,7 +399,11 @@ class CaseReader(CaseIO): for ts in xrange(nbOfTimeSteps): i=0 for field in fieldsInfo: - self.__convertField(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt) + if typeOfFile: + self.__convertField(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt); + else: + self.__convertFieldFortran(mlfields,m1,field[3],field[0],field[1],field[2],i,curIt) + pass i+=len(m1) pass curIt+=incrIt -- 2.39.2