src/ParaMEDMEM/MEDLoader/MEDLoader.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #include "MEDLoader.hxx"
  20 #include "CellModel.hxx"
  21 #include "ParaMESH.hxx"
  22 #include "BlockTopology.hxx"
  23 #include "MEDCouplingUMesh.hxx"
  24
  25 extern "C"
  26 {
  27 #include "med.h"
  28 }
  29
  30 #include <string>
  31 #include <cstring>
  32 #include <sstream>
  33 #include <fstream>
  34
  35 med_geometrie_element typmai[MED_NBR_GEOMETRIE_MAILLE] = { MED_POINT1,
  36                                                            MED_SEG2,
  37                                                            MED_SEG3,
  38                                                            MED_TRIA3,
  39                                                            MED_TRIA6,
  40                                                            MED_QUAD4,
  41                                                            MED_QUAD8,
  42                                                            MED_TETRA4,
  43                                                            MED_TETRA10,
  44                                                            MED_HEXA8,
  45                                                            MED_HEXA20,
  46                                                            MED_PENTA6,
  47                                                            MED_PENTA15,
  48                                                            MED_PYRA5,
  49                                                            MED_PYRA13 };
  50
  51 INTERP_KERNEL::NormalizedCellType typmai2[MED_NBR_GEOMETRIE_MAILLE] = { INTERP_KERNEL::NORM_ERROR,
  52                                                                         INTERP_KERNEL::NORM_SEG2,
  53                                                                         INTERP_KERNEL::NORM_SEG3,
  54                                                                         INTERP_KERNEL::NORM_TRI3,
  55                                                                         INTERP_KERNEL::NORM_TRI6,
  56                                                                         INTERP_KERNEL::NORM_QUAD4,
  57                                                                         INTERP_KERNEL::NORM_QUAD8,
  58                                                                         INTERP_KERNEL::NORM_TETRA4,
  59                                                                         INTERP_KERNEL::NORM_TETRA10,
  60                                                                         INTERP_KERNEL::NORM_HEXA8,
  61                                                                         INTERP_KERNEL::NORM_HEXA20,
  62                                                                         INTERP_KERNEL::NORM_PENTA6,
  63                                                                         INTERP_KERNEL::NORM_PENTA15,
  64                                                                         INTERP_KERNEL::NORM_PYRA5,
  65                                                                         INTERP_KERNEL::NORM_PYRA13 };
  66
  67 using namespace ParaMEDMEM;
  68
  69 const char WHITE_SPACES[]=" \n";
  70
  71 MEDLoader::MEDConnOfOneElemType::MEDConnOfOneElemType(INTERP_KERNEL::NormalizedCellType type, int *conn, int lgth):_lgth(lgth),
  72                                                                                                                    _conn(conn),_global(0),
  73                                                                                                                    _type(type)
  74 {
  75 }
  76
  77 void MEDLoader::MEDConnOfOneElemType::setGlobal(int *global)
  78 {
  79   if(_global!=global)
  80     {
  81       if(_global)
  82         delete [] _global;
  83       _global=global;
  84     }
  85 }
  86
  87 void MEDLoader::MEDConnOfOneElemType::releaseArray()
  88 {
  89   delete [] _conn;
  90   delete [] _global;
  91 }
  92
  93 std::string buildStringFromFortran(const char *expr, int lgth)
  94 {
  95   std::string ret(expr,lgth);
  96   std::string whiteSpaces(WHITE_SPACES);
  97   std::size_t lgthReal=strlen(ret.c_str());
  98   std::string ret2=ret.substr(0,lgthReal);
  99   std::size_t found=ret2.find_last_not_of(whiteSpaces);
 100   if (found!=std::string::npos)
 101     ret2.erase(found+1);
 102   else
 103     ret2.clear();//ret is all whitespace
 104   return ret2;
 105 }
 106
 107 namespace MEDLoader
 108 {
 109   med_int getIdFromMeshName(med_idt fid, const char *meshName) throw(INTERP_KERNEL::Exception)
 110   {
 111     if(meshName==0)
 112       return 1;
 113     med_int n=MEDnMaa(fid);
 114     if(n==0)
 115       throw INTERP_KERNEL::Exception("No mesh in file.");
 116     med_maillage type_maillage;
 117     char maillage_description[MED_TAILLE_DESC+1];
 118     med_int dim;
 119     char nommaa[MED_TAILLE_NOM+1];
 120     std::ostringstream os;
 121     for(med_int i=1;i<=n;i++)
 122       {
 123         MEDmaaInfo(fid,i,nommaa,&dim,&type_maillage,maillage_description);
 124         std::string cur=buildStringFromFortran(nommaa,sizeof(nommaa));
 125         if(cur==meshName)
 126           return i;
 127         os << "\'" << cur.c_str() << "\' ";
 128       }
 129     std::ostringstream os2;
 130     os2 << "MeshName '" << meshName << "' not in file : meshes available : " << os.str();
 131     throw INTERP_KERNEL::Exception(os2.str().c_str());
 132   }
 133
 134   void readUMeshDataInMedFile(med_idt fid, med_int meshId, double *&coords, int& nCoords, int& spaceDim, std::list<MEDLoader::MEDConnOfOneElemType>& conn)
 135   {
 136     char nommaa[MED_TAILLE_NOM+1];
 137     char maillage_description[MED_TAILLE_DESC+1];
 138     char comp[3*MED_TAILLE_PNOM+1];
 139     char unit[3*MED_TAILLE_PNOM+1];
 140     med_maillage type_maillage;
 141     med_int Mdim;
 142     MEDmaaInfo(fid,meshId,nommaa,&Mdim,&type_maillage,maillage_description);
 143     spaceDim=(int)Mdim;
 144     nCoords=MEDnEntMaa(fid,nommaa,MED_COOR,MED_NOEUD,(med_geometrie_element)0,(med_connectivite)0);
 145     coords=new double[nCoords*spaceDim];
 146     med_repere repere;
 147     MEDcoordLire(fid,nommaa,Mdim,coords,MED_FULL_INTERLACE,MED_ALL,NULL,0,&repere,comp,unit);
 148     med_booleen inoele, inuele;
 149     for(int i=0;i<MED_NBR_GEOMETRIE_MAILLE;i++)
 150       {
 151         med_geometrie_element curMedType=typmai[i];
 152         med_entite_maillage whichEntity;
 153         int curNbOfElemM=MEDnEntMaa(fid,nommaa,MED_CONN,MED_MAILLE,curMedType,MED_NOD);
 154         int curNbOfElemF=MEDnEntMaa(fid,nommaa,MED_CONN,MED_FACE,curMedType,MED_NOD);
 155         int curNbOfElem;
 156         if(curNbOfElemF>curNbOfElemM)
 157           {
 158             curNbOfElem=curNbOfElemF;
 159             whichEntity=MED_FACE;
 160           }
 161         else
 162           {
 163             curNbOfElem=curNbOfElemM;
 164             whichEntity=MED_MAILLE;
 165           }
 166         if(curNbOfElem>0)
 167           {
 168             int *connTab=new int[(curMedType%100)*curNbOfElem];
 169             MEDLoader::MEDConnOfOneElemType elem(typmai2[i],connTab,curNbOfElem);
 170             int *tmp=new int[curNbOfElem];
 171             int *fam=new int[curNbOfElem];
 172             char *noms=new char[MED_TAILLE_PNOM*curNbOfElem+1];
 173             MEDelementsLire(fid,nommaa,Mdim,connTab,MED_FULL_INTERLACE,noms,&inoele,tmp,&inuele,fam,curNbOfElem,whichEntity,curMedType,MED_NOD);
 174             delete [] tmp;
 175             delete [] fam;
 176             delete [] noms;
 177             //trying to read global numbering
 178             int *globArr=new int[curNbOfElem];
 179             if(MEDglobalNumLire(fid,nommaa,globArr,curNbOfElem,whichEntity,curMedType)==0)
 180               elem.setGlobal(globArr);
 181             else
 182               delete [] globArr;
 183             conn.push_back(elem);
 184           }
 185       }
 186   }
 187 }
 188
 189 unsigned MEDLoader::calculateHighestMeshDim(const std::list<MEDLoader::MEDConnOfOneElemType>& conn)
 190 {
 191   unsigned ret=0;
 192   for(std::list<MEDLoader::MEDConnOfOneElemType>::const_iterator iter=conn.begin();iter!=conn.end();iter++)
 193     {
 194       unsigned curDim=INTERP_KERNEL::CellModel::getCellModel((*iter).getType()).getDimension();
 195       if(ret<curDim)
 196         ret=curDim;
 197     }
 198   return ret;
 199 }
 200
 201 void MEDLoader::keepSpecifiedMeshDim(std::list<MEDLoader::MEDConnOfOneElemType>& conn, unsigned meshDim)
 202 {
 203   for(std::list<MEDLoader::MEDConnOfOneElemType>::iterator iter=conn.begin();iter!=conn.end();)
 204     {
 205       unsigned curDim=INTERP_KERNEL::CellModel::getCellModel((*iter).getType()).getDimension();
 206       if(curDim!=meshDim)
 207         {
 208           (*iter).releaseArray();
 209           iter=conn.erase(iter);
 210         }
 211       else
 212         iter++;
 213     }
 214 }
 215
 216 void MEDLoader::tradMEDFileCoreFrmt2MEDCouplingUMesh(const std::list<MEDLoader::MEDConnOfOneElemType>& medConnFrmt,
 217                                                      DataArrayInt* &conn,
 218                                                      DataArrayInt* &connIndex)
 219 {
 220   if(medConnFrmt.empty())
 221     {
 222       conn=0;
 223       connIndex=0;
 224       return ;
 225     }
 226   std::list<MEDLoader::MEDConnOfOneElemType>::const_iterator iter=medConnFrmt.begin();
 227   int totalNbOfCells=0;
 228   int totalNbOfMedConn=0;
 229   for(;iter!=medConnFrmt.end();iter++)
 230     {
 231       const INTERP_KERNEL::CellModel& cellMod=INTERP_KERNEL::CellModel::getCellModel((*iter).getType());
 232       totalNbOfCells+=(*iter).getLength();
 233       if(!cellMod.isDynamic())
 234         totalNbOfMedConn+=(*iter).getLength()*cellMod.getNumberOfNodes();
 235       else
 236         throw INTERP_KERNEL::Exception("Polyg/polh not implemented yet !");
 237     }
 238   connIndex=DataArrayInt::New();
 239   conn=DataArrayInt::New();
 240   connIndex->alloc(totalNbOfCells+1,1);
 241   int *connIdxPtr=connIndex->getPointer();
 242   int connFillId=0;
 243   conn->alloc(totalNbOfMedConn+totalNbOfCells,1);
 244   int *connPtr=conn->getPointer();
 245   for(iter=medConnFrmt.begin();iter!=medConnFrmt.end();iter++)
 246     {
 247       INTERP_KERNEL::NormalizedCellType type=(*iter).getType();
 248       int *sourceConn=(*iter).getArray();
 249       const INTERP_KERNEL::CellModel& cellMod=INTERP_KERNEL::CellModel::getCellModel(type);
 250       int nbOfCellsInCurType;
 251       int nbOfNodesIn1Cell=cellMod.getNumberOfNodes();
 252       if(!cellMod.isDynamic())
 253         nbOfCellsInCurType=(*iter).getLength();
 254       else
 255         throw INTERP_KERNEL::Exception("Polyg/polh not implemented yet !");
 256       if(!cellMod.isDynamic())
 257         {
 258           for(int i=0;i<nbOfCellsInCurType;i++,connIdxPtr++)
 259             {
 260               *connIdxPtr=connFillId;
 261               *connPtr++=type;
 262               connPtr=std::transform(sourceConn,sourceConn+nbOfNodesIn1Cell,connPtr,std::bind2nd(std::minus<int>(),1));
 263               connFillId+=nbOfNodesIn1Cell+1;
 264               sourceConn+=nbOfNodesIn1Cell;
 265             }
 266           *connIdxPtr=connFillId;
 267         }
 268     }
 269 }
 270
 271 void MEDLoader::releaseMEDFileCoreFrmt(std::list<MEDLoader::MEDConnOfOneElemType>& medConnFrmt)
 272 {
 273   for(std::list<MEDLoader::MEDConnOfOneElemType>::iterator iter=medConnFrmt.begin();iter!=medConnFrmt.end();iter++)
 274     (*iter).releaseArray();
 275   medConnFrmt.clear();
 276 }
 277
 278 /*!
 279  * This method builds a sub set of connectivity for a given type 'type'.
 280  * @param conn input containing connectivity with MEDCoupling format.
 281  * @param connIndex input containing connectivity index in MEDCoupling format.
 282  * @param type input specifying which cell types will be extracted in conn4MEDFile.
 283  * @param conn4MEDFile output containing the connectivity directly understandable by MEDFile; conn4MEDFile has to be empty before this method called.
 284  * @return nb of elements extracted.
 285  */
 286 int MEDLoader::buildMEDSubConnectivityOfOneType(DataArrayInt *conn, DataArrayInt *connIndex, INTERP_KERNEL::NormalizedCellType type, std::vector<int>& conn4MEDFile)
 287 {
 288   int ret=0;
 289   int nbOfElem=connIndex->getNbOfElems()-1;
 290   const int *connPtr=conn->getPointer();
 291   const int *connIdxPtr=connIndex->getPointer();
 292   for(int i=0;i<nbOfElem;i++)
 293     {
 294       int delta=connIdxPtr[1]-connIdxPtr[0];
 295       if(*connPtr==type)
 296         {
 297           conn4MEDFile.insert(conn4MEDFile.end(),connPtr+1,connPtr+delta);
 298           ret++;
 299         }
 300       connIdxPtr++;
 301       connPtr+=delta;
 302     }
 303   std::transform(conn4MEDFile.begin(),conn4MEDFile.end(),conn4MEDFile.begin(),std::bind2nd(std::plus<int>(),1));
 304   return ret;
 305 }
 306
 307 MEDCouplingUMesh *MEDLoader::ReadUMeshFromFile(const char *fileName, const char *meshName, int meshDimRelToMax) throw(INTERP_KERNEL::Exception)
 308 {
 309   //Extraction data from MED file.
 310   med_idt fid=MEDouvrir((char *)fileName,MED_LECTURE);
 311   med_int mid=getIdFromMeshName(fid,meshName);
 312   unsigned meshDimExtract;
 313   double *coords;
 314   int nCoords;
 315   int spaceDim;
 316   std::list<MEDLoader::MEDConnOfOneElemType> conn;
 317   readUMeshDataInMedFile(fid,mid,coords,nCoords,spaceDim,conn);
 318   meshDimExtract=calculateHighestMeshDim(conn);
 319   meshDimExtract=meshDimExtract+meshDimRelToMax;
 320   keepSpecifiedMeshDim(conn,meshDimExtract);
 321   MEDfermer(fid);
 322   //Put data in returned data structure.
 323   MEDCouplingUMesh *ret=MEDCouplingUMesh::New();
 324   ret->setName(meshName);
 325   ret->setMeshDimension(meshDimExtract);
 326   //
 327   DataArrayDouble *coordsArr=DataArrayDouble::New();
 328   coordsArr->useArray(coords,true,ParaMEDMEM::CPP_DEALLOC,nCoords,spaceDim);
 329   ret->setCoords(coordsArr);
 330   coordsArr->decrRef();
 331   //
 332   DataArrayInt *connArr,*connIndexArr;
 333   tradMEDFileCoreFrmt2MEDCouplingUMesh(conn,connArr,connIndexArr);
 334   ret->setConnectivity(connArr,connIndexArr);
 335   //clean-up
 336   if(connArr)
 337     connArr->decrRef();
 338   if(connIndexArr)
 339     connIndexArr->decrRef();
 340   releaseMEDFileCoreFrmt(conn);
 341   return ret;
 342 }
 343
 344 void MEDLoader::writeUMesh(const char *fileName, ParaMEDMEM::MEDCouplingUMesh *mesh)
 345 {
 346   med_idt fid=MEDouvrir((char *)fileName,MED_CREATION);
 347   char maa[MED_TAILLE_NOM+1];
 348   std::fill(maa,maa+MED_TAILLE_NOM+1,'\0');
 349   const char *meshName=mesh->getName();
 350   strcpy(maa,meshName);
 351   MEDmaaCr(fid,maa,mesh->getMeshDimension(),MED_NON_STRUCTURE,maa);
 352   std::set<INTERP_KERNEL::NormalizedCellType> allTypes(mesh->getAllTypes());
 353   DataArrayInt *conn=mesh->getNodalConnectivity();
 354   DataArrayInt *connIndex=mesh->getNodalConnectivityIndex();
 355   char familyName[MED_TAILLE_NOM+1];
 356   std::fill(familyName,familyName+MED_TAILLE_NOM+1,'\0');
 357   const char DftFamilyName[]="DftFamily";
 358   std::copy(DftFamilyName,DftFamilyName+sizeof(DftFamilyName),familyName);
 359   for(int i=0;i<MED_NBR_GEOMETRIE_MAILLE;i++)
 360     {
 361       med_geometrie_element curMedType=typmai[i];
 362       INTERP_KERNEL::NormalizedCellType curType=typmai2[i];
 363       if(allTypes.find(curType)!=allTypes.end())
 364         {
 365           std::vector<int> medConn;
 366           int nbOfElt=buildMEDSubConnectivityOfOneType(conn,connIndex,curType,medConn);
 367           MEDconnEcr(fid,maa,mesh->getMeshDimension(),&medConn[0],MED_FULL_INTERLACE,nbOfElt,MED_MAILLE,curMedType,MED_NOD);
 368         }
 369     }
 370   MEDfamCr(fid,maa,familyName,0,0,0,0,0,0,0);
 371   DataArrayDouble *arr=mesh->getCoords();
 372   char comp[2*MED_TAILLE_PNOM+1];
 373   char unit[2*MED_TAILLE_PNOM+1];
 374   std::fill(comp,comp+2*MED_TAILLE_PNOM,' ');
 375   comp[2*MED_TAILLE_PNOM]='\0';
 376   char *work=comp;
 377   for(int i=0;i<mesh->getSpaceDimension();i++,work+=3)
 378     *work='X'+i;
 379   std::fill(unit,unit+2*MED_TAILLE_PNOM+1,'\0');
 380   MEDcoordEcr(fid,maa,mesh->getSpaceDimension(),arr->getPointer(),MED_FULL_INTERLACE,mesh->getNumberOfNodes(),MED_CART,comp,unit);
 381   MEDfermer(fid);
 382 }
 383
 384 /*!
 385  * This method builds the master file 'fileName' of a parallel MED file defined in 'fileNames'.
 386  */
 387 void MEDLoader::writeMasterFile(const char *fileName, const std::vector<std::string>& fileNames, const char *meshName)
 388 {
 389   int nbOfDom=fileNames.size();
 390   std::ofstream fs(fileName);
 391   fs << "#MED Fichier V 2.3" << " " << std::endl;
 392   fs << "#"<<" " << std::endl;
 393   fs << nbOfDom <<" " << std::endl;
 394   for(int i=0;i<nbOfDom;i++)
 395     fs << meshName << " " << i+1 << " " << meshName << "_" << i+1 << " localhost " << fileNames[i] << " " << std::endl;
 396 }
 397
 398 void MEDLoader::writeParaMesh(const char *fileName, ParaMEDMEM::ParaMESH *mesh)
 399 {
 400   if(!mesh->getBlockTopology()->getProcGroup()->containsMyRank())
 401     return ;
 402   int myRank=mesh->getBlockTopology()->getProcGroup()->myRank();
 403   int nbDomains=mesh->getBlockTopology()->getProcGroup()->size();
 404   std::vector<std::string> fileNames(nbDomains);
 405   for(int i=0;i<nbDomains;i++)
 406     {
 407       std::ostringstream sstr;
 408       sstr << fileName << i+1 << ".med";
 409       fileNames[i]=sstr.str();
 410     }
 411   if(myRank==0)
 412     writeMasterFile(fileName,fileNames,mesh->getCellMesh()->getName());
 413   writeUMesh(fileNames[myRank].c_str(),mesh->getCellMesh());
 414 }
 415
 416 void MEDLoader::writeParaField(const char *fileName, const char *meshName, ParaMEDMEM::ParaFIELD *f)
 417 {
 418 }