src/MEDPartitioner/MEDPARTITIONER_MeshCollectionDriver.cxx

   1 // Copyright (C) 2007-2015  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, or (at your option) any later version.
   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
  20 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
  21
  22 #include "MEDPARTITIONER_ConnectZone.hxx"
  23 #include "MEDPARTITIONER_MeshCollection.hxx"
  24 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
  25 #include "MEDPARTITIONER_ParallelTopology.hxx"
  26 #include "MEDPARTITIONER_Utils.hxx"
  27
  28 #include "MEDCouplingFieldDouble.hxx"
  29 #include "MEDCouplingRefCountObject.hxx"
  30 #include "MEDCouplingUMesh.hxx"
  31 #include "MEDFileField.hxx"
  32 #include "MEDFileJoint.hxx"
  33 #include "MEDFileMesh.hxx"
  34 #include "MEDLoader.hxx"
  35
  36 #include <map>
  37 #include <set>
  38 #include <vector>
  39 #include <string>
  40 #include <fstream>
  41 #include <iostream>
  42
  43 #include <libxml/tree.h>
  44 #include <libxml/parser.h>
  45 #include <libxml/xpath.h>
  46 #include <libxml/xpathInternals.h>
  47
  48 #include "med.h"
  49
  50 using namespace MEDPARTITIONER;
  51
  52 MeshCollectionDriver::MeshCollectionDriver(MeshCollection* collection):_collection(collection)
  53 {
  54 }
  55
  56 /*!reads a unique MED File v>=2.1
  57  * and mounts the corresponding mesh in memory
  58  *\param filename binary file
  59  *\param meshname mesh name in the MED file
  60  * */
  61 int MeshCollectionDriver::readSeq(const char* filename, const char* meshname)
  62 {
  63   std::cout << "readSeq" << std::endl;
  64   MyGlobals::_File_Names.resize(1);
  65   MyGlobals::_File_Names[0]=std::string(filename);
  66
  67   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(filename,meshname);
  68   //puts the only mesh in the mesh vector
  69   (_collection->getMesh()).push_back(mfm->getLevel0Mesh(false));
  70   (_collection->getFaceMesh()).push_back(mfm->getLevelM1Mesh(false));
  71
  72   //reading family ids
  73   ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
  74   ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
  75   (_collection->getCellFamilyIds()).push_back(cellIds);
  76   (_collection->getFaceFamilyIds()).push_back(faceIds);
  77
  78   //reading groups
  79   (_collection->getFamilyInfo())=mfm->getFamilyInfo();
  80   (_collection->getGroupInfo())=mfm->getGroupInfo();
  81
  82   (_collection->getCZ()).clear();
  83
  84   ParallelTopology* aPT = new ParallelTopology((_collection->getMesh()));
  85   _collection->setTopology(aPT, true);
  86   _collection->setName(meshname);
  87   _collection->setDomainNames(meshname);
  88   return 0;
  89 }
  90
  91
  92 void MeshCollectionDriver::readMEDFileData(const ParaMEDMEM::MEDFileData* filedata)
  93 {
  94   const int nbDomains = filedata->getMeshes()->getNumberOfMeshes();
  95   _collection->getMesh()         .resize( nbDomains, 0 );
  96   _collection->getFaceMesh()     .resize( nbDomains, 0 );
  97   _collection->getCellFamilyIds().resize( nbDomains, 0 );
  98   _collection->getFaceFamilyIds().resize( nbDomains, 0 );
  99
 100   for (int i=0; i<nbDomains; i++)
 101     {
 102       ParaMEDMEM::MEDFileUMesh *mfm = dynamic_cast<ParaMEDMEM::MEDFileUMesh *>(filedata->getMeshes()->getMeshAtPos(i));
 103       readData(mfm,i);
 104       if ( mfm && mfm->getMeshDimension() > 0 )
 105         _collection->setNonEmptyMesh( i );
 106     }
 107
 108   ParallelTopology* aPT = new ParallelTopology(_collection->getMesh());
 109   _collection->setTopology(aPT, true);
 110   if ( nbDomains > 0 )
 111     {
 112       _collection->setName( filedata->getMeshes()->getMeshAtPos(0)->getName() );
 113       _collection->setDomainNames( _collection->getName() );
 114     }
 115   if ( ParaDomainSelector* domainSelector = _collection->getParaDomainSelector() )
 116     if ( _collection->isParallelMode() )
 117       {
 118         //to know nb of cells on each proc to compute global cell ids from locally global
 119         domainSelector->gatherNbOf(_collection->getMesh());
 120       }
 121 }
 122
 123 void MeshCollectionDriver::readFileData(std::string file,std::string meshname,int idomain) const
 124 {
 125   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file,meshname);
 126   readData(mfm,idomain);
 127   mfm->decrRef();
 128 }
 129
 130 void MeshCollectionDriver::readData(ParaMEDMEM::MEDFileUMesh* mfm, int idomain) const
 131 {
 132   std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
 133   try
 134     {
 135       (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false);
 136       //reading families groups
 137       ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
 138       (_collection->getCellFamilyIds())[idomain]=cellIds;
 139     }
 140   catch(...)
 141     {
 142       (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
 143       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
 144       empty->alloc(0,1);
 145       (_collection->getCellFamilyIds())[idomain]=empty;
 146       std::cout<<"\nNO Level0Mesh (Cells)\n";
 147     }
 148   try
 149     {
 150       if (nonEmpty.size()>1 && nonEmpty[1]==-1)
 151         {
 152           (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
 153           //reading families groups
 154           ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
 155           (_collection->getFaceFamilyIds())[idomain]=faceIds;
 156           if (MyGlobals::_Verbose>10)
 157             std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
 158         }
 159       else
 160         {
 161           throw INTERP_KERNEL::Exception("no faces");
 162         }
 163     }
 164   catch(...)
 165     {
 166       (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
 167       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
 168       (_collection->getFaceFamilyIds())[idomain]=empty;
 169       if (MyGlobals::_Verbose>10)
 170         std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
 171     }
 172   //reading groups
 173   _collection->getFamilyInfo()=mfm->getFamilyInfo();
 174   _collection->getGroupInfo()=mfm->getGroupInfo();
 175 }
 176
 177 void MeshCollectionDriver::readSubdomain(int idomain)
 178 {
 179   std::string meshname=MyGlobals::_Mesh_Names[idomain];
 180   std::string file=MyGlobals::_File_Names[idomain];
 181   readFileData(file,meshname,idomain);
 182
 183   std::vector<std::string> localInformation;
 184   std::string str;
 185   localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
 186   localInformation.push_back(str+"meshName="+meshname);
 187   MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
 188   std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
 189   if (localFields.size()>0)
 190     MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
 191 }
 192
 193 ParaMEDMEM::MEDFileMesh* MeshCollectionDriver::getMesh(int idomain) const
 194 {
 195   ParaMEDMEM::MEDFileUMesh* mfm = ParaMEDMEM::MEDFileUMesh::New();
 196
 197   ParaMEDMEM::MEDCouplingUMesh* cellMesh=_collection->getMesh(idomain);
 198   ParaMEDMEM::MEDCouplingUMesh* faceMesh=_collection->getFaceMesh(idomain);
 199   // std::string cleFilter=Cle1ToStr("filterFaceOnCell",idomain);
 200   // ParaMEDMEM::DataArrayInt* filter=0;
 201   // if (_collection->getMapDataArrayInt().find(cleFilter)!=_collection->getMapDataArrayInt().end())
 202   //   {
 203   //     filter=_collection->getMapDataArrayInt().find(cleFilter)->second;
 204   //     int* index=filter->getPointer();
 205   //     faceMeshFilter=(ParaMEDMEM::MEDCouplingUMesh *) faceMesh->buildPartOfMySelf(index,index+filter->getNbOfElems(),true);
 206   //     faceMesh=faceMeshFilter;
 207   //   }
 208   // if (faceMeshFilter!=0)
 209   //   faceMeshFilter->decrRef();
 210   std::string finalMeshName="";
 211   if (MyGlobals::_General_Informations.size()!=0)
 212     {
 213       std::size_t found=MyGlobals::_General_Informations[0].find("finalMeshName=");
 214       if ((found!=std::string::npos) && (found>0))
 215         {
 216           finalMeshName=ExtractFromDescription(MyGlobals::_General_Informations[0], "finalMeshName=");
 217         }
 218     }
 219   if (finalMeshName.empty())
 220     {
 221       finalMeshName=_collection->getName();
 222     }
 223   cellMesh->setName(finalMeshName);
 224   mfm->setMeshAtLevel( 0, cellMesh );
 225
 226   faceMesh->checkCoherency();
 227   if (faceMesh->getNumberOfCells()>0)
 228     {
 229       faceMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10);
 230       faceMesh->setName(finalMeshName);
 231       mfm->setMeshAtLevel( -1, faceMesh );
 232     }
 233
 234   // ParaMEDMEM::MEDCouplingUMesh* boundaryMesh=0;
 235   // if (MyGlobals::_Creates_Boundary_Faces>0)
 236   //   {
 237   //     //try to write Boundary meshes
 238   //     bool keepCoords=false; //TODO or true
 239   //     boundaryMesh=(ParaMEDMEM::MEDCouplingUMesh *) cellMesh->buildBoundaryMesh(keepCoords);
 240   //     boundaryMesh->setName("boundaryMesh");
 241   // if (boundaryMesh!=0)
 242   //   {
 243   //     //doing that testMesh becomes second mesh sorted by alphabetical order of name
 244   //     MEDLoader::WriteUMesh(distfilename, boundaryMesh, false);
 245   //     boundaryMesh->decrRef();
 246   //   }
 247
 248   mfm->setFamilyInfo(_collection->getFamilyInfo());
 249   mfm->setGroupInfo(_collection->getGroupInfo());
 250   std::string key=Cle1ToStr("faceFamily_toArray",idomain);
 251   if ( faceMesh->getNumberOfCells()>0 && _collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
 252     mfm->setFamilyFieldArr(-1,_collection->getMapDataArrayInt().find(key)->second);
 253   key=Cle1ToStr("cellFamily_toArray",idomain);
 254   if (_collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
 255     mfm->setFamilyFieldArr(0,_collection->getMapDataArrayInt().find(key)->second);
 256
 257   // add joints
 258
 259   using ParaMEDMEM::MEDCouplingAutoRefCountObjectPtr;
 260   using ParaMEDMEM::MEDCouplingSkyLineArray;
 261   using ParaMEDMEM::MEDFileJoint;
 262   using ParaMEDMEM::MEDFileJointCorrespondence;
 263   using ParaMEDMEM::MEDFileJointOneStep;
 264   using ParaMEDMEM::MEDFileJoints;
 265   using ParaMEDMEM::MEDFileJoints;
 266
 267   if ( _collection->getCZ().size() > 0 )
 268     {
 269       MEDCouplingAutoRefCountObjectPtr< MEDFileJoints > joints = MEDFileJoints::New();
 270
 271       for ( size_t i = 0; i < _collection->getCZ().size(); ++i )
 272         {
 273           ConnectZone* cz = _collection->getCZ()[i];
 274           if ( !cz ||
 275                cz->getLocalDomainNumber() != idomain )
 276             continue;
 277           {
 278             std::ostringstream oss;
 279             oss << "joint_" << cz->getDistantDomainNumber();
 280             cz->setName( oss.str() );
 281           }
 282           {
 283             std::ostringstream oss;
 284             oss << "connect_zone_" << i;
 285             cz->setDescription( oss.str() );
 286           }
 287
 288           MEDCouplingAutoRefCountObjectPtr< MEDFileJoint>
 289             joint = MEDFileJoint::New( cz->getName(), finalMeshName,
 290                                        finalMeshName, cz->getDistantDomainNumber() );
 291           joint->setDescription( cz->getDescription() );
 292           joints->pushJoint( joint );
 293
 294           MEDCouplingAutoRefCountObjectPtr< MEDFileJointOneStep> j1st = MEDFileJointOneStep::New();
 295           joint->pushStep( j1st );
 296
 297           const MEDCouplingSkyLineArray * nodeCorr = cz->getNodeCorresp();
 298           if ( nodeCorr )
 299             {
 300               MEDCouplingAutoRefCountObjectPtr< MEDFileJointCorrespondence >
 301                 corr = MEDFileJointCorrespondence::New( nodeCorr->getValueArray() );
 302               j1st->pushCorrespondence( corr );
 303             }
 304
 305           std::vector< std::pair< int,int > > types = cz->getEntities();
 306           INTERP_KERNEL::NormalizedCellType t1, t2;
 307           for ( size_t it = 0; it < types.size(); ++it )
 308             {
 309               const MEDCouplingSkyLineArray * cellCorr =
 310                 cz->getEntityCorresp( types[it].first, types[it].second );
 311               if ( cellCorr && cellCorr->getNumberOf() > 0 )
 312                 {
 313                   t1 = INTERP_KERNEL::NormalizedCellType( types[it].first );
 314                   t2 = INTERP_KERNEL::NormalizedCellType( types[it].second );
 315                   MEDCouplingAutoRefCountObjectPtr< MEDFileJointCorrespondence>
 316                     corr = MEDFileJointCorrespondence::New( cellCorr->getValueArray(), t1, t2 );
 317                   j1st->pushCorrespondence( corr );
 318                 }
 319             }
 320         }
 321       mfm->setJoints( joints );
 322     }
 323
 324   return mfm;
 325 }
 326
 327 ParaMEDMEM::MEDCouplingFieldDouble* MeshCollectionDriver::getField(std::string key, std::string description, ParaMEDMEM::DataArrayDouble* data, ParaMEDMEM::MEDFileMesh* mfm, int idomain) const
 328 {
 329   std::string desc=description;
 330   if (MyGlobals::_Verbose>20)
 331     std::cout << "proc " << MyGlobals::_Rank << " : write field " << desc << std::endl;
 332   std::string meshName, fieldName;
 333   int typeField, DT, IT, entity;
 334   FieldShortDescriptionToData(desc, fieldName, typeField, entity, DT, IT);
 335   double time=StrToDouble(ExtractFromDescription(desc, "time="));
 336   int typeData=StrToInt(ExtractFromDescription(desc, "typeData="));
 337   std::string entityName=ExtractFromDescription(desc, "entityName=");
 338   ParaMEDMEM::MEDCouplingFieldDouble* field=0;
 339   if (typeData!=6)
 340     {
 341       std::cout << "WARNING : writeMedFile : typeData " << typeData << " not implemented for fields\n";
 342     }
 343   if (entityName=="MED_CELL")
 344     {
 345       //there is a field of idomain to write
 346       field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME);
 347     }
 348   if (entityName=="MED_NODE_ELEMENT")
 349     {
 350       //there is a field of idomain to write
 351       field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_GAUSS_NE,ParaMEDMEM::ONE_TIME);
 352     }
 353   if (!field)
 354     {
 355       std::cout << "WARNING : writeMedFile : entityName " << entityName << " not implemented for fields\n";
 356     }
 357   if (field && typeData==6)
 358     {
 359       field->setName(fieldName);
 360       field->setMesh(mfm->getGenMeshAtLevel(0));
 361       ParaMEDMEM::DataArrayDouble *da=data;
 362       //get information for components etc..
 363       std::vector<std::string> r1;
 364       r1=SelectTagsInVectorOfString(MyGlobals::_General_Informations,"fieldName="+fieldName);
 365       r1=SelectTagsInVectorOfString(r1,"typeField="+IntToStr(typeField));
 366       r1=SelectTagsInVectorOfString(r1,"DT="+IntToStr(DT));
 367       r1=SelectTagsInVectorOfString(r1,"IT="+IntToStr(IT));
 368       //not saved in file? field->setDescription(ExtractFromDescription(r1[0], "fieldDescription="));
 369       int nbc=StrToInt(ExtractFromDescription(r1[0], "nbComponents="));
 370       if (nbc==da->getNumberOfComponents())
 371         {
 372           for (int i=0; i<nbc; i++)
 373             da->setInfoOnComponent(i,ExtractFromDescription(r1[0], "componentInfo"+IntToStr(i)+"="));
 374         }
 375       else
 376         {
 377           std::cerr << "Problem On field " << fieldName << " : number of components unexpected " << da->getNumberOfComponents() << std::endl;
 378         }
 379       field->setArray(da);
 380       field->setTime(time,DT,IT);
 381       field->checkCoherency();
 382     }
 383   return field;
 384 }
 385
 386 void MeshCollectionDriver::writeMedFile(int idomain, const std::string& distfilename) const
 387 {
 388   ParaMEDMEM::MEDFileMesh* mfm = getMesh( idomain );
 389   mfm->write(distfilename,2);
 390
 391   std::string key="/inewFieldDouble="+IntToStr(idomain)+"/";
 392   std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it;
 393   int nbfFieldFound=0;
 394   for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
 395     {
 396       size_t found=(*it).first.find(key);
 397       if (found==std::string::npos)
 398         continue;
 399       ParaMEDMEM::MEDCouplingFieldDouble* field=0;
 400       field = getField(key, (*it).first, (*it).second, mfm, idomain);
 401       nbfFieldFound++;
 402       try
 403         {
 404           MEDLoader::WriteField(distfilename,field,false);
 405         }
 406       catch(INTERP_KERNEL::Exception& e)
 407         {
 408           //cout trying rewrite all data, only one field defined
 409           std::string tmp,newName=distfilename;
 410           std::string fieldName;
 411           fieldName=field->getName();
 412           tmp+="_"+fieldName+"_"+IntToStr(nbfFieldFound)+".med";
 413           newName.replace(newName.find(".med"),4,tmp);
 414           std::cout << "WARNING : writeMedFile : create a new file name with only one field because MEDLoader::WriteField throw:" << newName << std::endl;
 415           MEDLoader::WriteField(newName,field,true);
 416         }
 417     }
 418   mfm->decrRef();
 419 }
 420
 421 ParaMEDMEM::MEDFileData* MeshCollectionDriver::getMEDFileData()
 422 {
 423   ParaMEDMEM::MEDFileData* newdata = ParaMEDMEM::MEDFileData::New();
 424
 425   ParaMEDMEM::MEDCouplingAutoRefCountObjectPtr<ParaMEDMEM::MEDFileMeshes> meshes;
 426   ParaMEDMEM::MEDCouplingAutoRefCountObjectPtr<ParaMEDMEM::MEDFileFields> fields;
 427   meshes = ParaMEDMEM::MEDFileMeshes::New();
 428   fields = ParaMEDMEM::MEDFileFields::New();
 429
 430   for (size_t i=0; i<_collection->getMesh().size(); i++)
 431     {
 432       ParaMEDMEM::MEDFileMesh* mfm = getMesh( i );
 433       meshes->pushMesh(mfm);
 434
 435       std::string key="/inewFieldDouble="+IntToStr(i)+"/";
 436       std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it;
 437       ParaMEDMEM::MEDFileFieldMultiTS* fieldsMTS = ParaMEDMEM::MEDFileFieldMultiTS::New();
 438       for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
 439         {
 440           size_t found=(*it).first.find(key);
 441           if (found==std::string::npos)
 442             continue;
 443           ParaMEDMEM::MEDCouplingFieldDouble* field=0;
 444           field=getField(key, (*it).first, (*it).second, mfm, i);
 445           ParaMEDMEM::MEDFileField1TS* f1ts = ParaMEDMEM::MEDFileField1TS::New();
 446           f1ts->setFieldNoProfileSBT(field);
 447           fieldsMTS->pushBackTimeStep(f1ts);
 448
 449           field->decrRef();
 450           f1ts->decrRef();
 451         }
 452       fields->pushField(fieldsMTS);
 453
 454       fieldsMTS->decrRef();
 455       mfm->decrRef();
 456     }
 457   newdata->setMeshes(meshes);
 458   newdata->setFields(fields);
 459   return newdata;
 460 }