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