src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx

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