src/MEDLoader/Swig/CaseWriter.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 import numpy as np
  23 from CaseIO import CaseIO
  24 from MEDLoader import *
  25 import sys,re,os,mmap
  26
  27 ### www-vis.lbl.gov/NERSC/Software/ensight/doc/OnlineHelp/UM-C11.pdf
  28
  29 class CaseWriter(CaseIO):
  30     """ Converting MED file format in memory to a the Case file format (Ensight).
  31     A new file with the same base name and the .case extension is created with its dependencies (.geo ...).
  32     """
  33
  34     header="""FORMAT
  35 type: ensight gold
  36 GEOMETRY
  37 model: %(geofilewithoutpath)s
  38 """
  39     header_varpart="""VARIABLE"""
  40     header_timepart="""TIME
  41 time set:               1
  42 number of steps:        %(NbTimeSteps)i
  43 filename start number:  0
  44 filename increment:     1
  45 time values:
  46 %(TimeValues)s
  47 """
  48
  49     @classmethod
  50     def New(cls):
  51         """ Static constructor. """
  52         return CaseWriter()
  53         pass
  54
  55     def __init__(self):
  56         """ Constructor """
  57         self.__export_groups=False
  58         pass
  59
  60     def setMEDFileDS(self,medData):
  61         """ Input should be MEDFileData instance  """
  62         self._med_data=medData
  63         pass
  64
  65     def isExportingGroups(self):
  66         """ return the status of exporting groups policy """
  67         return self.__export_groups
  68
  69     def setExportingGroups(self,status):
  70         assert(isinstance(status,bool))
  71         self.__export_groups=status
  72         pass
  73
  74
  75     def write(self,fileName):
  76         """ Write into the specified fileName series the result """
  77         self._file_name=fileName
  78         self._base_name_without_dir=os.path.splitext(os.path.basename(self._file_name))[0]
  79         self._l=self._file_name.split(os.path.sep) ; self._l[-1]=os.path.splitext(self._l[-1])[0]
  80         self._base_name_with_dir=os.path.sep.join(self._l)
  81         self._real_written_file_name=[]
  82         self._dico={}
  83         for mesh in self._med_data.getMeshes():
  84             additionnalFileNamePart=""
  85             if len(self._med_data.getMeshes())!=1:
  86                 additionnalFileNamePart="_%s"%(mesh.getName())
  87                 pass
  88             self._dico["geofilewithoutpath"]="%s%s.geo"%(self._base_name_without_dir,additionnalFileNamePart)
  89             h0=self.header%self._dico
  90             self.__writeMeshesPart(mesh,"%s%s.geo"%(self._base_name_with_dir,additionnalFileNamePart))
  91             #
  92             h2=self.__writeFieldsPart(self._med_data.getFields().partOfThisLyingOnSpecifiedMeshName(mesh.getName()))
  93             realWrittenCaseFileNameForCurMesh="%s%s.case"%(self._base_name_with_dir,additionnalFileNamePart)
  94             fheader=open(realWrittenCaseFileNameForCurMesh,"w") ; fheader.write((h0+h2)%self._dico)
  95             self._real_written_file_name.append(realWrittenCaseFileNameForCurMesh)
  96             pass
  97         return self._real_written_file_name
  98
  99     def __writeMeshesPart(self,mdm,meshfn):
 100         try:
 101             os.remove(meshfn)
 102         except:
 103             pass
 104         f=open(meshfn,"w+b")
 105         sz=5*80
 106         #
 107         assert(isinstance(mdm,MEDFileUMesh))
 108         ms2=[[mdm.getMeshAtLevel(lev) for lev in mdm.getNonEmptyLevels()[:1]]]
 109         if self.__export_groups:
 110             for grpnm in mdm.getGroupsNames():
 111                 ms3=[]
 112                 for lev in mdm.getGrpNonEmptyLevels(grpnm)[:1]:
 113                     m=mdm.getGroup(lev,grpnm) ; m.zipCoords()
 114                     ms3.append(m)
 115                     pass
 116                 ms2.append(ms3)
 117                 pass
 118             pass
 119         for ms in ms2:
 120             nn=ms[0].getNumberOfNodes()
 121             sz+=self.__computeSizeOfGeoFile(ms,nn)
 122             pass
 123         pass
 124         a=np.memmap(f,dtype='byte',mode='w+',offset=0,shape=(sz,)) ; a.flush() # truncate to set the size of the file
 125         mm=mmap.mmap(f.fileno(),offset=0,length=0)
 126         mm.write(self.__str80("C Binary"))
 127         mm.write(self.__str80("Exported from MEDCoupling/MEDLoader version %s"%(MEDCouplingVersionStr())))
 128         mm.write(self.__str80("Conversion using CaseWriter class"))
 129         mm.write(self.__str80("node id off"))
 130         mm.write(self.__str80("element id off"))
 131         for iii,ms in enumerate(ms2):
 132             nn=ms[0].getNumberOfNodes()
 133             mm.write(self.__str80("part"))
 134             a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(1,))
 135             a[0]=iii+1 ; a.flush() ; mm.seek(mm.tell()+4) # part number maybe to change ?
 136             name=ms[0].getName()
 137             if iii>0:
 138                 name="%s_%s"%(ms2[0][0].getName(),name)
 139                 pass
 140             mm.write(self.__str80(name))
 141             mm.write(self.__str80("coordinates"))
 142             a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(1,))
 143             a[0]=nn ; a.flush() # number of nodes
 144             mm.seek(mm.tell()+4)
 145             coo=ms[0].getCoords()
 146             spaceDim=coo.getNumberOfComponents()
 147             if spaceDim!=3:
 148                 coo=coo.changeNbOfComponents(3,0.)
 149                 pass
 150             a=np.memmap(f,dtype='float32',mode='w+',offset=mm.tell(),shape=(3,nn))
 151             c=coo.toNoInterlace() ; c.rearrange(1) ; cnp=c.toNumPyArray() ; cnp=cnp.reshape(3,nn)
 152             a[:]=cnp ; a.flush() ; mm.seek(mm.tell()+3*nn*4)
 153             for m in ms:
 154                 i=0
 155                 for typ2,nbelem,dummy in m.getDistributionOfTypes():
 156                     typ=typ2
 157                     if typ not in self.dictMCTyp:
 158                         typ=MEDCouplingMesh.GetCorrespondingPolyType(typ)
 159                         pass
 160                     mp=m[i:i+nbelem]
 161                     mm.write(self.__str80(self.dictMCTyp_str[typ]))
 162                     a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(1,))
 163                     a[0]=nbelem ; a.flush() ; mm.seek(mm.tell()+4)
 164                     if typ!=NORM_POLYHED and typ!=NORM_POLYGON:
 165                         nbNodesPerElem=MEDCouplingMesh.GetNumberOfNodesOfGeometricType(typ)
 166                         c=mp.getNodalConnectivity() ; c.rearrange(nbNodesPerElem+1) ; c=c[:,1:] ; c.rearrange(1) ; c+=1
 167                         a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(nbNodesPerElem*nbelem,))
 168                         a[:]=c.toNumPyArray() ; a.flush() ; mm.seek(mm.tell()+nbNodesPerElem*nbelem*4)
 169                         pass
 170                     elif typ==NORM_POLYHED:
 171                         mp.orientCorrectlyPolyhedrons()
 172                         c=mp.computeNbOfFacesPerCell()
 173                         a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(nbelem,))
 174                         a[:]=c.toNumPyArray(); a.flush() ; mm.seek(mm.tell()+nbelem*4)
 175                         c=mp.getNodalConnectivity()[:] ; c.pushBackSilent(-1) ; c[mp.getNodalConnectivityIndex()[:-1]]=-1 ; ids=c.findIdsEqual(-1) ; nbOfNodesPerFace=ids.deltaShiftIndex()-1
 176                         a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(len(nbOfNodesPerFace),))
 177                         a[:]=nbOfNodesPerFace.toNumPyArray() ; a.flush() ; mm.seek(mm.tell()+len(nbOfNodesPerFace)*4)
 178                         ids2=ids.buildComplement(ids.back()+1)
 179                         c2=mp.getNodalConnectivity()[ids2]+1
 180                         a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(len(c2),))
 181                         a[:]=c2.toNumPyArray() ; a.flush() ; mm.seek(mm.tell()+len(c2)*4)
 182                         pass
 183                     else:
 184                         nbOfNodesPerCell=mp.getNodalConnectivityIndex().deltaShiftIndex()-1
 185                         a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(len(nbOfNodesPerCell),))
 186                         a[:]=nbOfNodesPerCell.toNumPyArray() ; a.flush() ; mm.seek(mm.tell()+len(nbOfNodesPerCell)*4)
 187                         ids2=mp.getNodalConnectivityIndex().buildComplement(mp.getNodalConnectivityIndex().back()+1)
 188                         c2=mp.getNodalConnectivity()[ids2]+1
 189                         a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(len(c2),))
 190                         a[:]=c2.toNumPyArray() ; a.flush() ; mm.seek(mm.tell()+len(c2)*4)
 191                         pass
 192                     i+=nbelem
 193                     pass
 194                 pass
 195             pass
 196         pass
 197
 198     def __writeFieldsPart(self,mdfs):
 199         if not mdfs:
 200             return ""
 201         self._ze_top_dict={}
 202         its,areForgottenTS=mdfs.getCommonIterations()
 203         if areForgottenTS:
 204             print(("WARNING : some iterations are NOT present in all fields ! Kept iterations are : %s !"%(str(its))))
 205             pass
 206         TimeValues=""
 207         for it in its:
 208             TimeValues+="%s\n"%(str(mdfs[0][it].getTime()[-1]))
 209             pass
 210         dictVars={}
 211         for mdf in mdfs:
 212             nbCompo=mdf.getNumberOfComponents()
 213             if nbCompo not in self.dictCompo:
 214                 l = [x for x in self.dictCompo if x - nbCompo > 0]
 215                 if len(l)==0:
 216                     print(("Field \"%s\" will be ignored because number of components (%i) is too big to be %s supported by case files !"%(mdf.getName(),nbCompo,str(list(self.dictCompo.keys())))))
 217                     continue
 218                     pass
 219                 print(("WARNING : Field \"%s\" will have its number of components (%i) set to %i, in order to be supported by case files (must be in %s) !"%(mdf.getName(),nbCompo,l[0],str(list(self.dictCompo.keys())))))
 220                 nbCompo=l[0]
 221                 pass
 222             if nbCompo in dictVars:
 223                 dictVars[nbCompo].append(mdf)
 224                 pass
 225             else:
 226                 dictVars[nbCompo]=[mdf]
 227                 pass
 228             pass
 229         for mdf in mdfs:
 230             nbCompo=mdf.getNumberOfComponents()
 231             if nbCompo not in self.dictCompo:
 232                 l = [x for x in self.dictCompo if x - nbCompo > 0]
 233                 if len(l)==0:
 234                     continue;
 235                 nbCompo=l[0]
 236                 pass
 237             for iii,it in enumerate(its):
 238                 ff=mdf[it]
 239                 isMultiDisc=len(ff.getTypesOfFieldAvailable())>1
 240                 for typ in ff.getTypesOfFieldAvailable():
 241                     l=self._l[:] ; l[-1]="%s%s.%s"%(self._base_name_without_dir,str(iii).rjust(4,"0"),ff.getName())
 242                     if isMultiDisc:
 243                         l[-1]="%s_%s"(l[-1],MEDCouplingFieldDiscretization.New(typ).getStringRepr())
 244                         pass
 245                     fffn=l[-1]
 246                     try:
 247                         os.remove(os.path.sep.join(l))
 248                     except:
 249                         pass
 250                     f=open(os.path.sep.join(l),"w+b")
 251                     summ=0
 252                     for geo,[(curTyp,(bg,end),pfl,loc)] in ff.getFieldSplitedByType():
 253                         if typ==curTyp:
 254                             summ+=4*nbCompo*(end-bg)+80
 255                             pass
 256                         pass
 257                     a=np.memmap(f,dtype='byte',mode='w+',offset=0,shape=(2*80+4+summ,)) ; a.flush() # truncate to set the size of the file
 258                     mm=mmap.mmap(f.fileno(),offset=0,length=0)
 259                     k1=ff.getName()
 260                     if isMultiDisc:
 261                         k1="%s_%s"%(k1,MEDCouplingFieldDiscretization.New(typ).getStringRepr())
 262                         pass
 263                     mm.write(self.__str80(k1))
 264                     mm.write(self.__str80("part"))
 265                     a=np.memmap(f,dtype='int32',mode='w+',offset=mm.tell(),shape=(1,))
 266                     a[0]=1 ; a.flush() ; mm.seek(mm.tell()+4) # part number maybe to change ?
 267                     for geo,[(curTyp,(bg,end),pfl,loc)] in ff.getFieldSplitedByType():
 268                         if pfl!="":
 269                             raise Exception("Field \"%s\" contains profiles ! Profiles are not supported yet !"%(mdf.getName()))
 270                         if typ==curTyp:
 271                             arr=ff.getUndergroundDataArray()[bg:end].changeNbOfComponents(nbCompo,0.) ; arr=arr.toNoInterlace()
 272                             if typ==ON_CELLS:
 273                                 mm.write(self.__str80(self.dictMCTyp_str[geo]))
 274                                 pass
 275                             elif typ==ON_NODES:
 276                                 mm.write(self.__str80("coordinates"))
 277                                 pass
 278                             else:
 279                                 print(("UnManaged type of field for field \"%s\" !"%(mdf.getName())))
 280                                 pass
 281                             a=np.memmap(f,dtype='float32',mode='w+',offset=mm.tell(),shape=(nbCompo,end-bg))
 282                             b=arr.toNumPyArray() ; b=b.reshape(nbCompo,end-bg)
 283                             a[:]=b
 284                             a.flush() ; mm.seek(mm.tell()+nbCompo*(end-bg)*4)
 285                             pass
 286                         pass
 287                     k="%s per %s"%(self.dictCompo[nbCompo],self.discSpatial[typ])
 288                     if k in self._ze_top_dict:
 289                         if k1 in self._ze_top_dict[k]:
 290                             self._ze_top_dict[k][k1].append(fffn)
 291                             pass
 292                         else:
 293                             self._ze_top_dict[k][k1]=[fffn]
 294                             pass
 295                     else:
 296                         self._ze_top_dict[k]={k1:[fffn]}
 297                         pass
 298                     pass
 299                 pass
 300             pass
 301         headerPart=""
 302         if len(self._ze_top_dict)!=0:
 303             hvp=self.header_varpart[:]
 304             for k in self._ze_top_dict:
 305                 for k1 in self._ze_top_dict[k]:
 306                     hvp+="\n%s: %s %s"%(k,k1,re.sub("([\d]{4})",4*"*",self._ze_top_dict[k][k1][0]))
 307                     pass
 308                 pass
 309             hvp+="\n"
 310             headerPart+=hvp
 311             #
 312             ddd={"NbTimeSteps":len(its),"TimeValues":TimeValues}
 313             htp=self.header_timepart%ddd
 314             headerPart+=htp
 315             pass
 316         return headerPart
 317
 318     @classmethod
 319     def __str80(cls,st):
 320         if len(st)>79:
 321             raise Exception("String \"%s\" is too long (>79) !"%(st))
 322         return bytes(str(st).ljust(79)+"\n", "ascii")
 323
 324     def __computeSizeOfGeoFile(self,listOfMeshes,nn):
 325         sz=0
 326         for m in listOfMeshes:
 327             distribTypes=m.getDistributionOfTypes()
 328             sz+=80+4+2*80+4+nn*3*4
 329             i=0
 330             for typ2,nbelem,dummy in distribTypes:
 331                 typ=typ2
 332                 if typ not in self.dictMCTyp:
 333                     typ=MEDCouplingMesh.GetCorrespondingPolyType()
 334                     pass
 335                 if typ!=NORM_POLYHED and typ!=NORM_POLYGON:
 336                     sz+=80+4+MEDCouplingMesh.GetNumberOfNodesOfGeometricType(typ)*nbelem*4
 337                     pass
 338                 elif typ==NORM_POLYHED:
 339                     mplh=m[i:i+nbelem] ; delta=len(mplh.getNodalConnectivity())+nbelem
 340                     sz+=80+4+delta*4
 341                     pass
 342                 else:
 343                     mplh=m[i:i+nbelem] ; delta=len(mplh.getNodalConnectivity())
 344                     sz+=80+4+delta*4
 345                     pass
 346                 i+=nbelem
 347                 pass
 348             pass
 349         return sz