src/MEDPartitioner/MEDPARTITIONER_MeshCollection.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_MeshCollection.hxx"
  21
  22 #include "MEDPARTITIONER_ConnectZone.hxx"
  23 #include "MEDPARTITIONER_Graph.hxx"
  24 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
  25 #include "MEDPARTITIONER_MeshCollectionMedAsciiDriver.hxx"
  26 #include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx"
  27 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
  28 #include "MEDPARTITIONER_ParallelTopology.hxx"
  29 #include "MEDPARTITIONER_Topology.hxx"
  30 #include "MEDPARTITIONER_UserGraph.hxx"
  31 #include "MEDPARTITIONER_Utils.hxx"
  32
  33 #ifdef HAVE_MPI
  34 #include "MEDPARTITIONER_JointFinder.hxx"
  35 #endif
  36
  37 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
  38 #include "MEDCouplingFieldDouble.hxx"
  39 #include "MEDCouplingMemArray.hxx"
  40 #include "MEDCouplingNormalizedUnstructuredMesh.hxx"
  41 #include "MEDCouplingUMesh.hxx"
  42 #include "MEDLoader.hxx"
  43 #include "MEDLoaderBase.hxx"
  44
  45 #ifdef HAVE_MPI
  46 #include <mpi.h>
  47 #endif
  48
  49 #ifdef MED_ENABLE_PARMETIS
  50 #include "MEDPARTITIONER_ParMetisGraph.hxx"
  51 #endif
  52 #ifdef MED_ENABLE_METIS
  53 #include "MEDPARTITIONER_MetisGraph.hxx"
  54 #endif
  55 #ifdef MED_ENABLE_SCOTCH
  56 #include "MEDPARTITIONER_ScotchGraph.hxx"
  57 #endif
  58
  59 #include <set>
  60 #include <vector>
  61 #include <string>
  62 #include <limits>
  63 #include <iostream>
  64 #include <fstream>
  65
  66 MEDPARTITIONER::MeshCollection::MeshCollection()
  67   : _topology(0),
  68     _owns_topology(false),
  69     _driver(0),
  70     _domain_selector( 0 ),
  71     _i_non_empty_mesh(-1),
  72     _driver_type(MEDPARTITIONER::MedXml),
  73     _subdomain_boundary_creates( MyGlobals::_Creates_Boundary_Faces ),
  74     _family_splitting(false),
  75     _create_empty_groups(false),
  76     _joint_finder(0)
  77 {
  78 }
  79
  80 /*!constructor creating a new mesh collection (mesh series + topology)
  81  *from an old collection and a new topology
  82  *
  83  * On output, the constructor has built the meshes corresponding to the new mesh collection.
  84  * The new topology has been updated so that face and node mappings are included.
  85  * The families have been cast to their projections in the new topology.
  86  *
  87  * \param initial_collection collection from which the data (coordinates, connectivity) are taken
  88  * \param topology topology containing the cell mappings
  89  */
  90
  91 MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection,
  92                                                Topology* topology,
  93                                                bool family_splitting,
  94                                                bool create_empty_groups)
  95   : _topology(topology),
  96     _owns_topology(false),
  97     _driver(0),
  98     _domain_selector( initialCollection._domain_selector ),
  99     _i_non_empty_mesh(-1),
 100     _name(initialCollection._name),
 101     _driver_type(MEDPARTITIONER::MedXml),
 102     _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
 103     _family_splitting(family_splitting),
 104     _create_empty_groups(create_empty_groups),
 105     _joint_finder(0)
 106 {
 107   std::vector<std::vector<std::vector<int> > > new2oldIds(initialCollection.getTopology()->nbDomain());
 108   std::vector<ParaMEDMEM::DataArrayInt*> o2nRenumber;
 109
 110   castCellMeshes(initialCollection, new2oldIds, o2nRenumber );
 111
 112   //defining the name for the collection and the underlying meshes
 113   setName(initialCollection.getName());
 114
 115   /////////////////
 116   //treating faces
 117   /////////////////
 118
 119 #ifdef HAVE_MPI
 120   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
 121     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
 122 #endif
 123   if (MyGlobals::_Is0verbose)
 124     std::cout<<"treating faces"<<std::endl;
 125   NodeMapping nodeMapping;
 126   //nodeMapping contains the mapping between old nodes and new nodes
 127   // (iolddomain,ioldnode)->(inewdomain,inewnode)
 128   createNodeMapping(initialCollection, nodeMapping);
 129   std::vector<std::vector<std::vector<int> > > new2oldFaceIds;
 130   castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds);
 131
 132   ////////////////////
 133   //treating families
 134   ////////////////////
 135 #ifdef HAVE_MPI
 136   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
 137     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
 138 #endif
 139   if (MyGlobals::_Is0verbose)
 140     {
 141       if (isParallelMode())
 142         std::cout << "ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
 143       else
 144         std::cout << "NOT ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
 145     }
 146   if (MyGlobals::_Is0verbose>10)
 147     std::cout<<"treating cell and face families"<<std::endl;
 148
 149   castIntField(initialCollection.getMesh(),
 150                this->getMesh(),
 151                initialCollection.getCellFamilyIds(),
 152                "cellFamily");
 153   castIntField(initialCollection.getFaceMesh(),
 154                this->getFaceMesh(),
 155                initialCollection.getFaceFamilyIds(),
 156                "faceFamily");
 157
 158   //treating groups
 159 #ifdef HAVE_MPI
 160   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
 161     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
 162 #endif
 163   if (MyGlobals::_Is0verbose)
 164     std::cout << "treating groups" << std::endl;
 165   _family_info=initialCollection.getFamilyInfo();
 166   _group_info=initialCollection.getGroupInfo();
 167
 168 #ifdef HAVE_MPI
 169   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
 170     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
 171 #endif
 172   if (MyGlobals::_Is0verbose)
 173     std::cout << "treating fields" << std::endl;
 174   castAllFields(initialCollection,"cellFieldDouble");
 175   if (_i_non_empty_mesh<0)
 176     {
 177       for (size_t i=0; i<_mesh.size(); i++)
 178         {
 179           if (_mesh[i])
 180             {
 181               _i_non_empty_mesh=i; //first existing one local
 182               break;
 183             }
 184         }
 185     }
 186
 187   // find faces common with neighbor domains and put them in groups
 188   buildBoundaryFaces();
 189
 190   //building the connect zones necessary for writing joints
 191   buildConnectZones( nodeMapping, o2nRenumber, initialCollection.getTopology()->nbDomain() );
 192
 193   // delete o2nRenumber
 194   for ( size_t i = 0; i < o2nRenumber.size(); ++i )
 195     if ( o2nRenumber[i] )
 196       o2nRenumber[i]->decrRef();
 197 }
 198
 199 /*!
 200   Creates the meshes using the topology underlying he mesh collection and the mesh data
 201   coming from the ancient collection
 202   \param initialCollection collection from which the data is extracted to create the new meshes
 203   \param [out] o2nRenumber returns for each new domain a permutation array returned by sortCellsInMEDFileFrmt()
 204 */
 205
 206 void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection,
 207                                                     std::vector<std::vector<std::vector<int> > >& new2oldIds,
 208                                                     std::vector<ParaMEDMEM::DataArrayInt*> & o2nRenumber)
 209 {
 210   if (MyGlobals::_Verbose>10)
 211     std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl;
 212   if (_topology==0)
 213     throw INTERP_KERNEL::Exception("Topology has not been defined on call to castCellMeshes");
 214
 215   int nbNewDomain=_topology->nbDomain();
 216   int nbOldDomain=initialCollection.getTopology()->nbDomain();
 217
 218   _mesh.resize(nbNewDomain);
 219   o2nRenumber.resize(nbNewDomain,0);
 220   int rank=MyGlobals::_Rank;
 221   //splitting the initial domains into smaller bits
 222   std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
 223   splitMeshes.resize(nbNewDomain);
 224   for (int inew=0; inew<nbNewDomain; inew++)
 225     {
 226       splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
 227     }
 228
 229   for (int iold=0; iold<nbOldDomain; iold++)
 230     {
 231       if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
 232         {
 233           int size=(initialCollection._mesh)[iold]->getNumberOfCells();
 234           std::vector<int> globalids(size);
 235           initialCollection.getTopology()->getCellList(iold, &globalids[0]);
 236           std::vector<int> ilocalnew(size); //local
 237           std::vector<int> ipnew(size);     //idomain old
 238           _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
 239
 240           new2oldIds[iold].resize(nbNewDomain);
 241           for (int i=0; i<(int)ilocalnew.size(); i++)
 242             {
 243               new2oldIds[iold][ipnew[i]].push_back(i);
 244             }
 245           for (int inew=0; inew<nbNewDomain; inew++)
 246             {
 247               splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)
 248                 (initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
 249                                                                        &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
 250                                                                        true);
 251               if (MyGlobals::_Verbose>400)
 252                 std::cout<< "proc " << rank << " : a splitMesh iold inew NbCells " << iold << " " << inew << " "
 253                          << splitMeshes[inew][iold]->getNumberOfCells() << std::endl;
 254             }
 255         }
 256     }
 257 #ifdef HAVE_MPI
 258   if (isParallelMode())
 259     {
 260       //if (MyGlobals::_Verbose>300) std::cout<<"proc "<<rank<<" : castCellMeshes send/receive"<<std::endl;
 261       for (int iold=0; iold<nbOldDomain; iold++)
 262         for(int inew=0; inew<nbNewDomain; inew++)
 263           {
 264             if (initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
 265               continue;
 266
 267             if(initialCollection._domain_selector->isMyDomain(iold))
 268               _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew));
 269
 270             if (_domain_selector->isMyDomain(inew))
 271               _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold));
 272
 273           }
 274     }
 275 #endif
 276
 277   //fusing the split meshes
 278   if (MyGlobals::_Verbose>200)
 279     std::cout << "proc " << rank << " : castCellMeshes fusing" << std::endl;
 280   for (int inew=0; inew<nbNewDomain ;inew++)
 281     {
 282       std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
 283
 284       for (int i=0; i<(int)splitMeshes[inew].size(); i++)
 285         if (splitMeshes[inew][i]!=0)
 286           if (splitMeshes[inew][i]->getNumberOfCells()>0)
 287             meshes.push_back(splitMeshes[inew][i]);
 288
 289       if (!isParallelMode()||_domain_selector->isMyDomain(inew))
 290         {
 291           if (meshes.size()==0)
 292             {
 293               _mesh[inew]=CreateEmptyMEDCouplingUMesh();
 294               std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl;
 295             }
 296           else
 297             {
 298               _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes);
 299               o2nRenumber[inew]=_mesh[inew]->sortCellsInMEDFileFrmt();
 300               bool areNodesMerged;
 301               int nbNodesMerged;
 302               if (meshes.size()>1)
 303                 {
 304                   ParaMEDMEM::DataArrayInt* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
 305                   array->decrRef(); // array is not used in this case
 306                 }
 307               _mesh[inew]->zipCoords();
 308             }
 309         }
 310       for (int i=0;i<(int)splitMeshes[inew].size();i++)
 311         if (splitMeshes[inew][i]!=0)
 312           splitMeshes[inew][i]->decrRef();
 313     }
 314   if (MyGlobals::_Verbose>300)
 315     std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl;
 316 }
 317
 318 /*!
 319   \param initialCollection source mesh collection
 320   \param nodeMapping structure containing the correspondency between nodes in the initial collection and the node(s) in the new collection
 321 */
 322 void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping)
 323 {
 324   using std::vector;
 325   using std::make_pair;
 326   //  NodeMapping reverseNodeMapping;
 327   for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
 328     {
 329
 330       double* bbox;
 331       BBTreeOfDim* tree = 0;
 332       int dim = 3;
 333       if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
 334         {
 335           //      std::map<pair<double,pair<double, double> >, int > nodeClassifier;
 336           ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
 337           double* coordsPtr=coords->getPointer();
 338           dim = coords->getNumberOfComponents();
 339           int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
 340           bbox=new double[nvertices*2*dim];
 341
 342           for (int i=0; i<nvertices*dim;i++)
 343             {
 344               bbox[i*2]=coordsPtr[i]-1e-8;
 345               bbox[i*2+1]=coordsPtr[i]+1e-8;
 346             }
 347           tree=new BBTreeOfDim( dim, bbox,0,0,nvertices,1e-9);
 348         }
 349
 350       for (int inew=0; inew<_topology->nbDomain(); inew++)
 351         {
 352 #ifdef HAVE_MPI
 353           //sending meshes for parallel computation
 354           if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold))
 355             _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold));
 356           else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold))
 357             {
 358               ParaMEDMEM::MEDCouplingUMesh* mesh;
 359               _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
 360               ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords();
 361               for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
 362                 {
 363                   double* coordsPtr=coords->getPointer()+inode*dim;
 364                   vector<int> elems;
 365                   tree->getElementsAroundPoint(coordsPtr,elems);
 366                   if (elems.size()==0) continue;
 367                   nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
 368                 }
 369               mesh->decrRef();
 370             }
 371           else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
 372 #else
 373             if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
 374 #endif
 375               {
 376                 ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords();
 377                 for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
 378                   {
 379                     double* coordsPtr=coords->getPointer()+inode*dim;
 380                     vector<int> elems;
 381                     tree->getElementsAroundPoint(coordsPtr,elems);
 382                     if (elems.size()==0) continue;
 383                     nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
 384                   }
 385               }
 386         }
 387       if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
 388         {
 389           delete tree;
 390           delete[] bbox;
 391         }
 392     }
 393
 394 }
 395
 396 void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
 397 {
 398   using std::vector;
 399   using MEDPARTITIONER::BBTreeOfDim;
 400   if (!&meshOne || !&meshTwo) return;  //empty or not existing
 401   double* bbox;
 402   BBTreeOfDim* tree = 0;
 403   int nv1=meshOne.getNumberOfNodes();
 404   ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords();
 405   int dim = coords->getNumberOfComponents();
 406
 407   bbox=new double[nv1*2*dim];
 408   double* coordsPtr=coords->getPointer();
 409   for (int i=0; i<nv1*dim; i++)
 410     {
 411       bbox[i*2]=coordsPtr[i]-1e-8;
 412       bbox[i*2+1]=coordsPtr[i]+1e-8;
 413     }
 414   tree=new BBTreeOfDim( dim, bbox,0,0,nv1,1e-9);
 415
 416   int nv2=meshTwo.getNumberOfNodes();
 417   nodeIds.resize(nv2,-1);
 418   coords=meshTwo.getCoords();
 419   for (int inode=0; inode<nv2; inode++)
 420     {
 421       double* coordsPtr2=coords->getPointer()+inode*dim;
 422       vector<int> elems;
 423       tree->getElementsAroundPoint(coordsPtr2,elems);
 424       if (elems.size()==0) continue;
 425       nodeIds[inode]=elems[0];
 426     }
 427   delete tree;
 428   delete[] bbox;
 429 }
 430
 431 /*!
 432   creates the face meshes on the new domains from the faces on the old domain and the node mapping
 433   faces at the interface are duplicated
 434 */
 435 void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialCollection,
 436                                                     const std::multimap<std::pair<int,int>, std::pair<int,int> >& nodeMapping,
 437                                                     std::vector<std::vector<std::vector<int> > >& new2oldIds)
 438 {
 439   //splitMeshes structure will contain the partition of
 440   //the old faces on the new ones
 441   //splitMeshes[4][2] contains the faces from old domain 2
 442   //that have to be added to domain 4
 443
 444   using std::vector;
 445   using std::map;
 446   using std::multimap;
 447   using std::pair;
 448   using std::make_pair;
 449
 450   if (MyGlobals::_Verbose>10)
 451     std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl;
 452   if (_topology==0)
 453     throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes");
 454
 455   int nbNewDomain=_topology->nbDomain();
 456   int nbOldDomain=initialCollection.getTopology()->nbDomain();
 457
 458   vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
 459   vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
 460
 461   vector< vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
 462
 463   splitMeshes.resize(nbNewDomain);
 464   for (int inew=0; inew<nbNewDomain; inew++)
 465     {
 466       splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
 467     }
 468   new2oldIds.resize(nbOldDomain);
 469   for (int iold=0; iold<nbOldDomain; iold++) new2oldIds[iold].resize(nbNewDomain);
 470
 471   //init null pointer for empty meshes
 472   for (int inew=0; inew<nbNewDomain; inew++)
 473     {
 474       for (int iold=0; iold<nbOldDomain; iold++)
 475         {
 476           splitMeshes[inew][iold]=0; //null for empty meshes
 477           new2oldIds[iold][inew].clear();
 478         }
 479     }
 480
 481   //loop over the old domains to analyse the faces and decide
 482   //on which new domain they belong
 483   for (int iold=0; iold<nbOldDomain; iold++)
 484     {
 485       if (isParallelMode() && !initialCollection._domain_selector->isMyDomain(iold)) continue;
 486       if (MyGlobals::_Verbose>400)
 487         std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes iodDomain "<<iold<<std::endl;
 488       //initial face mesh known : in my domain
 489       if (meshesCastFrom[iold] != 0)
 490         {
 491           for (int ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
 492             {
 493               vector<int> nodes;
 494               meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
 495
 496               map <int,int> faces;
 497
 498               //analysis of element ielem
 499               //counters are set for the element
 500               //for each source node, the mapping is interrogated and the domain counters
 501               //are incremented for each target node
 502               //the face is considered as going to target domains if the counter of the domain
 503               //is equal to the number of nodes
 504               for (int inode=0; inode<(int)nodes.size(); inode++)
 505                 {
 506                   typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
 507                   int mynode=nodes[inode];
 508
 509                   pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
 510                   for (MI iter=myRange.first; iter!=myRange.second; iter++)
 511                     {
 512                       int inew=iter->second.first;
 513                       if (faces.find(inew)==faces.end())
 514                         faces[inew]=1;
 515                       else
 516                         faces[inew]++;
 517                     }
 518                 }
 519
 520               for (map<int,int>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
 521                 {
 522                   if (iter->second==(int)nodes.size())
 523                     //cvw eligible but may be have to be face of a cell of this->getMesh()[inew]?
 524                     //it is not sure here...
 525                     //done before writeMedfile on option?... see filterFaceOnCell()
 526                     new2oldIds[iold][iter->first].push_back(ielem);
 527                 }
 528             }
 529
 530           //creating the splitMeshes from the face ids
 531           for (int inew=0; inew<nbNewDomain; inew++)
 532             {
 533               if (meshesCastFrom[iold]->getNumberOfCells() > 0)
 534                 {
 535                   splitMeshes[inew][iold]=
 536                     (ParaMEDMEM::MEDCouplingUMesh*)
 537                     ( meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
 538                                                               &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
 539                                                               true)
 540                       );
 541                   splitMeshes[inew][iold]->zipCoords();
 542                 }
 543               else
 544                 {
 545                   splitMeshes[inew][iold]=CreateEmptyMEDCouplingUMesh();
 546                 }
 547             }
 548         }
 549       else
 550         {
 551           std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes empty mesh from iodDomain "<<iold<<std::endl;
 552         }
 553     }
 554
 555 #ifdef HAVE_MPI
 556   //send/receive stuff
 557   if (isParallelMode())
 558     {
 559       ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
 560       for (int iold=0; iold<nbOldDomain; iold++)
 561         for (int inew=0; inew<nbNewDomain; inew++)
 562           {
 563             if (initialCollection._domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
 564               {
 565                 if (splitMeshes[inew][iold] != 0)
 566                   {
 567                     _domain_selector->sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew));
 568                     //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes send " <<inew<<" "<<iold<<" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
 569                   }
 570                 else
 571                   {
 572                     _domain_selector->sendMesh(*(empty), _domain_selector->getProcessorID(inew));
 573                     //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes sen0 " <<inew<<" "<<iold<<std::endl;
 574                   }
 575               }
 576             if (!initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
 577               _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold));
 578               //int nb=0;
 579               //if (splitMeshes[inew][iold])
 580               //  nb=splitMeshes[inew][iold]->getNumberOfCells();
 581               //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<<inew<<" "<<iold<<" "<<nb<<std::endl;//" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
 582           }
 583       empty->decrRef();
 584     }
 585 #endif
 586
 587   //fusing the split meshes
 588   if (MyGlobals::_Verbose>200)
 589     std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes fusing" << std::endl;
 590   meshesCastTo.resize(nbNewDomain);
 591   for (int inew=0; inew<nbNewDomain; inew++)
 592     {
 593       vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
 594       for (int iold=0; iold<nbOldDomain; iold++)
 595         {
 596           ParaMEDMEM::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
 597           if (umesh!=0)
 598             if (umesh->getNumberOfCells()>0)
 599                 myMeshes.push_back(umesh);
 600         }
 601
 602       ParaMEDMEM::MEDCouplingUMesh *bndMesh = 0;
 603       if ( _subdomain_boundary_creates &&
 604            _mesh[inew] &&
 605            _mesh[inew]->getNumberOfCells()>0 )
 606         {
 607           bndMesh =
 608             ((ParaMEDMEM::MEDCouplingUMesh *)_mesh[inew]->buildBoundaryMesh(/*keepCoords=*/true));
 609           if (bndMesh->getNumberOfCells()>0)
 610             myMeshes.push_back( bndMesh );
 611         }
 612
 613       if (myMeshes.size()>0)
 614         {
 615           meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
 616           meshesCastTo[inew]->sortCellsInMEDFileFrmt()->decrRef();
 617         }
 618       else
 619         {
 620           ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
 621           meshesCastTo[inew]=empty;
 622         }
 623       for (int iold=0; iold<nbOldDomain; iold++)
 624         if (splitMeshes[inew][iold]!=0)
 625           splitMeshes[inew][iold]->decrRef();
 626       if ( bndMesh )
 627         bndMesh->decrRef();
 628     }
 629   if (MyGlobals::_Verbose>300)
 630     std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
 631 }
 632
 633
 634
 635 void MEDPARTITIONER::MeshCollection::castIntField(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
 636                                                   std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
 637                                                   std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
 638                                                   std::string nameArrayTo)
 639 {
 640   using std::vector;
 641
 642   int ioldMax=meshesCastFrom.size();
 643   int inewMax=meshesCastTo.size();
 644
 645
 646   //preparing bounding box trees for accelerating source-target node identifications
 647   if (MyGlobals::_Verbose>99)
 648     std::cout<<"making accelerating structures"<<std::endl;
 649   std::vector<BBTreeOfDim* > acceleratingStructures(ioldMax);
 650   std::vector<ParaMEDMEM::DataArrayDouble*>bbox(ioldMax);
 651   for (int iold =0; iold< ioldMax; iold++)
 652     if (isParallelMode() && _domain_selector->isMyDomain(iold))
 653       {
 654         ParaMEDMEM::DataArrayDouble* sourceCoords=meshesCastFrom[iold]->getBarycenterAndOwner();
 655         bbox[iold]=sourceCoords->computeBBoxPerTuple(1.e-6);
 656         acceleratingStructures[iold]=new BBTreeOfDim( sourceCoords->getNumberOfComponents(), bbox[iold]->getConstPointer(),0,0,bbox[iold]->getNumberOfTuples());
 657         sourceCoords->decrRef();
 658       }
 659   // send-recv operations
 660 #ifdef HAVE_MPI
 661   for (int inew=0; inew<inewMax; inew++)
 662     {
 663       for (int iold=0; iold<ioldMax; iold++)
 664         {
 665           //sending arrays for distant domains
 666           if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
 667             {
 668               //send mesh
 669               _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
 670               //send vector
 671               int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
 672               vector<int>sendIds;
 673               if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
 674               if (size>0) //no empty
 675                 {
 676                   sendIds.resize(size);
 677                   std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
 678                 }
 679               else //empty
 680                 {
 681                   size=0;
 682                   sendIds.resize(size);
 683                 }
 684               SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
 685             }
 686           //receiving arrays from distant domains
 687           if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
 688             {
 689               //receive mesh
 690               vector<int> recvIds;
 691               ParaMEDMEM::MEDCouplingUMesh* recvMesh;
 692               _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
 693               //receive vector
 694               if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
 695               RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
 696               remapIntField(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo,0);
 697               recvMesh->decrRef(); //cww is it correct?
 698             }
 699         }
 700     }
 701 #endif
 702
 703   //local contributions and aggregation
 704   for (int inew=0; inew<inewMax; inew++)
 705     {
 706       for (int iold=0; iold<ioldMax; iold++)
 707         if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
 708           {
 709             remapIntField(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo,acceleratingStructures[iold]);
 710           }
 711     }
 712   for (int iold=0; iold<ioldMax;iold++)
 713     if (isParallelMode() && _domain_selector->isMyDomain(iold))
 714       {
 715         bbox[iold]->decrRef();
 716         delete acceleratingStructures[iold];
 717       }
 718 }
 719
 720 void MEDPARTITIONER::MeshCollection::remapIntField(int inew, int iold,
 721                                                     const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
 722                                                     const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
 723                                                     const int* fromArray,
 724                                                     std::string nameArrayTo,
 725                                                     const BBTreeOfDim* myTree)
 726 {
 727
 728   if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
 729   ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
 730   const double*  tc=targetCoords->getConstPointer();
 731   int targetSize=targetMesh.getNumberOfCells();
 732   int sourceSize=sourceMesh.getNumberOfCells();
 733   if (MyGlobals::_Verbose>200)
 734     std::cout<<"remap vers target de taille "<<targetSize<<std::endl;
 735   std::vector<int> ccI;
 736   std::string str,cle;
 737   str=nameArrayTo+"_toArray";
 738   cle=Cle1ToStr(str,inew);
 739   int* toArray;
 740
 741   const BBTreeOfDim* tree;
 742   bool cleantree=false;
 743   ParaMEDMEM::DataArrayDouble* sourceBBox=0;
 744   int dim = targetCoords->getNumberOfComponents();
 745   if (myTree==0)
 746     {
 747       sourceBBox=sourceMesh.getBarycenterAndOwner()->computeBBoxPerTuple(1e-8);
 748       tree=new BBTreeOfDim( dim, sourceBBox->getConstPointer(),0,0, sourceBBox->getNumberOfTuples(),1e-10);
 749       cleantree=true;
 750     }
 751   else tree=myTree;
 752   //first time iold : create and initiate
 753   if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
 754     {
 755       if (MyGlobals::_Is0verbose>100)
 756         std::cout << "create " << cle << " size " << targetSize << std::endl;
 757       ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
 758       p->alloc(targetSize,1);
 759       p->fillWithZero();
 760       toArray=p->getPointer();
 761       _map_dataarray_int[cle]=p;
 762     }
 763   else //other times iold: refind and complete
 764     {
 765       toArray=_map_dataarray_int.find(cle)->second->getPointer();
 766     }
 767
 768   std::map< int, int > isource2nb; // count coincident elements
 769   std::map<int,int>::iterator i2nb;
 770
 771   for (int itargetnode=0; itargetnode<targetSize; itargetnode++)
 772     {
 773       std::vector<int> intersectingElems;
 774       tree->getElementsAroundPoint(tc+itargetnode*dim,intersectingElems);
 775       if (intersectingElems.size()!=0)
 776         {
 777           int isourcenode=intersectingElems[0];
 778           if ( intersectingElems.size() > 1 )
 779             {
 780               i2nb = isource2nb.insert( std::make_pair( isourcenode, 0 )).first;
 781               isourcenode = intersectingElems[ i2nb->second++ ];
 782             }
 783           if ( isourcenode < sourceSize ) // protection from invalid elements
 784             {
 785               toArray[itargetnode]=fromArray[isourcenode];
 786               ccI.push_back(itargetnode);
 787               ccI.push_back(isourcenode);
 788             }
 789         }
 790     }
 791   if (MyGlobals::_Verbose>200)
 792     std::cout<<"nb points trouves"<<ccI.size()/2<<std::endl;
 793   //memories intersection for future same job on fields (if no existing cle=no intersection)
 794   str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
 795   if (MyGlobals::_Verbose>700)
 796     std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
 797
 798   _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
 799
 800   targetCoords->decrRef();
 801   if (cleantree) delete tree;
 802   if (sourceBBox !=0) sourceBBox->decrRef();
 803 }
 804
 805 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
 806 {
 807   if (nameArrayTo!="cellFieldDouble")
 808     throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
 809
 810   std::string nameTo="typeData=6"; //resume the type of field casted
 811   // send-recv operations
 812   int ioldMax=initialCollection.getMesh().size();
 813   int inewMax=this->getMesh().size();
 814   int iFieldMax=initialCollection.getFieldDescriptions().size();
 815   if (MyGlobals::_Verbose>10)
 816     std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
 817   //see collection.prepareFieldDescriptions()
 818   for (int ifield=0; ifield<iFieldMax; ifield++)
 819     {
 820       std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
 821       if (descriptionField.find(nameTo)==std::string::npos)
 822         continue; //only nameTo accepted in Fields name description
 823 #ifdef HAVE_MPI
 824       for (int inew=0; inew<inewMax; inew++)
 825         {
 826           for (int iold=0; iold<ioldMax; iold++)
 827             {
 828               //sending arrays for distant domains
 829               if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
 830                 {
 831                   int target=_domain_selector->getProcessorID(inew);
 832                   ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
 833                   if (MyGlobals::_Verbose>10)
 834                     std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
 835                   SendDataArrayDouble(field, target);
 836                 }
 837               //receiving arrays from distant domains
 838               if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
 839                 {
 840                   int source=_domain_selector->getProcessorID(iold);
 841                   //receive vector
 842                   if (MyGlobals::_Verbose>10)
 843                     std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
 844                   ParaMEDMEM::DataArrayDouble* field=RecvDataArrayDouble(source);
 845                   remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
 846                 }
 847             }
 848         }
 849 #endif
 850       //local contributions and aggregation
 851       for (int inew=0; inew<inewMax; inew++)
 852         {
 853           for (int iold=0; iold<ioldMax; iold++)
 854             if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
 855               {
 856                 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
 857                 remapDoubleField(inew,iold,field,nameArrayTo,descriptionField);
 858               }
 859         }
 860     }
 861 }
 862
 863 void MEDPARTITIONER::MeshCollection::remapDoubleField(int inew, int iold,
 864                                                        ParaMEDMEM::DataArrayDouble* fromArray,
 865                                                        std::string nameArrayTo,
 866                                                        std::string descriptionField)
 867 //here we use 'cellFamily_ccI inew iold' set in remapIntField
 868 {
 869   if (nameArrayTo!="cellFieldDouble")
 870     throw INTERP_KERNEL::Exception("Error remapDoubleField only on cellFieldDouble");
 871   std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
 872
 873   std::map<std::string,ParaMEDMEM::DataArrayInt*>::iterator it1;
 874   it1=_map_dataarray_int.find(key);
 875   if (it1==_map_dataarray_int.end())
 876     {
 877       std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found" << std::endl;
 878       std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
 879       return;
 880     }
 881   //create ccI in remapIntField
 882   ParaMEDMEM::DataArrayInt *ccI=it1->second;
 883   if (MyGlobals::_Verbose>300)
 884     std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField " << key << " size " << ccI->getNbOfElems() << std::endl;
 885
 886   int nbcell=this->getMesh()[inew]->getNumberOfCells();
 887   int nbcomp=fromArray->getNumberOfComponents();
 888   int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
 889
 890   std::string tag="inewFieldDouble="+IntToStr(inew);
 891   key=descriptionField+SerializeFromString(tag);
 892   int fromArrayNbOfElem=fromArray->getNbOfElems();
 893   int fromArrayNbOfComp=fromArray->getNumberOfComponents();
 894   int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
 895
 896   if (MyGlobals::_Verbose>1000)
 897     {
 898       std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
 899         " fromArray nbOfElems " << fromArrayNbOfElem <<
 900         " nbTuples " << fromArray->getNumberOfTuples() <<
 901         " nbcells " << fromArrayNbOfCell <<
 902         " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
 903     }
 904
 905   ParaMEDMEM::DataArrayDouble* field=0;
 906   std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it2;
 907   it2=_map_dataarray_double.find(key);
 908   if (it2==_map_dataarray_double.end())
 909     {
 910       if (MyGlobals::_Verbose>300)
 911         std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField key '" << key << "' not found and create it" << std::endl;
 912       field=ParaMEDMEM::DataArrayDouble::New();
 913       _map_dataarray_double[key]=field;
 914       field->alloc(nbcell*nbPtGauss,nbcomp);
 915       field->fillWithZero();
 916     }
 917   else
 918     {
 919       field=it2->second;
 920       if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
 921         {
 922           std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
 923             " trying remap of field double on cells : \n" << descriptionField << std::endl;
 924           return;
 925         }
 926     }
 927
 928   if (nbPtGauss==1)
 929     {
 930       field->setPartOfValuesAdv(fromArray,ccI);
 931     }
 932   else
 933     {
 934       //replaced by setPartOfValuesAdv if nbPtGauss==1
 935       int iMax=ccI->getNbOfElems();
 936       int* pccI=ccI->getPointer();
 937       double* pField=field->getPointer();
 938       double* pFrom=fromArray->getPointer();
 939       int itarget, isource, delta=nbPtGauss*nbcomp;
 940       for (int i=0; i<iMax; i=i+2)  //cell
 941         {
 942           itarget=pccI[i];
 943           isource=pccI[i+1];
 944           if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
 945             throw INTERP_KERNEL::Exception("Error field override");
 946           int ita=itarget*delta;
 947           int iso=isource*delta;
 948           for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
 949         }
 950     }
 951 }
 952
 953 namespace
 954 {
 955   using namespace ParaMEDMEM;
 956   //================================================================================
 957   /*!
 958    * \brief Sort correspondence ids of one domain and permute ids of the other accordingly
 959    *  \param [in,out] ids1 - ids of one domain
 960    *  \param [in,out] ids2 - ids of another domain
 961    *  \param [in] delta - a delta to change all ids
 962    *  \param [in] removeEqual - to remove equal ids
 963    *  \return DataArrayInt* - array of ids joined back
 964    */
 965   //================================================================================
 966
 967   DataArrayInt* sortCorrespondences( DataArrayInt* ids1,
 968                                      DataArrayInt* ids2,
 969                                      int           delta,
 970                                      bool removeEqual = false)
 971   {
 972     // sort
 973     MEDCouplingAutoRefCountObjectPtr< DataArrayInt > renumN2O = ids1->buildPermArrPerLevel();
 974     ids1->renumberInPlaceR( renumN2O->begin() );
 975     ids2->renumberInPlaceR( renumN2O->begin() );
 976
 977     if ( removeEqual )
 978       {
 979         ids1 = ids1->buildUnique();
 980         ids2 = ids2->buildUnique();
 981       }
 982     if ( delta != 0 )
 983       {
 984         int * id = ids1->getPointer();
 985         for ( ; id < ids1->end(); ++id )
 986           ++(*id);
 987         id = ids2->getPointer();
 988         for ( ; id < ids2->end(); ++id )
 989           ++(*id);
 990       }
 991
 992     // join
 993     DataArrayInt* ids12 = DataArrayInt::Meld( ids1, ids2 ); // two components
 994     ids12->rearrange( 1 ); // make one component
 995     return ids12;
 996   }
 997
 998   //================================================================================
 999   /*!
1000    * \brief Renumber ids according to mesh->sortCellsInMEDFileFrmt()
1001    *  \param [in,out] ids - cell ids to renumber
1002    *  \param [in] o2nRenumber - renumbering array in "Old to New" mode
1003    */
1004   //================================================================================
1005
1006   void renumber( DataArrayInt* ids, const DataArrayInt* o2nRenumber )
1007   {
1008     if ( !ids || !o2nRenumber )
1009       return;
1010     int *        id = ids->getPointer();
1011     const int * o2n = o2nRenumber->getConstPointer();
1012     for ( ; id < ids->end(); ++id )
1013       {
1014         *id = o2n[ *id ];
1015       }
1016   }
1017 }
1018
1019 //================================================================================
1020 /*!
1021  * \brief Fill up ConnectZone's stored in _topology with nodal correspondences
1022  *  \param [in] nodeMapping - mapping between old nodes and new nodes
1023  *              (iolddomain,ioldnode)->(inewdomain,inewnode)
1024  *  \param [in] o2nRenumber - renumbering array returned by mesh->sortCellsInMEDFileFrmt()
1025  *              per a new domain
1026  *  \param [in] nbInitialDomains - nb of old domains
1027  */
1028 //================================================================================
1029
1030 void MEDPARTITIONER::MeshCollection::buildConnectZones( const NodeMapping& nodeMapping,
1031                                                         const std::vector<ParaMEDMEM::DataArrayInt*> & o2nRenumber,
1032                                                         int                nbInitialDomains)
1033 {
1034   if ( !MyGlobals::_Create_Joints || _topology->nbDomain() < 2 )
1035     return;
1036
1037   if ( MyGlobals::_World_Size > 1 )
1038     {
1039       _topology->getCZ().clear();
1040       return; // not implemented for parallel mode
1041     }
1042
1043   //  at construction, _topology creates cell correspondences basing on Graph information,
1044   // and here we
1045   // 1) add node correspondences,
1046   // 2) split cell correspondences by cell geometry types
1047   // 3) sort ids to be in ascending order
1048
1049   const int nb_domains = _topology->nbDomain();
1050
1051   // ==================================
1052   // 1) add node correspondences
1053   // ==================================
1054
1055   std::vector< std::vector< std::vector< int > > > nodeCorresp( nb_domains );
1056   for ( int idomain = 0; idomain < nb_domains; ++idomain )
1057     {
1058       nodeCorresp[ idomain ].resize( nb_domains );
1059     }
1060
1061   NodeMapping::const_iterator nmIt1, nmIt2 = nodeMapping.begin();
1062   for ( nmIt1 = nmIt2; nmIt1 != nodeMapping.end(); nmIt1 = nmIt2 )
1063     {
1064       // look for an "old" node mapped into several "new" nodes in different domains
1065       int nbSameOld = 0;
1066       while ( ++nmIt2 != nodeMapping.end() && nmIt2->first == nmIt1->first )
1067         nbSameOld += ( nmIt2->second != nmIt1->second );
1068
1069       if ( nbSameOld > 0 )
1070         {
1071           NodeMapping::const_iterator nmEnd = nmIt2;
1072           for ( ; true; ++nmIt1 )
1073             {
1074               nmIt2 = nmIt1;
1075               if ( ++nmIt2 == nmEnd )
1076                 break;
1077               int dom1  = nmIt1->second.first;
1078               int node1 = nmIt1->second.second;
1079               for ( ; nmIt2 != nmEnd; ++nmIt2 )
1080                 {
1081                   int dom2  = nmIt2->second.first;
1082                   int node2 = nmIt2->second.second;
1083                   if ( dom1 != dom2 )
1084                     {
1085                       nodeCorresp[ dom1 ][ dom2 ].push_back( node1 );
1086                       nodeCorresp[ dom1 ][ dom2 ].push_back( node2 );
1087                       nodeCorresp[ dom2 ][ dom1 ].push_back( node2 );
1088                       nodeCorresp[ dom2 ][ dom1 ].push_back( node1 );
1089                     }
1090                 }
1091             }
1092         }
1093     }
1094
1095   // add nodeCorresp to czVec
1096
1097   std::vector<MEDPARTITIONER::ConnectZone*>& czVec = _topology->getCZ();
1098
1099   for ( int idomain = 0; idomain < nb_domains; ++idomain )
1100     {
1101       for ( int idomainNear = 0; idomainNear < nb_domains; ++idomainNear )
1102         {
1103           std::vector< int > & corresp = nodeCorresp[ idomain ][ idomainNear ];
1104           if ( corresp.empty() )
1105             continue;
1106
1107           MEDPARTITIONER::ConnectZone* cz = 0;
1108           for ( size_t i = 0; i < czVec.size() && !cz; ++i )
1109             if ( czVec[i] &&
1110                  czVec[i]->getLocalDomainNumber  () == idomain &&
1111                  czVec[i]->getDistantDomainNumber() == idomainNear )
1112               cz = czVec[i];
1113
1114           if ( !cz )
1115             {
1116               cz = new MEDPARTITIONER::ConnectZone();
1117               cz->setName( "Nodal Connect Zone defined by MEDPARTITIONER" );
1118               cz->setLocalDomainNumber  ( idomain );
1119               cz->setDistantDomainNumber( idomainNear );
1120               czVec.push_back(cz);
1121             }
1122
1123           cz->setNodeCorresp( &corresp[0], corresp.size()/2 );
1124         }
1125     }
1126
1127   // ==========================================================
1128   // 2) split cell correspondences by cell geometry types
1129   // ==========================================================
1130
1131   for ( size_t i = 0; i < czVec.size(); ++i )
1132     {
1133       MEDPARTITIONER::ConnectZone* cz = czVec[i];
1134       if ( !cz                                         ||
1135            cz->getEntityCorrespNumber( 0,0 ) == 0      ||
1136            cz->getLocalDomainNumber  () > (int)_mesh.size() ||
1137            cz->getDistantDomainNumber() > (int)_mesh.size() )
1138         continue;
1139       ParaMEDMEM::MEDCouplingUMesh* mesh1 = _mesh[ cz->getLocalDomainNumber  () ];
1140       ParaMEDMEM::MEDCouplingUMesh* mesh2 = _mesh[ cz->getDistantDomainNumber() ];
1141
1142       // separate ids of two domains
1143       const ParaMEDMEM::MEDCouplingSkyLineArray *corrArray = cz->getEntityCorresp( 0, 0 );
1144       const DataArrayInt* ids12 = corrArray->getValueArray();
1145       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids1, ids2, ids12Sorted;
1146       ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 );
1147       ids2 = ids12->selectByTupleId2( 1, corrArray->getLength(), 2 );
1148
1149       // renumber cells according to mesh->sortCellsInMEDFileFrmt()
1150       renumber( ids1, o2nRenumber[ cz->getLocalDomainNumber() ]);
1151       renumber( ids2, o2nRenumber[ cz->getDistantDomainNumber() ]);
1152
1153       // check nb cell types
1154       std::set<INTERP_KERNEL::NormalizedCellType> types1, types2;
1155       types1 = mesh1->getTypesOfPart( ids1->begin(), ids1->end() );
1156       types2 = mesh2->getTypesOfPart( ids2->begin(), ids2->end() );
1157       if ( types1.size() < 1 || types2.size() < 1 )
1158         continue; // parallel mode?
1159
1160       MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1161       for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1162         if ( czVec[j] &&
1163              czVec[j]->getLocalDomainNumber  () == cz->getDistantDomainNumber() &&
1164              czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1165           cz21 = czVec[j];
1166
1167       if ( types1.size() == 1 && types2.size() == 1 ) // split not needed, only sort
1168         {
1169           ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1170           cz->setEntityCorresp( *types1.begin(), *types2.begin(),
1171                                 ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1172
1173           if ( cz21 )// set 2->1 correspondence
1174           {
1175             ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1176             cz21->setEntityCorresp( *types2.begin(), *types1.begin(),
1177                                     ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1178           }
1179         }
1180       else // split and sort
1181         {
1182           typedef std::pair< std::vector< int >, std::vector< int > > T2Vecs;
1183           T2Vecs idsByType[ INTERP_KERNEL::NORM_MAXTYPE ][ INTERP_KERNEL::NORM_MAXTYPE ];
1184           int t1, t2;
1185
1186           const int nbIds = ids1->getNbOfElems();
1187           const int * p1 = ids1->begin(), * p2 = ids2->begin();
1188           for ( int i = 0; i < nbIds; ++i )
1189             {
1190               t1 = mesh1->getTypeOfCell( p1[ i ]);
1191               t2 = mesh2->getTypeOfCell( p2[ i ]);
1192               T2Vecs & ids = idsByType[ t1 ][ t2 ];
1193               ids.first .push_back( p1[ i ]);
1194               ids.second.push_back( p1[ i ]);
1195             }
1196
1197           const int maxType = int( INTERP_KERNEL::NORM_MAXTYPE );
1198           for ( t1 = 0; t1 < maxType; ++t1 )
1199             for ( t2 = 0; t2 < maxType; ++t2 )
1200               {
1201                 T2Vecs & ids = idsByType[ t1 ][ t2 ];
1202                 if ( ids.first.empty() ) continue;
1203                 p1 = & ids.first[0];
1204                 p2 = & ids.second[0];
1205                 ids1->desallocate();
1206                 ids1->pushBackValsSilent( p1, p1+ids.first.size() );
1207                 ids2->desallocate();
1208                 ids2->pushBackValsSilent( p2, p2+ids.first.size() );
1209                 ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/1 );
1210
1211                 cz->setEntityCorresp( t1, t2,
1212                                       ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1213
1214                 if ( cz21 )// set 2->1 correspondence
1215                   {
1216                     ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0 );
1217                     cz21->setEntityCorresp( t2, t1,
1218                                             ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1219                     break;
1220                   }
1221               }
1222         }// split and sort
1223
1224       cz->setEntityCorresp( 0, 0, 0, 0 ); // erase ids computed by _topology
1225       if ( cz21 )
1226         cz21->setEntityCorresp( 0, 0, 0, 0 );
1227
1228     } // loop on czVec
1229
1230
1231   // ==========================================
1232   // 3) sort node ids to be in ascending order
1233   // ==========================================
1234
1235   const bool removeEqual = ( nbInitialDomains > 1 );
1236
1237   for ( size_t i = 0; i < czVec.size(); ++i )
1238     {
1239       MEDPARTITIONER::ConnectZone* cz = czVec[i];
1240       if ( !cz || cz->getNodeNumber() < 1 )
1241         continue;
1242       if ( cz->getDistantDomainNumber() < cz->getLocalDomainNumber() )
1243         continue; // treat a pair of domains once
1244
1245       MEDPARTITIONER::ConnectZone* cz21 = 0; // zone 2 -> 1
1246       for ( size_t j = 0; j < czVec.size() && !cz21; ++j )
1247         if ( czVec[j] &&
1248              czVec[j]->getLocalDomainNumber  () == cz->getDistantDomainNumber() &&
1249              czVec[j]->getDistantDomainNumber() == cz->getLocalDomainNumber() )
1250           cz21 = czVec[j];
1251
1252       // separate ids of two domains
1253       const ParaMEDMEM::MEDCouplingSkyLineArray *corrArray = cz->getNodeCorresp();
1254       const DataArrayInt *ids12 = corrArray->getValueArray();
1255       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids1, ids2, ids12Sorted;
1256       ids1 = ids12->selectByTupleId2( 0, corrArray->getLength(), 2 );
1257       ids2 = ids12->selectByTupleId2( 1, corrArray->getLength(), 2 );
1258
1259       ids12Sorted = sortCorrespondences( ids1, ids2, /*delta=*/0, removeEqual );
1260       cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1261
1262       if ( cz21 )// set 2->1 correspondence
1263         {
1264           ids12Sorted = sortCorrespondences( ids2, ids1, /*delta=*/0, false );
1265           cz->setNodeCorresp( ids12Sorted->begin(), ids12Sorted->getNbOfElems() / 2 );
1266         }
1267     }
1268 }
1269
1270 //================================================================================
1271 /*!
1272  * \brief Find faces common with neighbor domains and put them in "JOINT_n_p_Faces"
1273  *        group (where "n" and "p" are domain IDs)
1274  */
1275 //================================================================================
1276
1277 void MEDPARTITIONER::MeshCollection::buildBoundaryFaces()
1278 {
1279   if (_topology->nbDomain() < 2 || !_subdomain_boundary_creates )
1280     return;
1281
1282   if ( getMeshDimension() < 2 )
1283     return;
1284
1285   using ParaMEDMEM::MEDCouplingUMesh;
1286   using ParaMEDMEM::DataArrayDouble;
1287   using ParaMEDMEM::DataArrayInt;
1288
1289   std::vector<MEDCouplingUMesh*>& faceMeshes = getFaceMesh();
1290   int nbMeshes = faceMeshes.size();
1291
1292   //preparing bounding box trees for accelerating search of coincident faces
1293   std::vector<BBTreeOfDim* >   bbTrees(nbMeshes);
1294   std::vector<DataArrayDouble*>bbox   (nbMeshes);
1295   for (int inew = 0; inew < nbMeshes-1; inew++)
1296     if ( !isParallelMode() || _domain_selector->isMyDomain(inew) )
1297       {
1298         DataArrayDouble* bcCoords = faceMeshes[inew]->getBarycenterAndOwner();
1299         bbox   [inew] = bcCoords->computeBBoxPerTuple(1.e-6);
1300         bbTrees[inew] = new BBTreeOfDim( bcCoords->getNumberOfComponents(),
1301                                          bbox[inew]->getConstPointer(),0,0,
1302                                          bbox[inew]->getNumberOfTuples());
1303         bcCoords->decrRef();
1304       }
1305
1306   // loop on domains to find joint faces between them
1307   for (int inew1 = 0; inew1 < nbMeshes; inew1++ )
1308     {
1309       for (int inew2 = inew1+1; inew2 < nbMeshes; inew2++ )
1310         {
1311           MEDCouplingUMesh* mesh1 = 0;
1312           MEDCouplingUMesh* mesh2 = 0;
1313           //MEDCouplingUMesh* recvMesh = 0;
1314           bool mesh1Here = true, mesh2Here = true;
1315           if (isParallelMode())
1316             {
1317 #ifdef HAVE_MPI
1318               mesh1Here = _domain_selector->isMyDomain(inew1);
1319               mesh2Here = _domain_selector->isMyDomain(inew2);
1320               if ( !mesh1Here && mesh2Here )
1321                 {
1322                   //send mesh2 to domain of mesh1
1323                   _domain_selector->sendMesh(*faceMeshes[inew2],
1324                                              _domain_selector->getProcessorID(inew1));
1325                 }
1326               else if ( mesh1Here && !mesh2Here )
1327                 {
1328                   //receiving mesh2 from a distant domain
1329                   _domain_selector->recvMesh(mesh2,_domain_selector->getProcessorID(inew2));
1330                   if ( faceMeshes[ inew2 ] )
1331                     faceMeshes[ inew2 ]->decrRef();
1332                   faceMeshes[ inew2 ] = mesh2;
1333                 }
1334 #endif
1335             }
1336           if ( mesh1Here && !mesh1 ) mesh1 = faceMeshes[ inew1 ];
1337           if ( mesh2Here && !mesh2 ) mesh2 = faceMeshes[ inew2 ];
1338
1339           // find coincident faces
1340           std::vector< int > faces1, faces2;
1341           if ( mesh1 && mesh2 )
1342             {
1343               const DataArrayDouble* coords2 = mesh2->getBarycenterAndOwner();
1344               const double*   c2 = coords2->getConstPointer();
1345               const int      dim = coords2->getNumberOfComponents();
1346               const int nbFaces2 = mesh2->getNumberOfCells();
1347               const int nbFaces1 = mesh1->getNumberOfCells();
1348
1349               for (int i2 = 0; i2 < nbFaces2; i2++)
1350                 {
1351                   std::vector<int> coincFaces;
1352                   bbTrees[inew1]->getElementsAroundPoint( c2+i2*dim, coincFaces );
1353                   if (coincFaces.size()!=0)
1354                     {
1355                       int i1 = coincFaces[0];
1356                       // if ( coincFaces.size() > 1 )
1357                       //   {
1358                       //     i2nb = isource2nb.insert( std::make_pair( i1 , 0 )).first;
1359                       //     i1  = coincFaces[ i2nb->second++ ];
1360                       //   }
1361                       if ( i1  < nbFaces1 ) // protection from invalid elements
1362                         {
1363                           faces1.push_back( i1 );
1364                           faces2.push_back( i2 );
1365                         }
1366                     }
1367                 }
1368               coords2->decrRef();
1369             }
1370
1371           if ( isParallelMode())
1372             {
1373 #ifdef HAVE_MPI
1374               if ( mesh1Here && !mesh2Here )
1375                 {
1376                   //send faces2 to domain of recvMesh
1377                   SendIntVec(faces2, _domain_selector->getProcessorID(inew2));
1378                 }
1379               else if ( !mesh1Here && mesh2Here )
1380                 {
1381                   //receiving ids of faces from a domain of mesh1
1382                   RecvIntVec(faces2, _domain_selector->getProcessorID(inew1));
1383                 }
1384 #endif
1385             }
1386           // if ( recvMesh )
1387           //   recvMesh->decrRef();
1388
1389           // Create group "JOINT_inew1_inew2_Faces" and corresponding families
1390           for ( int is2nd = 0; is2nd < 2; ++is2nd )
1391             {
1392               createJointGroup( is2nd ? faces2 : faces1,
1393                                 inew1 , inew2, is2nd );
1394             }
1395
1396         } // loop on the 2nd domains (inew2)
1397     } // loop on the 1st domains (inew1)
1398
1399
1400   // delete bounding box trees
1401   for (int inew = 0; inew < nbMeshes-1; inew++)
1402     if (isParallelMode() && _domain_selector->isMyDomain(inew))
1403       {
1404         bbox[inew]->decrRef();
1405         delete bbTrees[inew];
1406       }
1407 }
1408
1409 //================================================================================
1410 /*!
1411  * \brief Create group "JOINT_inew1_inew2_Faces" and corresponding families
1412  *  \param faces - face ids to include into the group
1413  *  \param inew1 - index of the 1st domain
1414  *  \param inew2 - index of the 2nd domain
1415  *  \param is2nd - in which (1st or 2nd) domain to create the group
1416  */
1417 //================================================================================
1418
1419 void MEDPARTITIONER::MeshCollection::createJointGroup( const std::vector< int >& faces,
1420                                                        const int                 inew1,
1421                                                        const int                 inew2,
1422                                                        const bool                is2nd )
1423 {
1424   // get the name of JOINT group
1425   std::string groupName;
1426   {
1427     std::ostringstream oss;
1428     oss << "JOINT_"
1429         << (is2nd ? inew2 : inew1 ) << "_"
1430         << (is2nd ? inew1 : inew2 ) << "_"
1431         << ( getMeshDimension()==2 ? "Edge" : "Face" );
1432     groupName = oss.str();
1433   }
1434
1435   // remove existing "JOINT_*" group
1436   _group_info.erase( groupName );
1437
1438   // get family IDs array
1439   int* famIDs = 0;
1440   int inew = (is2nd ? inew2 : inew1 );
1441   int totalNbFaces =  _face_mesh[ inew ] ? _face_mesh[ inew ]->getNumberOfCells() : 0;
1442   std::string cle = Cle1ToStr( "faceFamily_toArray", inew );
1443   if ( !_map_dataarray_int.count(cle) )
1444     {
1445       if ( totalNbFaces > 0 )
1446         {
1447           ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
1448           p->alloc( totalNbFaces, 1 );
1449           p->fillWithZero();
1450           famIDs = p->getPointer();
1451           _map_dataarray_int[cle]=p;
1452         }
1453     }
1454   else
1455     {
1456       famIDs = _map_dataarray_int.find(cle)->second->getPointer();
1457     }
1458   // find a family ID of an existing JOINT group
1459   int familyID = 0;
1460   std::map<std::string, int>::iterator name2id = _family_info.find( groupName );
1461   if ( name2id != _family_info.end() )
1462     familyID = name2id->second;
1463
1464   // remove faces from the familyID-the family
1465   if ( familyID != 0 && famIDs )
1466     for ( int i = 0; i < totalNbFaces; ++i )
1467       if ( famIDs[i] == familyID )
1468         famIDs[i] = 0;
1469
1470   if ( faces.empty() )
1471     return;
1472
1473   if ( familyID == 0 )  // generate a family ID for JOINT group
1474     {
1475       std::set< int > familyIDs;
1476       for ( name2id = _family_info.begin(); name2id != _family_info.end(); ++name2id )
1477         familyIDs.insert( name2id->second );
1478       // find the next free family ID
1479       int freeIdCount = inew1 * getNbOfGlobalMeshes() + inew2 + is2nd;
1480       do
1481         {
1482           if ( !familyIDs.count( ++familyID ))
1483             --freeIdCount;
1484         }
1485       while ( freeIdCount > 0 );
1486     }
1487
1488   // push faces to familyID-th group
1489   if ( faces.back() >= totalNbFaces )
1490     throw INTERP_KERNEL::Exception("MeshCollection::createJointGroup(): to high face ID");
1491   for ( size_t i = 0; i < faces.size(); ++i )
1492     famIDs[ faces[i] ] = familyID;
1493
1494   // register JOINT group and family
1495   _family_info[ groupName ] = familyID; // name of the group and family is same
1496   _group_info [ groupName ].push_back( groupName );
1497 }
1498
1499 /*! constructing the MESH collection from a distributed file
1500  *
1501  * \param filename name of the master file containing the list of all the MED files
1502  */
1503 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
1504   : _topology(0),
1505     _owns_topology(true),
1506     _driver(0),
1507     _domain_selector( 0 ),
1508     _i_non_empty_mesh(-1),
1509     _driver_type(MEDPARTITIONER::Undefined),
1510     _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1511     _family_splitting(false),
1512     _create_empty_groups(false),
1513     _joint_finder(0)
1514 {
1515   try
1516     {
1517       _driver=new MeshCollectionMedXmlDriver(this);
1518       _driver->read (filename.c_str());
1519       _driver_type = MedXml;
1520     }
1521   catch(...)
1522     {  // Handle all exceptions
1523       if ( _driver ) delete _driver; _driver=0;
1524       try
1525         {
1526           _driver=new MeshCollectionMedAsciiDriver(this);
1527           _driver->read (filename.c_str());
1528           _driver_type=MedAscii;
1529         }
1530       catch(...)
1531         {
1532           delete _driver;
1533           _driver=0;
1534           throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
1535         }
1536     }
1537   for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1538     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1539       _i_non_empty_mesh = idomain;
1540 }
1541
1542 /*! Constructing the MESH collection from selected domains of a distributed file
1543  *
1544  * \param filename  - name of the master file containing the list of all the MED files
1545  * \param domainSelector - selector of domains to load
1546  */
1547 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
1548   : _topology(0),
1549     _owns_topology(true),
1550     _driver(0),
1551     _domain_selector( &domainSelector ),
1552     _i_non_empty_mesh(-1),
1553     _driver_type(MEDPARTITIONER::Undefined),
1554     _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1555     _family_splitting(false),
1556     _create_empty_groups(false),
1557     _joint_finder(0)
1558 {
1559   std::string myfile=filename;
1560   std::size_t found=myfile.find(".xml");
1561   if (found!=std::string::npos) //file .xml
1562     {
1563       try
1564         {
1565           _driver=new MeshCollectionMedXmlDriver(this);
1566           _driver->read ( filename.c_str(), _domain_selector );
1567           _driver_type = MedXml;
1568         }
1569       catch(...)
1570         {  // Handle all exceptions
1571           delete _driver;
1572           throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
1573         }
1574     }
1575   else
1576     {
1577       found=myfile.find(".med");
1578       if (found!=std::string::npos) //file .med single
1579         {
1580           //make a temporary file .xml and retry MedXmlDriver
1581           std::string xml="\
1582 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
1583 <root>\n \
1584   <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
1585   <description what=\"\" when=\"\"/>\n \
1586   <content>\n \
1587     <mesh name=\"$meshName\"/>\n \
1588   </content>\n \
1589   <splitting>\n \
1590     <subdomain number=\"1\"/>\n \
1591     <global_numbering present=\"no\"/>\n \
1592   </splitting>\n \
1593   <files>\n \
1594     <subfile id=\"1\">\n \
1595       <name>$fileName</name>\n \
1596       <machine>localhost</machine>\n \
1597     </subfile>\n \
1598   </files>\n \
1599   <mapping>\n \
1600     <mesh name=\"$meshName\">\n \
1601       <chunk subdomain=\"1\">\n \
1602         <name>$meshName</name>\n \
1603       </chunk>\n \
1604     </mesh>\n \
1605   </mapping>\n \
1606 </root>\n";
1607           std::vector<std::string> meshNames=MEDLoader::GetMeshNames(myfile);
1608           xml.replace(xml.find("$fileName"),9,myfile);
1609           xml.replace(xml.find("$meshName"),9,meshNames[0]);
1610           xml.replace(xml.find("$meshName"),9,meshNames[0]);
1611           xml.replace(xml.find("$meshName"),9,meshNames[0]);
1612           std::string nameFileXml(myfile);
1613           nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
1614           std::string nameFileXmlDN,nameFileXmlBN;
1615           MEDLoaderBase::getDirAndBaseName(nameFileXml,nameFileXmlDN,nameFileXmlBN);
1616           nameFileXml=MEDLoaderBase::joinPath(nameFileXmlDN,"medpartitioner_"+nameFileXmlBN);
1617           if (_domain_selector->rank()==0) //only on to write it
1618             {
1619               std::ofstream f(nameFileXml.c_str());
1620               f<<xml;
1621               f.close();
1622             }
1623 #ifdef HAVE_MPI
1624            if (MyGlobals::_World_Size>1)
1625              MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
1626 #endif
1627           try
1628             {
1629               _driver=new MeshCollectionMedXmlDriver(this);
1630               _driver->read ( nameFileXml.c_str(), _domain_selector );
1631               _driver_type = MedXml;
1632             }
1633           catch(...)
1634             {  // Handle all exceptions
1635               delete _driver; _driver=0;
1636               throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
1637             }
1638         }
1639       else //no extension
1640         {
1641           try
1642             {
1643               _driver=new MeshCollectionMedAsciiDriver(this);
1644               _driver->read ( filename.c_str(), _domain_selector );
1645               _driver_type=MedAscii;
1646             }
1647           catch(...)
1648             {
1649               delete _driver;
1650               _driver=0;
1651               throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
1652             }
1653         }
1654     }
1655   // find non-empty domain mesh
1656   for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
1657     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
1658       _i_non_empty_mesh = idomain;
1659
1660   try
1661     {
1662       //check for all proc/file compatibility of _field_descriptions
1663 #ifdef HAVE_MPI
1664       _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
1665 #else
1666       _field_descriptions=MyGlobals::_Field_Descriptions;
1667 #endif
1668     }
1669   catch(INTERP_KERNEL::Exception& e)
1670     {
1671       std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1672       throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
1673     }
1674 #ifdef HAVE_MPI
1675   try
1676     {
1677       //check for all proc/file compatibility of _family_info
1678       std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
1679       _family_info=DevectorizeToMapOfStringInt(v2);
1680     }
1681   catch(INTERP_KERNEL::Exception& e)
1682     {
1683       std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1684       throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
1685     }
1686
1687   try
1688     {
1689       //check for all proc/file compatibility of _group_info
1690       std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
1691       _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
1692     }
1693   catch(INTERP_KERNEL::Exception& e)
1694     {
1695       std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
1696       throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
1697     }
1698 #endif
1699 }
1700
1701 /*! constructing the MESH collection from a sequential MED-file
1702  *
1703  * \param filename MED file
1704  * \param meshname name of the mesh that is to be read
1705  */
1706 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
1707   : _topology(0),
1708     _owns_topology(true),
1709     _driver(0),
1710     _domain_selector( 0 ),
1711     _i_non_empty_mesh(-1),
1712     _name(meshname),
1713     _driver_type(MEDPARTITIONER::MedXml),
1714     _subdomain_boundary_creates(MyGlobals::_Creates_Boundary_Faces),
1715     _family_splitting(false),
1716     _create_empty_groups(false),
1717     _joint_finder(0)
1718 {
1719   try // avoid memory leak in case of inexistent filename
1720     {
1721       retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
1722     }
1723   catch (...)
1724     {
1725       delete _driver;
1726       _driver=0;
1727       throw INTERP_KERNEL::Exception("problem reading .med files");
1728     }
1729   if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
1730     _i_non_empty_mesh = 0;
1731 }
1732
1733 MEDPARTITIONER::MeshCollection::~MeshCollection()
1734 {
1735   for (int i=0; i<(int)_mesh.size();i++)
1736     if (_mesh[i]!=0) _mesh[i]->decrRef();
1737
1738   for (int i=0; i<(int)_cell_family_ids.size();i++)
1739     if (_cell_family_ids[i]!=0)
1740       _cell_family_ids[i]->decrRef();
1741
1742   for (int i=0; i<(int)_face_mesh.size();i++)
1743     if (_face_mesh[i]!=0)
1744       _face_mesh[i]->decrRef();
1745
1746   for (int i=0; i<(int)_face_family_ids.size();i++)
1747     if (_face_family_ids[i]!=0)
1748       _face_family_ids[i]->decrRef();
1749
1750   for (std::map<std::string, ParaMEDMEM::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
1751     if ((*it).second!=0)
1752       (*it).second->decrRef();
1753
1754   for (std::map<std::string, ParaMEDMEM::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
1755     if ((*it).second!=0)
1756       (*it).second->decrRef();
1757
1758   delete _driver;
1759   if (_topology!=0 && _owns_topology)
1760     delete _topology;
1761 #ifdef HAVE_MPI
1762   delete _joint_finder;
1763 #endif
1764 }
1765
1766 /*! constructing the MESH collection from a file
1767  *
1768  * The method creates as many MED-files as there are domains in the
1769  * collection. It also creates a master file that lists all the MED files.
1770  * The MED files created in ths manner contain joints that describe the
1771  * connectivity between subdomains.
1772  *
1773  * \param filename name of the master file that will contain the list of the MED files
1774  *
1775  */
1776 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
1777 {
1778   //suppresses link with driver so that it can be changed for writing
1779   delete _driver;
1780   _driver=0;
1781   retrieveDriver()->write (filename.c_str(), _domain_selector);
1782 }
1783
1784 /*! creates or gets the link to the collection driver
1785  */
1786 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
1787 {
1788   if (_driver==0)
1789     {
1790       switch(_driver_type)
1791         {
1792         case MedXml:
1793           _driver=new MeshCollectionMedXmlDriver(this);
1794           break;
1795         case MedAscii:
1796           _driver=new MeshCollectionMedAsciiDriver(this);
1797           break;
1798         default:
1799           throw INTERP_KERNEL::Exception("Unrecognized driver");
1800         }
1801     }
1802   return _driver;
1803 }
1804
1805
1806 /*! gets an existing driver
1807  *
1808  */
1809 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
1810 {
1811   return _driver;
1812 }
1813
1814 /*!
1815  * retrieves the mesh dimension
1816 */
1817 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
1818 {
1819   return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
1820 }
1821
1822 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
1823 {
1824   int nb=0;
1825   for (size_t i=0; i<_mesh.size(); i++)
1826     {
1827       if (_mesh[i]) nb++;
1828     }
1829   return nb;
1830 }
1831
1832 int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
1833 {
1834   int nb=0;
1835   for (size_t i=0; i<_mesh.size(); i++)
1836     {
1837       if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
1838     }
1839   return nb;
1840 }
1841
1842 int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
1843 {
1844   int nb=0;
1845   for (size_t i=0; i<_face_mesh.size(); i++)
1846     {
1847       if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
1848     }
1849   return nb;
1850 }
1851
1852 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
1853 {
1854   return _mesh;
1855 }
1856
1857 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
1858 {
1859   return _face_mesh;
1860 }
1861
1862 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
1863 {
1864   return _mesh[idomain];
1865 }
1866
1867 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
1868 {
1869   return _face_mesh[idomain];
1870 }
1871
1872 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
1873 {
1874   if ( _topology )
1875     return _topology->getCZ();
1876
1877   static std::vector<MEDPARTITIONER::ConnectZone*> noCZ;
1878   return noCZ;
1879 }
1880
1881 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
1882 {
1883   return _topology;
1884 }
1885
1886 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo, bool takeOwneship)
1887 {
1888   if (_topology!=0)
1889     {
1890       throw INTERP_KERNEL::Exception("topology is already set");
1891     }
1892   else
1893     {
1894       _topology = topo;
1895       _owns_topology = takeOwneship;
1896     }
1897 }
1898
1899 /*! Method creating the cell graph in serial mode
1900  *
1901  * \param array returns the pointer to the structure that contains the graph
1902  * \param edgeweight returns the pointer to the table that contains the edgeweights
1903  *        (only used if indivisible regions are required)
1904  */
1905 void MEDPARTITIONER::MeshCollection::buildCellGraph(ParaMEDMEM::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1906 {
1907
1908   using std::map;
1909   using std::vector;
1910   using std::make_pair;
1911   using std::pair;
1912
1913   if (_topology->nbDomain()>1) throw INTERP_KERNEL::Exception("buildCellGraph should be used for one domain only");
1914   const ParaMEDMEM::MEDCouplingUMesh* mesh=_mesh[0];
1915   if (MyGlobals::_Verbose>50)
1916     std::cout<<"getting nodal connectivity"<<std::endl;
1917
1918   //looking for reverse nodal connectivity i global numbering
1919   if (isParallelMode() && !_domain_selector->isMyDomain(0))
1920      {
1921         vector<int> value;
1922         vector<int> index(1,0);
1923
1924         array=new ParaMEDMEM::MEDCouplingSkyLineArray(index,value);
1925         return;
1926      }
1927   array=mesh->generateGraph();
1928 }
1929 /*! Method creating the cell graph in multidomain mode
1930  *
1931  * \param array returns the pointer to the structure that contains the graph
1932  * \param edgeweight returns the pointer to the table that contains the edgeweights
1933  *        (only used if indivisible regions are required)
1934  */
1935 void MEDPARTITIONER::MeshCollection::buildParallelCellGraph(ParaMEDMEM::MEDCouplingSkyLineArray* & array, int *& edgeweights )
1936 {
1937   using std::multimap;
1938   using std::vector;
1939   using std::make_pair;
1940   using std::pair;
1941
1942   std::multimap< int, int > node2cell;
1943   std::map< pair<int,int>, int > cell2cellcounter;
1944   std::multimap<int,int> cell2cell;
1945
1946   std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
1947   int nbdomain=_topology->nbDomain();
1948 #ifdef HAVE_MPI
1949   if (isParallelMode())
1950     {
1951       _joint_finder=new JointFinder(*this);
1952       _joint_finder->findCommonDistantNodes();
1953       commonDistantNodes=_joint_finder->getDistantNodeCell();
1954     }
1955
1956   if (MyGlobals::_Verbose>500)
1957     _joint_finder->print();
1958 #endif
1959
1960   if (MyGlobals::_Verbose>50)
1961     std::cout<<"getting nodal connectivity"<<std::endl;
1962   //looking for reverse nodal connectivity i global numbering
1963   int meshDim = 3;
1964   for (int idomain=0; idomain<nbdomain; idomain++)
1965     {
1966       if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
1967         continue;
1968       meshDim = _mesh[idomain]->getMeshDimension();
1969
1970       ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
1971       ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
1972       int nbNodes=_mesh[idomain]->getNumberOfNodes();
1973       _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
1974       //problem saturation over 1 000 000 nodes for 1 proc
1975       if (MyGlobals::_Verbose>100)
1976         std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
1977       int* index_ptr=index->getPointer();
1978       int* revConnPtr=revConn->getPointer();
1979       for (int i=0; i<nbNodes; i++)
1980         {
1981           for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
1982             {
1983               int globalNode=_topology->convertNodeToGlobal(idomain,i);
1984               int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
1985               node2cell.insert(make_pair(globalNode, globalCell));
1986             }
1987         }
1988       revConn->decrRef();
1989       index->decrRef();
1990 #ifdef HAVE_MPI
1991       for (int iother=0; iother<nbdomain; iother++)
1992         {
1993           std::multimap<int,int>::iterator it;
1994           int isource=idomain;
1995           int itarget=iother;
1996           for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin();
1997                it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
1998             {
1999               int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
2000               int globalCell=(*it).second;
2001               node2cell.insert(make_pair(globalNode, globalCell));
2002             }
2003         }
2004 #endif
2005     }  //endfor idomain
2006
2007   //creating graph arcs (cell to cell relations)
2008   //arcs are stored in terms of (index,value) notation
2009   // 0 3 5 6 6
2010   // 1 2 3 2 3 3
2011   // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
2012   // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
2013
2014   //warning here one node have less than or equal effective number of cell with it
2015   //but cell could have more than effective nodes
2016   //because other equals nodes in other domain (with other global inode)
2017   if (MyGlobals::_Verbose>50)
2018     std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
2019
2020   for (int inode=0;inode<_topology->nbNodes();inode++)
2021     {
2022       typedef multimap<int,int>::const_iterator MI;
2023       std::pair <MI,MI> nodeRange=node2cell.equal_range(inode);
2024       for (MI cell1=nodeRange.first;cell1!=nodeRange.second;cell1++)
2025         for (MI cell2=nodeRange.first;cell2!=cell1;cell2++)
2026           {
2027             int icell1=cell1->second;
2028             int icell2=cell2->second;
2029             if (icell1>icell2) {int tmp=icell1; icell1=icell2; icell2=tmp;}
2030             std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2031             if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2032             else (it->second)++;
2033           }
2034     }
2035   // for (int icell1=0; icell1<_topology->nbCells(); icell1++)  //on all nodes
2036   //   {
2037   //     typedef multimap<int,int>::const_iterator MI;
2038   //     std::pair <MI,MI> nodeRange=cell2node.equal_range(icell1);
2039   //     for (MI node1=nodeRange.first; node1!=nodeRange.second; node1++)  //on nodes with icell
2040   //       {
2041   //         std::pair<MI,MI> cellRange=node2cell.equal_range(node1->second);
2042   //         for (MI cell2=cellRange.first; cell2!=cellRange.second; cell2++)  //on one of these cell
2043   //           {
2044   //             int icell2=cell2->second;
2045   //             std::map<pair<int,int>,int>::iterator it=cell2cellcounter.find(make_pair(icell1,icell2));
2046   //             if (it==cell2cellcounter.end()) cell2cellcounter.insert(make_pair(make_pair(icell1,icell2),1));
2047   //             else (it->second)++;
2048   //           }
2049   //       }
2050   //   }
2051
2052
2053   //converting the counter to a multimap structure
2054   for (std::map<pair<int,int>,int>::const_iterator it=cell2cellcounter.begin();
2055        it!=cell2cellcounter.end();
2056        it++)
2057     if (it->second>=meshDim)
2058       {
2059         cell2cell.insert(std::make_pair(it->first.first,it->first.second));
2060         cell2cell.insert(std::make_pair(it->first.second, it->first.first));
2061       }
2062
2063
2064   if (MyGlobals::_Verbose>50)
2065     std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
2066   //filling up index and value to create skylinearray structure
2067   std::vector <int> index,value;
2068   index.push_back(0);
2069   int idep=0;
2070
2071   for (int idomain=0; idomain<nbdomain; idomain++)
2072     {
2073       if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
2074       int nbCells=_mesh[idomain]->getNumberOfCells();
2075       for (int icell=0; icell<nbCells; icell++)
2076         {
2077           int size=0;
2078           int globalCell=_topology->convertCellToGlobal(idomain,icell);
2079           multimap<int,int>::iterator it;
2080           pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
2081           ret=cell2cell.equal_range(globalCell);
2082           for (it=ret.first; it!=ret.second; ++it)
2083             {
2084               int ival=(*it).second; //no adding one existing yet
2085               for (int i=idep ; i<idep+size ; i++)
2086                 {
2087                   if (value[i]==ival)
2088                     {
2089                       ival= -1; break;
2090                     }
2091                 }
2092               if (ival!= -1)
2093                 {
2094                   value.push_back(ival);
2095                   size++;
2096                 }
2097             }
2098           idep=index[index.size()-1]+size;
2099           index.push_back(idep);
2100         }
2101     }
2102
2103   array=new ParaMEDMEM::MEDCouplingSkyLineArray(index,value);
2104
2105   if (MyGlobals::_Verbose>100)
2106     {
2107       std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
2108         index.size()-1 << " " << value.size() << std::endl;
2109       int max=index.size()>15?15:index.size();
2110       if (index.size()>1)
2111         {
2112           for (int i=0; i<max; ++i)
2113             std::cout<<index[i]<<" ";
2114           std::cout << "... " << index[index.size()-1] << std::endl;
2115           for (int i=0; i<max; ++i)
2116             std::cout<< value[i] << " ";
2117           int ll=index[index.size()-1]-1;
2118           std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
2119         }
2120     }
2121
2122 }
2123
2124
2125 /*! Creates the partition corresponding to the cell graph and the partition number
2126  *
2127  * \param nbdomain number of subdomains for the newly created graph
2128  *
2129  * returns a topology based on the new graph
2130  */
2131 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
2132                                                                           Graph::splitter_type split,
2133                                                                           const std::string& options_string,
2134                                                                           int *user_edge_weights,
2135                                                                           int *user_vertices_weights)
2136 {
2137   if (MyGlobals::_Verbose>10)
2138     std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
2139
2140   if (nbdomain <1)
2141     throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
2142   ParaMEDMEM::MEDCouplingSkyLineArray* array=0;
2143   int* edgeweights=0;
2144
2145   if (_topology->nbDomain()>1 || isParallelMode())
2146     buildParallelCellGraph(array,edgeweights);
2147   else
2148     buildCellGraph(array,edgeweights);
2149
2150   Graph* cellGraph = 0;
2151   switch (split)
2152     {
2153     case Graph::METIS:
2154       if ( isParallelMode() && MyGlobals::_World_Size > 1 )
2155         {
2156 #ifdef MED_ENABLE_PARMETIS
2157           if (MyGlobals::_Verbose>10)
2158             std::cout << "ParMETISGraph" << std::endl;
2159           cellGraph=new ParMETISGraph(array,edgeweights);
2160 #endif
2161         }
2162       if ( !cellGraph )
2163         {
2164 #ifdef MED_ENABLE_METIS
2165           if (MyGlobals::_Verbose>10)
2166             std::cout << "METISGraph" << std::endl;
2167           cellGraph=new METISGraph(array,edgeweights);
2168 #endif
2169         }
2170       if ( !cellGraph )
2171         throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
2172       break;
2173
2174     case Graph::SCOTCH:
2175 #ifdef MED_ENABLE_SCOTCH
2176       if (MyGlobals::_Verbose>10)
2177         std::cout << "SCOTCHGraph" << std::endl;
2178       cellGraph=new SCOTCHGraph(array,edgeweights);
2179 #else
2180       throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
2181 #endif
2182       break;
2183     }
2184
2185   //!user-defined weights
2186   if (user_edge_weights!=0)
2187     cellGraph->setEdgesWeights(user_edge_weights);
2188   if (user_vertices_weights!=0)
2189     cellGraph->setVerticesWeights(user_vertices_weights);
2190
2191   if (MyGlobals::_Is0verbose>10)
2192     std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
2193   cellGraph->partGraph(nbdomain, options_string, _domain_selector);
2194
2195   if (MyGlobals::_Is0verbose>10)
2196     std::cout << "building new topology" << std::endl;
2197   //cellGraph is a shared pointer
2198   Topology *topology=0;
2199   topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2200   //cleaning
2201   delete [] edgeweights;
2202   delete cellGraph;
2203   if (MyGlobals::_Verbose>11)
2204     std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
2205   return topology;
2206 }
2207
2208 /*! Creates a topology for a partition specified by the user
2209  *
2210  * \param table user-specified partition (for each cell contains the domain number from 0 to n-1)
2211  *
2212  * returns a topology based on the new partition
2213  */
2214 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
2215 {
2216   ParaMEDMEM::MEDCouplingSkyLineArray* array=0;
2217   int* edgeweights=0;
2218
2219   if ( _topology->nbDomain()>1)
2220     buildParallelCellGraph(array,edgeweights);
2221   else
2222     buildCellGraph(array,edgeweights);
2223
2224   Graph* cellGraph;
2225   std::set<int> domains;
2226   for (int i=0; i<_topology->nbCells(); i++)
2227     {
2228       domains.insert(partition[i]);
2229     }
2230   cellGraph=new UserGraph(array, partition, _topology->nbCells());
2231
2232   //cellGraph is a shared pointer
2233   Topology *topology=0;
2234   int nbdomain=domains.size();
2235   topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
2236   //  if (array!=0) delete array;
2237   delete cellGraph;
2238   return topology;
2239 }
2240
2241 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
2242 {
2243   for (int i=0; i<_topology->nbDomain(); i++)
2244     {
2245       std::ostringstream oss;
2246       oss<<name<<"_"<<i;
2247       if (!isParallelMode() || _domain_selector->isMyDomain(i))
2248         _mesh[i]->setName(oss.str());
2249     }
2250 }
2251
2252 ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
2253 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
2254 //something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name,f1->getMesh()->getName(),0,f1->getName(),0,1);
2255 {
2256   int rank=MyGlobals::_Rank;
2257   std::string tag="ioldFieldDouble="+IntToStr(iold);
2258   std::string descriptionIold=descriptionField+SerializeFromString(tag);
2259   if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
2260     {
2261       if (MyGlobals::_Verbose>300)
2262         std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
2263       ParaMEDMEM::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
2264       return res;
2265     }
2266   if (MyGlobals::_Verbose>200)
2267     std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
2268   std::string description, fileName, meshName, fieldName;
2269   int typeField, DT, IT, entity;
2270   fileName=MyGlobals::_File_Names[iold];
2271   if (MyGlobals::_Verbose>10)
2272     std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
2273   FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
2274   meshName=MyGlobals::_Mesh_Names[iold];
2275
2276   ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
2277                                                               fileName, meshName, 0, fieldName, DT, IT);
2278
2279   ParaMEDMEM::DataArrayDouble* res=f2->getArray();
2280   //to know names of components
2281   std::vector<std::string> browse=BrowseFieldDouble(f2);
2282   std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
2283   if (MyGlobals::_Verbose>10)
2284     std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
2285   MyGlobals::_General_Informations.push_back(localFieldInformation);
2286   res->incrRef();
2287   f2->decrRef();
2288   _map_dataarray_double[descriptionIold]=res;
2289   return res;
2290 }
2291
2292 void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions()
2293 //to have unique valid fields names/pointers/descriptions for partitionning
2294 //filter _field_descriptions to be in all procs compliant and equal
2295 {
2296   int nbfiles=MyGlobals::_File_Names.size(); //nb domains
2297   if (nbfiles==0)
2298     {
2299       nbfiles=_topology->nbDomain();
2300     }
2301   std::vector<std::string> r2;
2302   //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
2303   for (int i=0; i<(int)_field_descriptions.size(); i++)
2304     {
2305       std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
2306       for (int ii=0; ii<(int)r1.size(); ii++)
2307         r2.push_back(r1[ii]);
2308     }
2309   //here vector(procs*fields) of serialised vector(description) data
2310   _field_descriptions=r2;
2311   int nbfields=_field_descriptions.size(); //on all domains
2312   if ((nbfields%nbfiles)!=0)
2313     {
2314       if (MyGlobals::_Rank==0)
2315         {
2316           std::cerr<< "\nERROR : incoherent number of fields references in all files .med\n" << std::endl
2317               << "fileMedNames :" << std::endl
2318               << ReprVectorOfString(MyGlobals::_File_Names)
2319               << "field_descriptions :" << std::endl
2320               << ReprVectorOfString(MyGlobals::_Field_Descriptions);
2321         }
2322       throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
2323     }
2324   _field_descriptions.resize(nbfields/nbfiles);
2325   for (int i=0; i<(int)_field_descriptions.size(); i++)
2326     {
2327       std::string str=_field_descriptions[i];
2328       str=EraseTagSerialized(str,"idomain=");
2329       str=EraseTagSerialized(str,"fileName=");
2330       _field_descriptions[i]=str;
2331     }
2332 }
2333
2334 //returns true if inodes of a face are in inodes of a cell
2335 bool isFaceOncell(std::vector< int >& inodesFace, std::vector< int >&  inodesCell)
2336 {
2337   int ires=0;
2338   int nbok=inodesFace.size();
2339   for (int i=0; i<nbok; i++)
2340     {
2341       int ii=inodesFace[i];
2342       if (ii<0)
2343         std::cout << "isFaceOncell problem inodeface<0" << std::endl;
2344       for (int j=0; j<(int)inodesCell.size(); j++)
2345         {
2346           if (ii==inodesCell[j])
2347             {
2348               ires=ires+1;
2349               break; //inode of face found
2350             }
2351         }
2352       if (ires<i+1)
2353         break; //inode of face not found do not continue...
2354     }
2355   return (ires==nbok);
2356 }
2357
2358 void MEDPARTITIONER::MeshCollection::filterFaceOnCell()
2359 {
2360   for (int inew=0; inew<_topology->nbDomain(); inew++)
2361     {
2362       if (!isParallelMode() || _domain_selector->isMyDomain(inew))
2363         {
2364           if (MyGlobals::_Verbose>200)
2365             std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
2366           ParaMEDMEM::MEDCouplingUMesh* mcel=_mesh[inew];
2367           ParaMEDMEM::MEDCouplingUMesh* mfac=_face_mesh[inew];
2368
2369           //to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
2370           std::vector<int> nodeIds;
2371           getNodeIds(*mcel, *mfac, nodeIds);
2372           if (nodeIds.size()==0)
2373             continue;  //one empty mesh nothing to do
2374
2375           ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New();
2376           ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::DataArrayInt::New();
2377           mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
2378           int *revC=revNodalCel->getPointer();
2379           int *revIndxC=revNodalIndxCel->getPointer();
2380
2381           std::vector< int > faceOnCell;
2382           std::vector< int > faceNotOnCell;
2383           int nbface=mfac->getNumberOfCells();
2384           for (int iface=0; iface<nbface; iface++)
2385             {
2386               bool ok;
2387               std::vector< int > inodesFace;
2388               mfac->getNodeIdsOfCell(iface, inodesFace);
2389               int nbnodFace=inodesFace.size();
2390               if ( nbnodFace != mfac->getNumberOfNodesInCell( iface ))
2391                 continue; // invalid node ids
2392               //set inodesFace in mcel
2393               int nbok = 0;
2394               for (int i=0; i<nbnodFace; i++)
2395                 nbok += (( inodesFace[i]=nodeIds[inodesFace[i]] ) >= 0 );
2396               if ( nbok != nbnodFace )
2397                 continue;
2398               int inod=inodesFace[0];
2399               if (inod<0)
2400                 {
2401                   std::cout << "filterFaceOnCell problem 1" << std::endl;
2402                   continue;
2403                 }
2404               int nbcell=revIndxC[inod+1]-revIndxC[inod];
2405               for (int j=0; j<nbcell; j++) //look for each cell with inod
2406                 {
2407                   int icel=revC[revIndxC[inod]+j];
2408                   std::vector< int > inodesCell;
2409                   mcel->getNodeIdsOfCell(icel, inodesCell);
2410                   ok=isFaceOncell(inodesFace, inodesCell);
2411                   if (ok) break;
2412                 }
2413               if (ok)
2414                 {
2415                   faceOnCell.push_back(iface);
2416                 }
2417               else
2418                 {
2419                   faceNotOnCell.push_back(iface);
2420                   if (MyGlobals::_Is0verbose>300)
2421                     std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
2422                 }
2423             }
2424
2425           revNodalCel->decrRef();
2426           revNodalIndxCel->decrRef();
2427
2428           // std::string keyy;
2429           // keyy=Cle1ToStr("filterFaceOnCell",inew);
2430           // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
2431           // keyy=Cle1ToStr("filterNotFaceOnCell",inew);
2432           // _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
2433
2434           // filter the face mesh
2435           if ( faceOnCell.empty() )
2436             _face_mesh[inew] = CreateEmptyMEDCouplingUMesh();
2437           else
2438             _face_mesh[inew] = (ParaMEDMEM::MEDCouplingUMesh *)
2439               mfac->buildPartOfMySelf( &faceOnCell[0], &faceOnCell[0] + faceOnCell.size(),true);
2440           mfac->decrRef();
2441
2442           // filter the face families
2443           std::string key = Cle1ToStr("faceFamily_toArray",inew);
2444           if ( getMapDataArrayInt().count( key ))
2445             {
2446               ParaMEDMEM::DataArrayInt * &     fam = getMapDataArrayInt()[ key ];
2447               ParaMEDMEM::DataArrayInt * famFilter = ParaMEDMEM::DataArrayInt::New();
2448               famFilter->alloc(faceOnCell.size(),1);
2449               int* pfamFilter = famFilter->getPointer();
2450               int* pfam       = fam->getPointer();
2451               for ( size_t i=0; i<faceOnCell.size(); i++ )
2452                 pfamFilter[i]=pfam[faceOnCell[i]];
2453               fam->decrRef();
2454               fam = famFilter;
2455             }
2456         }
2457     }
2458 }