src/MEDPartitioner/MEDPARTITIONER_MeshCollection.cxx

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