src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx

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