src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx

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