#
# 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
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()
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()
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. """
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):
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))
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
