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