src/MEDPartitioner/MEDPARTITIONER_MeshCollectionDriver.cxx

   1 // Copyright (C) 2007-2013  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
  20 #include "MEDPARTITIONER_ParallelTopology.hxx"
  21 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
  22 #include "MEDPARTITIONER_MeshCollection.hxx"
  23 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
  24 #include "MEDPARTITIONER_Utils.hxx"
  25
  26 #include "MEDCouplingUMesh.hxx"
  27 #include "MEDCouplingFieldDouble.hxx"
  28 #include "MEDLoader.hxx"
  29 #include "MEDFileMesh.hxx"
  30
  31 #include <map>
  32 #include <set>
  33 #include <vector>
  34 #include <string>
  35 #include <fstream>
  36 #include <iostream>
  37
  38 #include <libxml/tree.h>
  39 #include <libxml/parser.h>
  40 #include <libxml/xpath.h>
  41 #include <libxml/xpathInternals.h>
  42
  43 #include "med.h"
  44
  45 using namespace MEDPARTITIONER;
  46
  47 MeshCollectionDriver::MeshCollectionDriver(MeshCollection* collection):_collection(collection)
  48 {
  49 }
  50
  51 /*!reads a unique MED File v>=2.1
  52  * and mounts the corresponding mesh in memory
  53  *\param filename binary file
  54  *\param meshname mesh name in the MED file
  55  * */
  56 int MeshCollectionDriver::readSeq(const char* filename, const char* meshname)
  57 {
  58   std::cout << "readSeq" << std::endl;
  59   MyGlobals::_File_Names.resize(1);
  60   MyGlobals::_File_Names[0]=std::string(filename);
  61
  62   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(filename,meshname);
  63   //puts the only mesh in the mesh vector
  64   (_collection->getMesh()).push_back(mfm->getLevel0Mesh(false));
  65   (_collection->getFaceMesh()).push_back(mfm->getLevelM1Mesh(false));
  66
  67   //reading family ids
  68   ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
  69   ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
  70   (_collection->getCellFamilyIds()).push_back(cellIds);
  71   (_collection->getFaceFamilyIds()).push_back(faceIds);
  72
  73   //reading groups
  74   (_collection->getFamilyInfo())=mfm->getFamilyInfo();
  75   (_collection->getGroupInfo())=mfm->getGroupInfo();
  76
  77   (_collection->getCZ()).clear();
  78
  79   ParallelTopology* aPT = new ParallelTopology((_collection->getMesh()));
  80   _collection->setTopology(aPT);
  81   _collection->setName(meshname);
  82   _collection->setDomainNames(meshname);
  83   return 0;
  84 }
  85
  86
  87 //================================================================================
  88 /*!
  89  * \brief Return mesh dimension from distributed med file had being read
  90  */
  91 //================================================================================
  92
  93 void MeshCollectionDriver::readSubdomain(std::vector<int*>& cellglobal,
  94                                          std::vector<int*>& faceglobal,
  95                                          std::vector<int*>& nodeglobal, int idomain)
  96 {
  97   std::string meshname=MyGlobals::_Mesh_Names[idomain];
  98   std::string file=MyGlobals::_File_Names[idomain];
  99
 100   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file.c_str(),meshname.c_str());
 101   std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
 102
 103   try
 104     {
 105       (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false);
 106       //reading families groups
 107       ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
 108       (_collection->getCellFamilyIds())[idomain]=cellIds;
 109     }
 110   catch(...)
 111     {
 112       (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
 113       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
 114       empty->alloc(0,1);
 115       (_collection->getCellFamilyIds())[idomain]=empty;
 116       std::cout << "\nNO Level0Mesh (Cells)\n";
 117     }
 118   try
 119     {
 120       if (nonEmpty.size()>1 && nonEmpty[1]==-1)
 121         {
 122           (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
 123           //reading families groups
 124           ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
 125           (_collection->getFaceFamilyIds())[idomain]=faceIds;
 126           if (MyGlobals::_Verbose>10)
 127             std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
 128
 129         }
 130       else
 131         {
 132           throw INTERP_KERNEL::Exception("no faces");
 133         }
 134     }
 135   catch(...)
 136     {
 137       (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
 138       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
 139       (_collection->getFaceFamilyIds())[idomain]=empty;
 140       if (MyGlobals::_Verbose>10)
 141         std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
 142     }
 143
 144   //reading groups
 145   _collection->getFamilyInfo()=mfm->getFamilyInfo();
 146   _collection->getGroupInfo()=mfm->getGroupInfo();
 147
 148   mfm->decrRef();
 149
 150   std::vector<std::string> localInformation;
 151   std::string str;
 152   localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
 153   localInformation.push_back(str+"meshName="+meshname);
 154   MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
 155   std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
 156   if (localFields.size()>0)
 157     MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
 158 }
 159
 160
 161 void MeshCollectionDriver::readSubdomain(int idomain)
 162 {
 163   std::string meshname=MyGlobals::_Mesh_Names[idomain];
 164   std::string file=MyGlobals::_File_Names[idomain];
 165
 166   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file.c_str(),meshname.c_str());
 167   std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
 168
 169   try
 170     {
 171       (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false);
 172       //reading families groups
 173       ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
 174       (_collection->getCellFamilyIds())[idomain]=cellIds;
 175     }
 176   catch(...)
 177     {
 178       (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
 179       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
 180       empty->alloc(0,1);
 181       (_collection->getCellFamilyIds())[idomain]=empty;
 182       std::cout<<"\nNO Level0Mesh (Cells)\n";
 183     }
 184   try
 185     {
 186       if (nonEmpty.size()>1 && nonEmpty[1]==-1)
 187         {
 188           (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
 189           //reading families groups
 190           ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
 191           (_collection->getFaceFamilyIds())[idomain]=faceIds;
 192           if (MyGlobals::_Verbose>10)
 193             std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
 194         }
 195       else
 196         {
 197           throw INTERP_KERNEL::Exception("no faces");
 198         }
 199     }
 200   catch(...)
 201     {
 202       (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
 203       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
 204       (_collection->getFaceFamilyIds())[idomain]=empty;
 205       if (MyGlobals::_Verbose>10)
 206         std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
 207     }
 208
 209   //reading groups
 210   _collection->getFamilyInfo()=mfm->getFamilyInfo();
 211   _collection->getGroupInfo()=mfm->getGroupInfo();
 212
 213   mfm->decrRef();
 214
 215   std::vector<std::string> localInformation;
 216   std::string str;
 217   localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
 218   localInformation.push_back(str+"meshName="+meshname);
 219   MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
 220   std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
 221   if (localFields.size()>0)
 222     MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
 223 }
 224
 225
 226 void MeshCollectionDriver::writeMedFile(int idomain, const std::string& distfilename) const
 227 {
 228   std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
 229   ParaMEDMEM::MEDCouplingUMesh* cellMesh=_collection->getMesh(idomain);
 230   ParaMEDMEM::MEDCouplingUMesh* faceMesh=_collection->getFaceMesh(idomain);
 231   //ParaMEDMEM::MEDCouplingUMesh* faceMeshFilter=0;
 232
 233   std::string finalMeshName=ExtractFromDescription(MyGlobals::_General_Informations[0], "finalMeshName=");
 234   // std::string cleFilter=Cle1ToStr("filterFaceOnCell",idomain);
 235   // ParaMEDMEM::DataArrayInt* filter=0;
 236   // if (_collection->getMapDataArrayInt().find(cleFilter)!=_collection->getMapDataArrayInt().end())
 237   //   {
 238   //     filter=_collection->getMapDataArrayInt().find(cleFilter)->second;
 239   //     int* index=filter->getPointer();
 240   //     faceMeshFilter=(ParaMEDMEM::MEDCouplingUMesh *) faceMesh->buildPartOfMySelf(index,index+filter->getNbOfElems(),true);
 241   //     faceMesh=faceMeshFilter;
 242   //   }
 243   cellMesh->setName(finalMeshName.c_str());
 244   meshes.push_back(cellMesh);
 245
 246   faceMesh->checkCoherency();
 247   if (faceMesh->getNumberOfCells()>0)
 248     {
 249       faceMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10);
 250       meshes.push_back(faceMesh);
 251     }
 252
 253   //ParaMEDMEM::MEDCouplingUMesh* boundaryMesh=0;
 254   // if (MyGlobals::_Creates_Boundary_Faces>0)
 255   //   {
 256   //     //try to write Boundary meshes
 257   //     bool keepCoords=false; //TODO or true
 258   //     boundaryMesh=(ParaMEDMEM::MEDCouplingUMesh *) cellMesh->buildBoundaryMesh(keepCoords);
 259   //     boundaryMesh->setName("boundaryMesh");
 260   //   }
 261
 262   MEDLoader::WriteUMeshes(distfilename.c_str(), meshes, true);
 263   // if (faceMeshFilter!=0)
 264   //   faceMeshFilter->decrRef();
 265
 266   // if (boundaryMesh!=0)
 267   //   {
 268   //     //doing that testMesh becomes second mesh sorted by alphabetical order of name
 269   //     MEDLoader::WriteUMesh(distfilename.c_str(), boundaryMesh, false);
 270   //     boundaryMesh->decrRef();
 271   //   }
 272   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(distfilename.c_str(), _collection->getMesh(idomain)->getName());
 273
 274   mfm->setFamilyInfo(_collection->getFamilyInfo());
 275   mfm->setGroupInfo(_collection->getGroupInfo());
 276
 277   std::string key=Cle1ToStr("faceFamily_toArray",idomain);
 278   if ( meshes.size() == 2 &&
 279       _collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
 280     {
 281       ParaMEDMEM::DataArrayInt *fam=_collection->getMapDataArrayInt().find(key)->second;
 282       mfm->setFamilyFieldArr(-1,fam);
 283     }
 284
 285   key=Cle1ToStr("cellFamily_toArray",idomain);
 286   if (_collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
 287     mfm->setFamilyFieldArr(0,_collection->getMapDataArrayInt().find(key)->second);
 288
 289   mfm->write(distfilename.c_str(),0);
 290   key="/inewFieldDouble="+IntToStr(idomain)+"/";
 291
 292   std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it;
 293   int nbfFieldFound=0;
 294   for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
 295     {
 296       std::string desc=(*it).first;
 297       size_t found=desc.find(key);
 298       if (found==std::string::npos)
 299         continue;
 300       if (MyGlobals::_Verbose>20)
 301         std::cout << "proc " << MyGlobals::_Rank << " : write field " << desc << std::endl;
 302       std::string meshName, fieldName;
 303       int typeField, DT, IT, entity;
 304       FieldShortDescriptionToData(desc, fieldName, typeField, entity, DT, IT);
 305       double time=StrToDouble(ExtractFromDescription(desc, "time="));
 306       int typeData=StrToInt(ExtractFromDescription(desc, "typeData="));
 307       std::string entityName=ExtractFromDescription(desc, "entityName=");
 308       ParaMEDMEM::MEDCouplingFieldDouble* field=0;
 309       if (typeData!=6)
 310         {
 311           std::cout << "WARNING : writeMedFile : typeData " << typeData << " not implemented for fields\n";
 312           continue;
 313         }
 314       if (entityName=="MED_CELL")
 315         {
 316           //there is a field of idomain to write
 317           field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME);
 318         }
 319       if (entityName=="MED_NODE_ELEMENT")
 320         {
 321           //there is a field of idomain to write
 322           field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_GAUSS_NE,ParaMEDMEM::ONE_TIME);
 323         }
 324       if (!field)
 325         {
 326           std::cout << "WARNING : writeMedFile : entityName " << entityName << " not implemented for fields\n";
 327           continue;
 328         }
 329       nbfFieldFound++;
 330       field->setName(fieldName.c_str());
 331       field->setMesh(mfm->getLevel0Mesh(false));
 332       ParaMEDMEM::DataArrayDouble *da=(*it).second;
 333
 334       //get information for components etc..
 335       std::vector<std::string> r1;
 336       r1=SelectTagsInVectorOfString(MyGlobals::_General_Informations,"fieldName="+fieldName);
 337       r1=SelectTagsInVectorOfString(r1,"typeField="+IntToStr(typeField));
 338       r1=SelectTagsInVectorOfString(r1,"DT="+IntToStr(DT));
 339       r1=SelectTagsInVectorOfString(r1,"IT="+IntToStr(IT));
 340       //not saved in file? field->setDescription(ExtractFromDescription(r1[0], "fieldDescription=").c_str());
 341       int nbc=StrToInt(ExtractFromDescription(r1[0], "nbComponents="));
 342       if (nbc==da->getNumberOfComponents())
 343         {
 344           for (int i=0; i<nbc; i++)
 345             da->setInfoOnComponent(i,ExtractFromDescription(r1[0], "componentInfo"+IntToStr(i)+"=").c_str());
 346         }
 347       else
 348         {
 349           std::cerr << "Problem On field " << fieldName << " : number of components unexpected " << da->getNumberOfComponents() << std::endl;
 350         }
 351
 352       field->setArray(da);
 353       field->setTime(time,DT,IT);
 354       field->checkCoherency();
 355       try
 356         {
 357           MEDLoader::WriteField(distfilename.c_str(),field,false);
 358         }
 359       catch(INTERP_KERNEL::Exception& e)
 360         {
 361           //cout trying rewrite all data, only one field defined
 362           std::string tmp,newName=distfilename;
 363           tmp+="_"+fieldName+"_"+IntToStr(nbfFieldFound)+".med";
 364           newName.replace(newName.find(".med"),4,tmp);
 365           std::cout << "WARNING : writeMedFile : create a new file name with only one field because MEDLoader::WriteField throw:" << newName << std::endl;
 366           MEDLoader::WriteField(newName.c_str(),field,true);
 367         }
 368     }
 369   mfm->decrRef();
 370 }