src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx

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