src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx

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