Back to index

salome-med  6.5.0
MEDPARTITIONER_MeshCollection.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDPARTITIONER_MeshCollection.hxx"
00021 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
00022 #include "MEDPARTITIONER_MeshCollectionMedXmlDriver.hxx"
00023 #include "MEDPARTITIONER_MeshCollectionMedAsciiDriver.hxx"
00024 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
00025 #include "MEDPARTITIONER_Topology.hxx"
00026 #include "MEDPARTITIONER_ParallelTopology.hxx"
00027 
00028 #ifdef HAVE_MPI2
00029 #include "MEDPARTITIONER_JointFinder.hxx"
00030 #endif
00031 
00032 #include "MEDPARTITIONER_Graph.hxx"
00033 #include "MEDPARTITIONER_UserGraph.hxx"
00034 #include "MEDPARTITIONER_Utils.hxx" 
00035 
00036 #include "MEDLoader.hxx"
00037 #include "MEDCouplingMemArray.hxx"
00038 #include "MEDCouplingUMesh.hxx"
00039 #include "MEDCouplingNormalizedUnstructuredMesh.hxx"
00040 #include "MEDCouplingFieldDouble.hxx"
00041 #include "PointLocator3DIntersectorP0P0.hxx"
00042 #include "BBTree.txx"
00043 
00044 #ifdef HAVE_MPI2
00045 #include <mpi.h>
00046 #endif
00047 
00048 #if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS)
00049 #include "MEDPARTITIONER_MetisGraph.hxx"
00050 #endif
00051 
00052 #ifdef MED_ENABLE_SCOTCH
00053 #include "MEDPARTITIONER_ScotchGraph.hxx"
00054 #endif
00055 
00056 #include <set>
00057 #include <vector>
00058 #include <string>
00059 #include <limits>
00060 #include <iostream>
00061 #include <fstream>
00062 
00063 MEDPARTITIONER::MeshCollection::MeshCollection()
00064   : _topology(0),
00065     _owns_topology(false),
00066     _driver(0),
00067     _domain_selector( 0 ),
00068     _i_non_empty_mesh(-1),
00069     _driver_type(MEDPARTITIONER::MedXml),
00070     _subdomain_boundary_creates(false),
00071     _family_splitting(false),
00072     _create_empty_groups(false),
00073     _joint_finder(0)
00074 {
00075 }
00076 
00088 MEDPARTITIONER::MeshCollection::MeshCollection(MeshCollection& initialCollection, 
00089                                                Topology* topology, 
00090                                                bool family_splitting, 
00091                                                bool create_empty_groups)
00092   : _topology(topology),
00093     _owns_topology(false),
00094     _driver(0),
00095     _domain_selector( initialCollection._domain_selector ),
00096     _i_non_empty_mesh(-1),
00097     _name(initialCollection._name),
00098     _driver_type(MEDPARTITIONER::MedXml),
00099     _subdomain_boundary_creates(false),
00100     _family_splitting(family_splitting),
00101     _create_empty_groups(create_empty_groups),
00102     _joint_finder(0)
00103 {
00104   std::vector<std::vector<std::vector<int> > > new2oldIds(initialCollection.getTopology()->nbDomain());
00105   castCellMeshes(initialCollection, new2oldIds);
00106 
00107   //defining the name for the collection and the underlying meshes
00108   setName(initialCollection.getName());
00109 
00111   //treating faces
00113 
00114 #ifdef HAVE_MPI2
00115   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
00116     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
00117 #endif
00118   if (MyGlobals::_Is0verbose)
00119     std::cout<<"treating faces"<<std::endl;
00120   NodeMapping nodeMapping;
00121   //nodeMapping contains the mapping between old nodes and new nodes
00122   // (iolddomain,ioldnode)->(inewdomain,inewnode)
00123   createNodeMapping(initialCollection, nodeMapping);
00124   std::vector<std::vector<std::vector<int> > > new2oldFaceIds;
00125   castFaceMeshes(initialCollection, nodeMapping, new2oldFaceIds);
00126 
00128   //treating families
00130 #ifdef HAVE_MPI2
00131   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
00132     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
00133 #endif
00134   if (MyGlobals::_Is0verbose)
00135     {
00136       if (isParallelMode())
00137         std::cout << "ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
00138       else
00139         std::cout << "NOT ParallelMode on " << topology->nbDomain() << " Domains" << std::endl;
00140     }
00141   if (MyGlobals::_Is0verbose>10)
00142     std::cout<<"treating cell and face families"<<std::endl;
00143   
00144   castIntField2(initialCollection.getMesh(),
00145                 this->getMesh(),
00146                 initialCollection.getCellFamilyIds(),
00147                 "cellFamily");
00148   castIntField2(initialCollection.getFaceMesh(), 
00149                 this->getFaceMesh(),
00150                 initialCollection.getFaceFamilyIds(),
00151                 "faceFamily");
00152 
00153   //treating groups
00154 #ifdef HAVE_MPI2
00155   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
00156     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
00157 #endif
00158   if (MyGlobals::_Is0verbose)
00159     std::cout << "treating groups" << std::endl;
00160   _family_info=initialCollection.getFamilyInfo();
00161   _group_info=initialCollection.getGroupInfo();
00162   
00163 #ifdef HAVE_MPI2
00164   if (MyGlobals::_Verbose>0 && MyGlobals::_World_Size>1)
00165     MPI_Barrier(MPI_COMM_WORLD); //synchronize verbose messages
00166 #endif
00167   if (MyGlobals::_Is0verbose)
00168     std::cout << "treating fields" << std::endl;
00169   castAllFields(initialCollection,"cellFieldDouble");
00170   if (_i_non_empty_mesh<0)
00171     {
00172       for (int i=0; i<_mesh.size(); i++)
00173         {
00174           if (_mesh[i])
00175             {
00176               _i_non_empty_mesh=i; //first existing one local
00177               break;
00178             }
00179         }
00180     }
00181 
00182 }
00183 
00190 void MEDPARTITIONER::MeshCollection::castCellMeshes(MeshCollection& initialCollection,
00191                                                     std::vector<std::vector<std::vector<int> > >& new2oldIds)
00192 {
00193   if (MyGlobals::_Verbose>10)
00194     std::cout << "proc " << MyGlobals::_Rank << " : castCellMeshes" << std::endl;
00195   if (_topology==0)
00196     throw INTERP_KERNEL::Exception("Topology has not been defined on call to castCellMeshes");
00197   
00198   int nbNewDomain=_topology->nbDomain();
00199   int nbOldDomain=initialCollection.getTopology()->nbDomain();
00200   
00201   _mesh.resize(nbNewDomain);
00202   int rank=MyGlobals::_Rank;
00203   //splitting the initial domains into smaller bits
00204   std::vector<std::vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
00205   splitMeshes.resize(nbNewDomain);
00206   for (int inew=0; inew<nbNewDomain; inew++)
00207     {
00208       splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
00209     }
00210 
00211   for (int iold=0; iold<nbOldDomain; iold++)
00212     {
00213       if (!isParallelMode() || initialCollection._domain_selector->isMyDomain(iold))
00214         {
00215           int size=(initialCollection._mesh)[iold]->getNumberOfCells();
00216           std::vector<int> globalids(size);
00217           initialCollection.getTopology()->getCellList(iold, &globalids[0]);
00218           std::vector<int> ilocalnew(size); //local
00219           std::vector<int> ipnew(size);     //idomain old
00220           _topology->convertGlobalCellList(&globalids[0],size,&ilocalnew[0],&ipnew[0]);
00221       
00222           new2oldIds[iold].resize(nbNewDomain);
00223           for (int i=0; i<(int)ilocalnew.size(); i++)
00224             {
00225               new2oldIds[iold][ipnew[i]].push_back(i);
00226             }
00227           for (int inew=0; inew<nbNewDomain; inew++)
00228             {
00229               splitMeshes[inew][iold]=(ParaMEDMEM::MEDCouplingUMesh*)
00230                 (initialCollection.getMesh())[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
00231                                                                        &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
00232                                                                        true);
00233               if (MyGlobals::_Verbose>400)
00234                 std::cout<< "proc " << rank << " : a splitMesh iold inew NbCells " << iold << " " << inew << " "
00235                          << splitMeshes[inew][iold]->getNumberOfCells() << std::endl;
00236             }
00237         }
00238     }
00239 #ifdef HAVE_MPI2
00240   if (isParallelMode())
00241     {
00242       //if (MyGlobals::_Verbose>300) std::cout<<"proc "<<rank<<" : castCellMeshes send/receive"<<std::endl;
00243       for (int iold=0; iold<nbOldDomain; iold++)
00244         for(int inew=0; inew<nbNewDomain; inew++)
00245           {
00246             if (initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
00247               continue;
00248 
00249             if(initialCollection._domain_selector->isMyDomain(iold))
00250               _domain_selector->sendMesh(*(splitMeshes[inew][iold]),_domain_selector->getProcessorID(inew));
00251 
00252             if (_domain_selector->isMyDomain(inew))
00253               _domain_selector->recvMesh(splitMeshes[inew][iold],_domain_selector->getProcessorID(iold));
00254 
00255           }
00256     }
00257 #endif
00258 
00259   //fusing the split meshes
00260   if (MyGlobals::_Verbose>200)
00261     std::cout << "proc " << rank << " : castCellMeshes fusing" << std::endl;
00262   for (int inew=0; inew<nbNewDomain ;inew++)
00263     {
00264       std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
00265     
00266       for (int i=0; i<(int)splitMeshes[inew].size(); i++)
00267         if (splitMeshes[inew][i]!=0) 
00268           if (splitMeshes[inew][i]->getNumberOfCells()>0)
00269             meshes.push_back(splitMeshes[inew][i]);
00270 
00271       if (!isParallelMode()||_domain_selector->isMyDomain(inew))
00272         {
00273           if (meshes.size()==0) 
00274             {
00275               _mesh[inew]=CreateEmptyMEDCouplingUMesh();
00276               std::cout << "WARNING : castCellMeshes fusing : no meshes try another number of processors" << std::endl;
00277             }
00278           else
00279             {
00280               _mesh[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(meshes);
00281               bool areNodesMerged;
00282               int nbNodesMerged;
00283               ParaMEDMEM::DataArrayInt* array=_mesh[inew]->mergeNodes(1e-12,areNodesMerged,nbNodesMerged);
00284               array->decrRef(); // array is not used in this case
00285               _mesh[inew]->zipCoords();
00286             }
00287         }
00288       for (int i=0;i<(int)splitMeshes[inew].size();i++)
00289         if (splitMeshes[inew][i]!=0)
00290           splitMeshes[inew][i]->decrRef();
00291     }  
00292   if (MyGlobals::_Verbose>300)
00293     std::cout << "proc " << rank << " : castCellMeshes end fusing" << std::endl;
00294 }
00295 
00300 void MEDPARTITIONER::MeshCollection::createNodeMapping( MeshCollection& initialCollection, NodeMapping& nodeMapping)
00301 {
00302   using std::vector;
00303   using std::make_pair;
00304   //  NodeMapping reverseNodeMapping;
00305   for (int iold=0; iold<initialCollection.getTopology()->nbDomain();iold++)
00306     {
00307 
00308       double* bbox;
00309       BBTree<3>* tree; 
00310       if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
00311         {
00312           //      std::map<pair<double,pair<double, double> >, int > nodeClassifier;
00313           int nvertices=initialCollection.getMesh(iold)->getNumberOfNodes();
00314           bbox=new double[nvertices*6];
00315           ParaMEDMEM::DataArrayDouble* coords = initialCollection.getMesh(iold)->getCoords();
00316           double* coordsPtr=coords->getPointer();
00317 
00318           for (int i=0; i<nvertices*3;i++)
00319             {
00320               bbox[i*2]=coordsPtr[i]-1e-6;
00321               bbox[i*2+1]=coordsPtr[i]+1e-6;
00322             }
00323           tree=new BBTree<3>(bbox,0,0,nvertices,1e-9);
00324         }
00325 
00326       for (int inew=0; inew<_topology->nbDomain(); inew++)
00327         {
00328 #ifdef HAVE_MPI2
00329           //sending meshes for parallel computation
00330           if (isParallelMode() && _domain_selector->isMyDomain(inew) && !_domain_selector->isMyDomain(iold))  
00331             _domain_selector->sendMesh(*(getMesh(inew)), _domain_selector->getProcessorID(iold));
00332           else if (isParallelMode() && !_domain_selector->isMyDomain(inew)&& _domain_selector->isMyDomain(iold))
00333             {
00334               ParaMEDMEM::MEDCouplingUMesh* mesh;
00335               _domain_selector->recvMesh(mesh, _domain_selector->getProcessorID(inew));
00336               ParaMEDMEM::DataArrayDouble* coords = mesh->getCoords();
00337               for (int inode=0; inode<mesh->getNumberOfNodes();inode++)
00338                 {
00339                   double* coordsPtr=coords->getPointer()+inode*3;
00340                   vector<int> elems;
00341                   tree->getElementsAroundPoint(coordsPtr,elems);
00342                   if (elems.size()==0) continue;
00343                   nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
00344                 }
00345               mesh->decrRef();
00346             }
00347           else if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
00348 #else
00349           if (!isParallelMode() || (_domain_selector->isMyDomain(inew) && _domain_selector->isMyDomain(iold)))
00350 #endif
00351             {
00352               ParaMEDMEM::DataArrayDouble* coords = getMesh(inew)->getCoords();
00353               for (int inode=0; inode<_mesh[inew]->getNumberOfNodes();inode++)
00354                 {
00355                   double* coordsPtr=coords->getPointer()+inode*3;
00356                   vector<int> elems;
00357                   tree->getElementsAroundPoint(coordsPtr,elems);
00358                   if (elems.size()==0) continue;
00359                   nodeMapping.insert(make_pair(make_pair(iold,elems[0]),make_pair(inew,inode)));
00360                 }
00361             }
00362         }
00363       if (!isParallelMode() || (_domain_selector->isMyDomain(iold)))
00364         {
00365           delete tree;
00366           delete[] bbox;
00367         }
00368     } 
00369 
00370 }
00371 
00372 void getNodeIds(ParaMEDMEM::MEDCouplingUMesh& meshOne, ParaMEDMEM::MEDCouplingUMesh& meshTwo, std::vector<int>& nodeIds)
00373 {
00374   using std::vector;
00375   if (!&meshOne || !&meshTwo) return;  //empty or not existing
00376   double* bbox;
00377   BBTree<3>* tree;
00378   int nv1=meshOne.getNumberOfNodes();
00379   bbox=new double[nv1*6];
00380   ParaMEDMEM::DataArrayDouble* coords=meshOne.getCoords();
00381   double* coordsPtr=coords->getPointer();
00382   for (int i=0; i<nv1*3; i++)
00383     {
00384       bbox[i*2]=coordsPtr[i]-1e-6;
00385       bbox[i*2+1]=coordsPtr[i]+1e-6;
00386     }
00387   tree=new BBTree<3>(bbox,0,0,nv1,1e-9);
00388   
00389   int nv2=meshTwo.getNumberOfNodes();
00390   nodeIds.resize(nv2,-1);
00391   coords=meshTwo.getCoords();
00392   for (int inode=0; inode<nv2; inode++)
00393     {
00394       double* coordsPtr2=coords->getPointer()+inode*3;
00395       vector<int> elems;
00396       tree->getElementsAroundPoint(coordsPtr2,elems);
00397       if (elems.size()==0) continue;
00398       nodeIds[inode]=elems[0];
00399     }
00400   delete tree;
00401   delete[] bbox;
00402 }
00403 
00408 void MEDPARTITIONER::MeshCollection::castFaceMeshes(MeshCollection& initialCollection,
00409                                                     const std::multimap<std::pair<int,int>, std::pair<int,int> >& nodeMapping,
00410                                                     std::vector<std::vector<std::vector<int> > >& new2oldIds)
00411 {
00412 
00413   //splitMeshes structure will contain the partition of 
00414   //the old faces on the new ones
00415   //splitMeshes[4][2] contains the faces from old domain 2
00416   //that have to be added to domain 4
00417   
00418   using std::vector;
00419   using std::map;
00420   using std::multimap;
00421   using std::pair;
00422   using std::make_pair;
00423   
00424   if (MyGlobals::_Verbose>10)
00425     std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes" << std::endl;
00426   if (_topology==0)
00427     throw INTERP_KERNEL::Exception("Topology has not been defined on call to castFaceMeshes");
00428   
00429   int nbNewDomain=_topology->nbDomain();
00430   int nbOldDomain=initialCollection.getTopology()->nbDomain();
00431   
00432   vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom=initialCollection.getFaceMesh();
00433   vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo=this->getFaceMesh();
00434   
00435   vector< vector<ParaMEDMEM::MEDCouplingUMesh*> > splitMeshes;
00436   
00437   splitMeshes.resize(nbNewDomain);
00438   for (int inew=0; inew<nbNewDomain; inew++)
00439     {
00440       splitMeshes[inew].resize(nbOldDomain, (ParaMEDMEM::MEDCouplingUMesh*)0);
00441     }
00442   new2oldIds.resize(nbOldDomain);
00443   for (int iold=0; iold<nbOldDomain; iold++) new2oldIds[iold].resize(nbNewDomain);
00444   
00445   //init null pointer for empty meshes
00446   for (int inew=0; inew<nbNewDomain; inew++)
00447     {
00448       for (int iold=0; iold<nbOldDomain; iold++)
00449         {
00450           splitMeshes[inew][iold]=0; //null for empty meshes
00451           new2oldIds[iold][inew].clear();
00452         }
00453     }
00454 
00455   //loop over the old domains to analyse the faces and decide 
00456   //on which new domain they belong
00457   for (int iold=0; iold<nbOldDomain; iold++)
00458     {
00459       if (isParallelMode() && !initialCollection._domain_selector->isMyDomain(iold)) continue;
00460       if (MyGlobals::_Verbose>400)
00461         std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes iodDomain "<<iold<<std::endl;
00462       //initial face mesh known : in my domain
00463       if (meshesCastFrom[iold] != 0)
00464         {
00465           for (int ielem=0; ielem<meshesCastFrom[iold]->getNumberOfCells(); ielem++)
00466             {
00467               vector<int> nodes;
00468               meshesCastFrom[iold]->getNodeIdsOfCell(ielem,nodes);
00469           
00470               map <int,int> faces;
00471 
00472               //analysis of element ielem
00473               //counters are set for the element
00474               //for each source node, the mapping is interrogated and the domain counters 
00475               //are incremented for each target node
00476               //the face is considered as going to target domains if the counter of the domain 
00477               //is equal to the number of nodes
00478               for (int inode=0; inode<(int)nodes.size(); inode++)
00479                 {
00480                   typedef multimap<pair<int,int>,pair<int,int> >::const_iterator MI;
00481                   int mynode=nodes[inode];
00482 
00483                   pair <MI,MI> myRange = nodeMapping.equal_range(make_pair(iold,mynode));
00484                   for (MI iter=myRange.first; iter!=myRange.second; iter++)
00485                     {
00486                       int inew=iter->second.first;
00487                       if (faces.find(inew)==faces.end())
00488                         faces[inew]=1;
00489                       else
00490                         faces[inew]++;
00491                     }
00492                 }
00493           
00494               for (map<int,int>::iterator iter=faces.begin(); iter!=faces.end(); iter++)
00495                 {
00496                   if (iter->second==(int)nodes.size())
00497                     //cvw eligible but may be have to be face of a cell of this->getMesh()[inew]?
00498                     //it is not sure here...
00499                     //done before writeMedfile on option?... see filterFaceOnCell()
00500                     new2oldIds[iold][iter->first].push_back(ielem);
00501                 }
00502             }
00503       
00504           //creating the splitMeshes from the face ids
00505           for (int inew=0; inew<nbNewDomain; inew++)
00506             {
00507               if (meshesCastFrom[iold]->getNumberOfCells() > 0)
00508                 {
00509                   splitMeshes[inew][iold]=
00510                     (ParaMEDMEM::MEDCouplingUMesh*) 
00511                     ( meshesCastFrom[iold]->buildPartOfMySelf(&new2oldIds[iold][inew][0],
00512                                                               &new2oldIds[iold][inew][0]+new2oldIds[iold][inew].size(),
00513                                                               true) 
00514                       );
00515                   splitMeshes[inew][iold]->zipCoords();
00516                 }
00517               else
00518                 {
00519                   splitMeshes[inew][iold]=CreateEmptyMEDCouplingUMesh();
00520                 }
00521             }
00522         }
00523       else
00524         {
00525           std::cout<<"proc "<<MyGlobals::_Rank<<" : castFaceMeshes empty mesh from iodDomain "<<iold<<std::endl;
00526         }
00527     }
00528   
00529 #ifdef HAVE_MPI2
00530   //send/receive stuff
00531   if (isParallelMode())
00532     {
00533       ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
00534       for (int iold=0; iold<nbOldDomain; iold++)
00535         for (int inew=0; inew<nbNewDomain; inew++)
00536           {
00537             if (initialCollection._domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
00538               {
00539                 if (splitMeshes[inew][iold] != 0)
00540                   {
00541                     _domain_selector->sendMesh(*(splitMeshes[inew][iold]), _domain_selector->getProcessorID(inew));
00542                     //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes send " <<inew<<" "<<iold<<" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
00543                   }
00544                 else
00545                   {
00546                     _domain_selector->sendMesh(*(empty), _domain_selector->getProcessorID(inew));
00547                     //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes sen0 " <<inew<<" "<<iold<<std::endl;
00548                   }
00549               }
00550             if (!initialCollection._domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
00551               _domain_selector->recvMesh(splitMeshes[inew][iold], _domain_selector->getProcessorID(iold));
00552               int nb=0;
00553               if (splitMeshes[inew][iold])
00554                 nb=splitMeshes[inew][iold]->getNumberOfCells();
00555               //std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes recv "<<inew<<" "<<iold<<" "<<nb<<std::endl;//" "<<splitMeshes[inew][iold]->getNumberOfCells()<<std::endl;
00556           }
00557       empty->decrRef();
00558     }
00559 #endif
00560 
00561   //fusing the split meshes
00562   if (MyGlobals::_Verbose>200)
00563     std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes fusing" << std::endl;
00564   meshesCastTo.resize(nbNewDomain);
00565   for (int inew=0; inew<nbNewDomain; inew++)
00566     {
00567       vector<const ParaMEDMEM::MEDCouplingUMesh*> myMeshes;
00568       for (int iold=0; iold<nbOldDomain; iold++)
00569         {
00570           ParaMEDMEM::MEDCouplingUMesh *umesh=splitMeshes[inew][iold];
00571           if (umesh!=0)
00572             if (umesh->getNumberOfCells()>0)
00573                 myMeshes.push_back(umesh);
00574         }
00575       
00576       if (myMeshes.size()>0)
00577         {
00578           meshesCastTo[inew]=ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes(myMeshes);
00579         }
00580       else
00581         {
00582           ParaMEDMEM::MEDCouplingUMesh *empty=CreateEmptyMEDCouplingUMesh();
00583           meshesCastTo[inew]=empty;
00584         }
00585       for (int iold=0; iold<nbOldDomain; iold++)
00586         if (splitMeshes[inew][iold]!=0)
00587           splitMeshes[inew][iold]->decrRef();
00588     }
00589   if (MyGlobals::_Verbose>300)
00590     std::cout << "proc " << MyGlobals::_Rank << " : castFaceMeshes end fusing" << std::endl;
00591 }
00592 
00593 void MEDPARTITIONER::MeshCollection::remapIntField(const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
00594                                    const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
00595                                    const int* fromArray,
00596                                    int* toArray)
00597 {
00598   using std::vector;
00599   if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
00600   //cvw std::cout<<"remapIntField "<<sourceMesh.getNumberOfCells()<<" "<<targetMesh.getNumberOfCells()<<std::endl;
00601   ParaMEDMEM::DataArrayDouble* sourceCoords=sourceMesh.getBarycenterAndOwner();
00602   ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
00603    
00604   ParaMEDMEM::MEDCouplingUMesh* tmpMesh=ParaMEDMEM::MEDCouplingUMesh::New();
00605   tmpMesh->setCoords(sourceCoords);
00606   vector<int> c;
00607   vector<int> cI;
00608   tmpMesh->getNodeIdsNearPoints(targetCoords->getConstPointer(),targetMesh.getNumberOfCells(),1e-10,c,cI);
00609   if ((int)cI.size()!= targetMesh.getNumberOfCells()+1) 
00610     throw INTERP_KERNEL::Exception("Error in source/target projection");
00611   for (int itargetnode=0; itargetnode<targetMesh.getNumberOfCells();itargetnode++)    
00612     {
00613       if (cI[itargetnode]==cI[itargetnode+1])
00614         continue;
00615       int isourcenode=c[cI[itargetnode]];
00616       toArray[itargetnode]=fromArray[isourcenode];
00617     } 
00618   sourceCoords->decrRef();
00619   targetCoords->decrRef();
00620   tmpMesh->decrRef();
00621 }
00622 
00623 void MEDPARTITIONER::MeshCollection::castIntField2(std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastFrom,
00624                                    std::vector<ParaMEDMEM::MEDCouplingUMesh*>& meshesCastTo,
00625                                    std::vector<ParaMEDMEM::DataArrayInt*>& arrayFrom,
00626                                    std::string nameArrayTo)
00627 {
00628   using std::vector;
00629   
00630   int ioldMax=meshesCastFrom.size();
00631   int inewMax=meshesCastTo.size();
00632   // send-recv operations
00633 #ifdef HAVE_MPI2
00634   for (int inew=0; inew<inewMax; inew++)
00635     {
00636       for (int iold=0; iold<ioldMax; iold++)
00637         {
00638           //sending arrays for distant domains
00639           if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
00640             {
00641               //send mesh
00642               _domain_selector->sendMesh(*meshesCastFrom[iold],_domain_selector->getProcessorID(inew));
00643               //send vector
00644               int size=arrayFrom[iold]->getNumberOfTuples(); //cvw may be -1!
00645               vector<int>sendIds;
00646               if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField SendIntVec size "<<size<<std::endl;
00647               if (size>0) //no empty
00648                 {
00649                   sendIds.resize(size);
00650                   std::copy(arrayFrom[iold]->getPointer(),arrayFrom[iold]->getPointer()+size,&sendIds[0]);
00651                 }
00652               else //empty
00653                 {
00654                   size=0;
00655                   sendIds.resize(size);
00656                 }
00657               SendIntVec(sendIds, _domain_selector->getProcessorID(inew));
00658             }
00659           //receiving arrays from distant domains
00660           if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
00661             {
00662               //receive mesh
00663               vector<int> recvIds;
00664               ParaMEDMEM::MEDCouplingUMesh* recvMesh;
00665               _domain_selector->recvMesh(recvMesh,_domain_selector->getProcessorID(iold));
00666               //receive vector
00667               if (MyGlobals::_Verbose>400) std::cout<<"proc "<<_domain_selector->rank()<<" : castIntField recIntVec "<<std::endl;
00668               RecvIntVec(recvIds, _domain_selector->getProcessorID(iold));
00669               remapIntField2(inew,iold,*recvMesh,*meshesCastTo[inew],&recvIds[0],nameArrayTo);
00670               recvMesh->decrRef(); //cww is it correct?
00671             }
00672         }
00673     }
00674 #endif
00675   
00676   //local contributions and aggregation
00677   for (int inew=0; inew<inewMax; inew++)    
00678     {
00679       for (int iold=0; iold<ioldMax; iold++)
00680         if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
00681           {
00682             remapIntField2(inew,iold,*meshesCastFrom[iold],*meshesCastTo[inew],arrayFrom[iold]->getConstPointer(),nameArrayTo);
00683           }
00684     }
00685 }
00686 
00687 void MEDPARTITIONER::MeshCollection::remapIntField2(int inew, int iold,
00688                                                     const ParaMEDMEM::MEDCouplingUMesh& sourceMesh,
00689                                                     const ParaMEDMEM::MEDCouplingUMesh& targetMesh,
00690                                                     const int* fromArray,
00691                                                     std::string nameArrayTo)
00692 //here we store ccI for next use in future call of castAllFields and remapDoubleField2
00693 {
00694   if (sourceMesh.getNumberOfCells()<=0) return; //empty mesh could exist
00695   ParaMEDMEM::DataArrayDouble* sourceCoords=sourceMesh.getBarycenterAndOwner();
00696   ParaMEDMEM::DataArrayDouble* targetCoords=targetMesh.getBarycenterAndOwner();
00697    
00698   ParaMEDMEM::MEDCouplingUMesh* tmpMesh=ParaMEDMEM::MEDCouplingUMesh::New();
00699   tmpMesh->setCoords(sourceCoords);
00700   std::vector<int> c;
00701   std::vector<int> cI;
00702   std::vector<int> ccI; //memorize intersection target<-source(inew,iold)
00703   std::string str,cle;
00704   str=nameArrayTo+"_toArray";
00705   cle=Cle1ToStr(str,inew);
00706   int* toArray;
00707   int targetSize=targetMesh.getNumberOfCells();
00708   //first time iold : create and initiate 
00709   if (_map_dataarray_int.find(cle)==_map_dataarray_int.end())
00710     {
00711       if (MyGlobals::_Is0verbose>100)
00712         std::cout << "create " << cle << " size " << targetSize << std::endl;
00713       ParaMEDMEM::DataArrayInt* p=ParaMEDMEM::DataArrayInt::New();
00714       p->alloc(targetSize,1);
00715       p->fillWithZero();
00716       toArray=p->getPointer();
00717       _map_dataarray_int[cle]=p;
00718     }
00719   else //other times iold: refind and complete
00720     {
00721       toArray=_map_dataarray_int.find(cle)->second->getPointer();
00722     }
00723   tmpMesh->getNodeIdsNearPoints(targetCoords->getConstPointer(),targetSize,1e-10,c,cI);
00724   if ((int)cI.size()!=targetSize+1) 
00725     throw INTERP_KERNEL::Exception("Error in source/target projection");
00726   for (int itargetnode=0; itargetnode<targetSize; itargetnode++)    
00727     {
00728       if (cI[itargetnode]==cI[itargetnode+1]) continue;
00729       int isourcenode=c[cI[itargetnode]];
00730       toArray[itargetnode]=fromArray[isourcenode];
00731       //memorize intersection 
00732       ccI.push_back(itargetnode); //next we'll do toArray[ccI[i]]=fromArray[ccI[i+1]]
00733       ccI.push_back(isourcenode);
00734     }
00735   //ccI.push_back(sourceMesh.getNumberOfCells()); //additionnal information at end??
00736   
00737   //memories intersection for future same job on fields (if no existing cle=no intersection)
00738   str=Cle2ToStr(nameArrayTo+"_ccI",inew,iold);
00739   if (MyGlobals::_Verbose>700)
00740     std::cout << "proc " << MyGlobals::_Rank << " : map memorize '" << str << "'\n";
00741   _map_dataarray_int[str]=CreateDataArrayIntFromVector(ccI, 2);
00742   sourceCoords->decrRef();
00743   targetCoords->decrRef();
00744   tmpMesh->decrRef();
00745 }
00746 
00747 void MEDPARTITIONER::MeshCollection::castAllFields(MeshCollection& initialCollection, std::string nameArrayTo)
00748 {
00749   if (nameArrayTo!="cellFieldDouble") 
00750     throw INTERP_KERNEL::Exception("Error castAllField only on cellFieldDouble");
00751 
00752   std::string nameTo="typeData=6"; //resume the type of field casted 
00753   // send-recv operations
00754   int ioldMax=initialCollection.getMesh().size();
00755   int inewMax=this->getMesh().size();
00756   int iFieldMax=initialCollection.getFieldDescriptions().size();
00757   if (MyGlobals::_Verbose>10)
00758     std::cout << "castAllFields with:\n" << ReprVectorOfString(initialCollection.getFieldDescriptions()) << std::endl;
00759   //see collection.prepareFieldDescriptions()
00760   for (int ifield=0; ifield<iFieldMax; ifield++)
00761     {
00762       std::string descriptionField=initialCollection.getFieldDescriptions()[ifield];
00763       if (descriptionField.find(nameTo)==std::string::npos)
00764         continue; //only nameTo accepted in Fields name description
00765 #ifdef HAVE_MPI2
00766       for (int inew=0; inew<inewMax; inew++)
00767         {
00768           for (int iold=0; iold<ioldMax; iold++)
00769             {
00770               //sending arrays for distant domains
00771               if (isParallelMode() && _domain_selector->isMyDomain(iold) && !_domain_selector->isMyDomain(inew))
00772                 {
00773                   int target=_domain_selector->getProcessorID(inew);
00774                   ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
00775                   if (MyGlobals::_Verbose>10) 
00776                     std::cout << "proc " << _domain_selector->rank() << " : castAllFields sendDouble" << std::endl;
00777                   SendDataArrayDouble(field, target);
00778                 }
00779               //receiving arrays from distant domains
00780               if (isParallelMode() && !_domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew))
00781                 {
00782                   int source=_domain_selector->getProcessorID(iold);
00783                   //receive vector
00784                   if (MyGlobals::_Verbose>10) 
00785                     std::cout << "proc " << _domain_selector->rank() << " : castAllFields recvDouble" << std::endl;
00786                   ParaMEDMEM::DataArrayDouble* field=RecvDataArrayDouble(source);
00787                   remapDoubleField3(inew,iold,field,nameArrayTo,descriptionField);
00788                 }
00789             }
00790         }
00791 #endif
00792       //local contributions and aggregation
00793       for (int inew=0; inew<inewMax; inew++)
00794         {
00795           for (int iold=0; iold<ioldMax; iold++)
00796             if (!isParallelMode() || ( _domain_selector->isMyDomain(iold) && _domain_selector->isMyDomain(inew)))
00797               {
00798                 ParaMEDMEM::DataArrayDouble* field=initialCollection.getField(descriptionField,iold);
00799                 remapDoubleField3(inew,iold,field,nameArrayTo,descriptionField);
00800               }
00801         }
00802     }
00803 }
00804 
00805 void MEDPARTITIONER::MeshCollection::remapDoubleField3(int inew, int iold, 
00806                                                        ParaMEDMEM::DataArrayDouble* fromArray,
00807                                                        std::string nameArrayTo,
00808                                                        std::string descriptionField)
00809 //here we use 'cellFamily_ccI inew iold' set in remapIntField2
00810 {
00811   if (nameArrayTo!="cellFieldDouble") 
00812     throw INTERP_KERNEL::Exception("Error remapDoubleField3 only on cellFieldDouble");
00813   std::string key=Cle2ToStr("cellFamily_ccI",inew,iold);
00814   
00815   std::map<std::string,ParaMEDMEM::DataArrayInt*>::iterator it1;
00816   it1=_map_dataarray_int.find(key);
00817   if (it1==_map_dataarray_int.end())
00818     {
00819       std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 key '" << key << "' not found" << std::endl;
00820       std::cerr << " trying remap of field double on cells : " << descriptionField << std::endl;
00821       return;
00822     }
00823   //create ccI in remapIntField2
00824   ParaMEDMEM::DataArrayInt *ccI=it1->second;
00825   if (MyGlobals::_Verbose>300)
00826     std::cout << "proc " << MyGlobals::_Rank << " : remapDoubleField3 " << key << " size " << ccI->getNbOfElems() << std::endl;
00827   
00828   int nbcell=this->getMesh()[inew]->getNumberOfCells();
00829   int nbcomp=fromArray->getNumberOfComponents();
00830   int nbPtGauss=StrToInt(ExtractFromDescription(descriptionField, "nbPtGauss="));
00831   
00832   std::string tag="inewFieldDouble="+IntToStr(inew);
00833   key=descriptionField+SerializeFromString(tag);
00834   int fromArrayNbOfElem=fromArray->getNbOfElems();
00835   int fromArrayNbOfComp=fromArray->getNumberOfComponents();
00836   int fromArrayNbOfCell=fromArrayNbOfElem/fromArrayNbOfComp/nbPtGauss;
00837   
00838   if (MyGlobals::_Verbose>1000)
00839     {
00840       std::cout<<"proc " << MyGlobals::_Rank << " nbcell " << nbcell << " nbcomp " << nbcomp << " nbPtGauss " << nbPtGauss <<
00841         " fromArray nbOfElems " << fromArrayNbOfElem <<
00842         " nbTuples " << fromArray->getNumberOfTuples() <<
00843         " nbcells " << fromArrayNbOfCell <<
00844         " nbComponents " << fromArray->getNumberOfComponents() << std::endl;
00845     }
00846 
00847   ParaMEDMEM::DataArrayDouble* field=0;
00848   std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it2;
00849   it2=_map_dataarray_double.find(key);
00850   if (it2==_map_dataarray_double.end())
00851     {
00852       if (MyGlobals::_Verbose>300)
00853         std::cout << "proc "<< MyGlobals::_Rank << " : remapDoubleField3 key '" << key << "' not found and create it" << std::endl;
00854       field=ParaMEDMEM::DataArrayDouble::New();
00855       _map_dataarray_double[key]=field;
00856       field->alloc(nbcell*nbPtGauss,nbcomp);
00857       field->fillWithZero();
00858     }
00859   else
00860     {
00861       field=it2->second;
00862       if (field->getNumberOfTuples()!=nbcell*nbPtGauss || field->getNumberOfComponents()!=nbcomp)
00863         {
00864           std::cerr << "proc " << MyGlobals::_Rank << " : remapDoubleField3 pb of size " <<
00865             " trying remap of field double on cells : \n" << descriptionField << std::endl;
00866           return;
00867         }
00868     }
00869   
00870   if (nbPtGauss==1)
00871     {
00872       field->setPartOfValuesAdv(fromArray,ccI);
00873     }
00874   else
00875     {
00876       //replaced by setPartOfValuesAdv if nbPtGauss==1
00877       int iMax=ccI->getNbOfElems();
00878       int* pccI=ccI->getPointer();
00879       double* pField=field->getPointer();
00880       double* pFrom=fromArray->getPointer();
00881       int itarget, isource, delta=nbPtGauss*nbcomp;
00882       for (int i=0; i<iMax; i=i+2)  //cell
00883         {
00884           itarget=pccI[i];
00885           isource=pccI[i+1];
00886           if ((itarget<0) || (itarget>=nbcell) || (isource<0) || (isource>=fromArrayNbOfCell))
00887             throw INTERP_KERNEL::Exception("Error field override");
00888           int ita=itarget*delta;
00889           int iso=isource*delta;
00890           for (int k=0; k<delta; k++) pField[ita+k]=pFrom[iso+k]; //components and gausspoints
00891         }
00892     }
00893 }
00894 
00899 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename)
00900   : _topology(0),
00901     _owns_topology(true),
00902     _driver(0),
00903     _domain_selector( 0 ),
00904     _i_non_empty_mesh(-1),
00905     _driver_type(MEDPARTITIONER::Undefined),
00906     _subdomain_boundary_creates(false),
00907     _family_splitting(false),
00908     _create_empty_groups(false),
00909     _joint_finder(0)
00910 {
00911   try
00912     {
00913       _driver=new MeshCollectionMedXmlDriver(this);
00914       _driver->read (filename.c_str());
00915       _driver_type = MedXml;
00916     }
00917   catch(...) 
00918     {  // Handle all exceptions
00919       if ( _driver ) delete _driver; _driver=0;
00920       try
00921         {
00922           _driver=new MeshCollectionMedAsciiDriver(this);
00923           _driver->read (filename.c_str());
00924           _driver_type=MedAscii;
00925         }
00926       catch(...)
00927         {
00928           delete _driver;
00929           _driver=0;
00930           throw INTERP_KERNEL::Exception("file does not comply with any recognized format");
00931         }
00932     }
00933   for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
00934     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
00935       _i_non_empty_mesh = idomain;
00936 }
00937 
00943 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, ParaDomainSelector& domainSelector)
00944   : _topology(0),
00945     _owns_topology(true),
00946     _driver(0),
00947     _domain_selector( &domainSelector ),
00948     _i_non_empty_mesh(-1),
00949     _driver_type(MEDPARTITIONER::Undefined),
00950     _subdomain_boundary_creates(false),
00951     _family_splitting(false),
00952     _create_empty_groups(false),
00953     _joint_finder(0)
00954 {
00955   std::string myfile=filename;
00956   std::size_t found=myfile.find(".xml");
00957   if (found!=std::string::npos) //file .xml
00958     {
00959       try
00960         {
00961           _driver=new MeshCollectionMedXmlDriver(this);
00962           _driver->read ( filename.c_str(), _domain_selector );
00963           _driver_type = MedXml;
00964         }
00965       catch(...)
00966         {  // Handle all exceptions
00967           delete _driver;
00968           throw INTERP_KERNEL::Exception("file .xml does not comply with any recognized format");
00969         }
00970     }
00971   else 
00972     {
00973       found=myfile.find(".med");
00974       if (found!=std::string::npos) //file .med single
00975         {
00976           //make a temporary file .xml and retry MedXmlDriver
00977           std::string xml="\
00978 <?xml version=\"1.0\" encoding=\"UTF-8\"?>\n \
00979 <root>\n \
00980   <version maj=\"2\" min=\"3\" ver=\"1\"/>\n \
00981   <description what=\"\" when=\"\"/>\n \
00982   <content>\n \
00983     <mesh name=\"$meshName\"/>\n \
00984   </content>\n \
00985   <splitting>\n \
00986     <subdomain number=\"1\"/>\n \
00987     <global_numbering present=\"no\"/>\n \
00988   </splitting>\n \
00989   <files>\n \
00990     <subfile id=\"1\">\n \
00991       <name>$fileName</name>\n \
00992       <machine>localhost</machine>\n \
00993     </subfile>\n \
00994   </files>\n \
00995   <mapping>\n \
00996     <mesh name=\"$meshName\">\n \
00997       <chunk subdomain=\"1\">\n \
00998         <name>$meshName</name>\n \
00999       </chunk>\n \
01000     </mesh>\n \
01001   </mapping>\n \
01002 </root>\n";
01003           std::vector<std::string> meshNames=MEDLoader::GetMeshNames(myfile.c_str());
01004           xml.replace(xml.find("$fileName"),9,myfile);
01005           xml.replace(xml.find("$meshName"),9,meshNames[0]);
01006           xml.replace(xml.find("$meshName"),9,meshNames[0]);
01007           xml.replace(xml.find("$meshName"),9,meshNames[0]);
01008           std::string nameFileXml=myfile;
01009           nameFileXml.replace(nameFileXml.find(".med"),4,".xml");
01010           nameFileXml="medpartitioner_"+nameFileXml;
01011           if (_domain_selector->rank()==0) //only on to write it
01012             {
01013               std::ofstream f(nameFileXml.c_str());
01014               f<<xml;
01015               f.close();
01016             }
01017 #ifdef HAVE_MPI2
01018            if (MyGlobals::_World_Size>1)
01019              MPI_Barrier(MPI_COMM_WORLD); //wait for creation of nameFileXml
01020 #endif
01021           try
01022             {
01023               _driver=new MeshCollectionMedXmlDriver(this);
01024               _driver->read ( nameFileXml.c_str(), _domain_selector );
01025               _driver_type = MedXml;
01026             }
01027           catch(...)
01028             {  // Handle all exceptions
01029               if ( _driver ) delete _driver; _driver=0;
01030               throw INTERP_KERNEL::Exception("file medpartitioner_xxx.xml does not comply with any recognized format");
01031             }
01032         }
01033       else //no extension
01034         {
01035           try
01036             {
01037               _driver=new MeshCollectionMedAsciiDriver(this);
01038               _driver->read ( filename.c_str(), _domain_selector );
01039               _driver_type=MedAscii;
01040             }
01041           catch(...)
01042             {
01043               delete _driver;
01044               _driver=0;
01045               throw INTERP_KERNEL::Exception("file name with no extension does not comply with any recognized format");
01046             }
01047         }
01048     }
01049   // find non-empty domain mesh
01050   for ( int idomain = 0; idomain < (int)_mesh.size(); ++idomain )
01051     if ( _mesh[idomain] && _mesh[idomain]->getNumberOfNodes() > 0 )
01052       _i_non_empty_mesh = idomain;
01053   
01054   try
01055     {
01056       //check for all proc/file compatibility of _field_descriptions
01057 #ifdef HAVE_MPI2
01058       _field_descriptions=AllgathervVectorOfString(MyGlobals::_Field_Descriptions);
01059 #else
01060       _field_descriptions=MyGlobals::_Field_Descriptions;
01061 #endif
01062     }
01063   catch(INTERP_KERNEL::Exception& e)
01064     {
01065       std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
01066       throw INTERP_KERNEL::Exception("Something wrong verifying coherency of files med ands fields");
01067     }
01068 #ifdef HAVE_MPI2
01069   try
01070     {
01071       //check for all proc/file compatibility of _family_info
01072       std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringInt(_family_info));
01073       _family_info=DevectorizeToMapOfStringInt(v2);
01074     }
01075   catch(INTERP_KERNEL::Exception& e)
01076     {
01077       std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
01078       throw INTERP_KERNEL::Exception("Something wrong merging all familyInfo");
01079     }
01080 
01081   try
01082     {
01083       //check for all proc/file compatibility of _group_info
01084       std::vector<std::string> v2=AllgathervVectorOfString(VectorizeFromMapOfStringVectorOfString(_group_info));
01085       _group_info=DeleteDuplicatesInMapOfStringVectorOfString(DevectorizeToMapOfStringVectorOfString(v2));
01086     }
01087   catch(INTERP_KERNEL::Exception& e)
01088     {
01089       std::cerr << "proc " << MyGlobals::_Rank << " : INTERP_KERNEL_Exception : " << e.what() << std::endl;
01090       throw INTERP_KERNEL::Exception("Something wrong merging all groupInfo");
01091     }
01092 #endif
01093 }
01094 
01100 MEDPARTITIONER::MeshCollection::MeshCollection(const std::string& filename, const std::string& meshname)
01101   : _topology(0),
01102     _owns_topology(true),
01103     _driver(0),
01104     _domain_selector( 0 ),
01105     _i_non_empty_mesh(-1),
01106     _name(meshname),
01107     _driver_type(MEDPARTITIONER::MedXml),
01108     _subdomain_boundary_creates(false),
01109     _family_splitting(false),
01110     _create_empty_groups(false),
01111     _joint_finder(0)
01112 {
01113   try // avoid memory leak in case of inexistent filename
01114     {
01115       retrieveDriver()->readSeq (filename.c_str(),meshname.c_str());
01116     }
01117   catch (...)
01118     {
01119       delete _driver;
01120       _driver=0;
01121       throw INTERP_KERNEL::Exception("problem reading .med files");
01122     }
01123   if ( _mesh[0] && _mesh[0]->getNumberOfNodes() > 0 )
01124     _i_non_empty_mesh = 0;
01125 }
01126 
01127 MEDPARTITIONER::MeshCollection::~MeshCollection()
01128 {
01129   for (int i=0; i<(int)_mesh.size();i++)
01130     if (_mesh[i]!=0) _mesh[i]->decrRef(); 
01131   
01132   for (int i=0; i<(int)_cell_family_ids.size();i++)
01133     if (_cell_family_ids[i]!=0)
01134       _cell_family_ids[i]->decrRef();
01135   
01136   for (int i=0; i<(int)_face_mesh.size();i++)
01137     if (_face_mesh[i]!=0)
01138       _face_mesh[i]->decrRef();
01139   
01140   for (int i=0; i<(int)_face_family_ids.size();i++)
01141     if (_face_family_ids[i]!=0)
01142       _face_family_ids[i]->decrRef();
01143   
01144   for (std::map<std::string, ParaMEDMEM::DataArrayInt*>::iterator it=_map_dataarray_int.begin() ; it!=_map_dataarray_int.end(); it++ )
01145     if ((*it).second!=0)
01146       (*it).second->decrRef();
01147   
01148   for (std::map<std::string, ParaMEDMEM::DataArrayDouble*>::iterator it=_map_dataarray_double.begin() ; it!=_map_dataarray_double.end(); it++ )
01149     if ((*it).second!=0)
01150       (*it).second->decrRef();
01151   
01152   delete _driver;
01153   if (_topology!=0 && _owns_topology)
01154     delete _topology;
01155 #ifdef HAVE_MPI2
01156   delete _joint_finder;
01157 #endif
01158 }
01159 
01170 void MEDPARTITIONER::MeshCollection::write(const std::string& filename)
01171 {
01172   //building the connect zones necessary for writing joints
01173   //   if (_topology->nbDomain()>1)
01174   //     buildConnectZones();
01175   //suppresses link with driver so that it can be changed for writing
01176   delete _driver;
01177   _driver=0;
01178   retrieveDriver()->write (filename.c_str(), _domain_selector);
01179 }
01180 
01183 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::retrieveDriver()
01184 {
01185   if (_driver==0)
01186     {
01187       switch(_driver_type)
01188         {
01189         case MedXml:
01190           _driver=new MeshCollectionMedXmlDriver(this);
01191           break;
01192         case MedAscii:
01193           _driver=new MeshCollectionMedAsciiDriver(this);
01194           break;
01195         default:
01196           throw INTERP_KERNEL::Exception("Unrecognized driver");
01197         }
01198     }
01199   return _driver;
01200 }
01201 
01202 
01206 MEDPARTITIONER::MeshCollectionDriver* MEDPARTITIONER::MeshCollection::getDriver() const
01207 {
01208   return _driver;
01209 }
01210 
01214 int MEDPARTITIONER::MeshCollection::getMeshDimension() const
01215 {
01216   return _i_non_empty_mesh < 0 ? -1 : _mesh[_i_non_empty_mesh]->getMeshDimension();
01217 }
01218 
01219 int MEDPARTITIONER::MeshCollection::getNbOfLocalMeshes() const
01220 {
01221   int nb=0;
01222   for (int i=0; i<_mesh.size(); i++)
01223     {
01224       if (_mesh[i]) nb++;
01225     }
01226   return nb;
01227 }
01228 
01229 int MEDPARTITIONER::MeshCollection::getNbOfLocalCells() const
01230 {
01231   int nb=0;
01232   for (int i=0; i<_mesh.size(); i++)
01233     {
01234       if (_mesh[i]) nb=nb+_mesh[i]->getNumberOfCells();
01235     }
01236   return nb;
01237 }
01238 
01239 int MEDPARTITIONER::MeshCollection::getNbOfLocalFaces() const
01240 {
01241   int nb=0;
01242   for (int i=0; i<_face_mesh.size(); i++)
01243     {
01244       if (_face_mesh[i]) nb=nb+_face_mesh[i]->getNumberOfCells();
01245     }
01246   return nb;
01247 }
01248 
01249 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getMesh()
01250 {
01251   return _mesh;
01252 }
01253 
01254 std::vector<ParaMEDMEM::MEDCouplingUMesh*>& MEDPARTITIONER::MeshCollection::getFaceMesh()
01255 {
01256   return _face_mesh;
01257 }
01258 
01259 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getMesh(int idomain) const
01260 {
01261   return _mesh[idomain];
01262 }
01263 
01264 ParaMEDMEM::MEDCouplingUMesh* MEDPARTITIONER::MeshCollection::getFaceMesh(int idomain)
01265 {
01266   return _face_mesh[idomain];
01267 }
01268 
01269 std::vector<MEDPARTITIONER::ConnectZone*>& MEDPARTITIONER::MeshCollection::getCZ()
01270 {
01271   return _connect_zones;
01272 }
01273 
01274 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::getTopology() const
01275 {
01276   return _topology;
01277 }
01278 
01279 void MEDPARTITIONER::MeshCollection::setTopology(Topology* topo)
01280 {
01281   if (_topology!=0)
01282     {
01283       throw INTERP_KERNEL::Exception("topology is already set");
01284     }
01285   else
01286     _topology = topo;
01287 }
01288 
01295 void MEDPARTITIONER::MeshCollection::buildCellGraph(MEDPARTITIONER::SkyLineArray* & array, int *& edgeweights )
01296 {
01297   using std::multimap;
01298   using std::vector;
01299   using std::make_pair;
01300   using std::pair;
01301   
01302   std::multimap< int, int > node2cell;
01303   std::multimap< int, int > cell2cell;
01304   std::multimap< int, int > cell2node;
01305 
01306   std::vector<std::vector<std::multimap<int,int> > > commonDistantNodes;
01307   int nbdomain=_topology->nbDomain();
01308 #ifdef HAVE_MPI2
01309   if (isParallelMode())
01310     {
01311       _joint_finder=new JointFinder(*this);
01312       _joint_finder->findCommonDistantNodes();
01313       commonDistantNodes=_joint_finder->getDistantNodeCell();
01314     }
01315 
01316   if (MyGlobals::_Verbose>500)
01317     _joint_finder->print();
01318 #endif
01319   //looking for reverse nodal connectivity i global numbering
01320   for (int idomain=0; idomain<nbdomain; idomain++)
01321     {
01322       if (isParallelMode() && !_domain_selector->isMyDomain(idomain))
01323         continue;
01324     
01325       ParaMEDMEM::DataArrayInt* index=ParaMEDMEM::DataArrayInt::New();
01326       ParaMEDMEM::DataArrayInt* revConn=ParaMEDMEM::DataArrayInt::New();
01327       int nbNodes=_mesh[idomain]->getNumberOfNodes();
01328       _mesh[idomain]->getReverseNodalConnectivity(revConn,index);
01329       //problem saturation over 1 000 000 nodes for 1 proc
01330       if (MyGlobals::_Verbose>100)
01331         std::cout << "proc " << MyGlobals::_Rank << " : getReverseNodalConnectivity done on " << nbNodes << " nodes" << std::endl;
01332       int* index_ptr=index->getPointer();
01333       int* revConnPtr=revConn->getPointer();
01334       for (int i=0; i<nbNodes; i++)
01335         {
01336           for (int icell=index_ptr[i]; icell<index_ptr[i+1]; icell++)
01337             {
01338               int globalNode=_topology->convertNodeToGlobal(idomain,i);
01339               int globalCell=_topology->convertCellToGlobal(idomain,revConnPtr[icell]);
01340               node2cell.insert(make_pair(globalNode, globalCell));
01341               cell2node.insert(make_pair(globalCell, globalNode));
01342             }
01343         }
01344       revConn->decrRef();
01345       index->decrRef();
01346 #ifdef HAVE_MPI2
01347       for (int iother=0; iother<nbdomain; iother++)
01348         {
01349           std::multimap<int,int>::iterator it;
01350           int isource=idomain;
01351           int itarget=iother;
01352           for (it=_joint_finder->getDistantNodeCell()[isource][itarget].begin(); 
01353                it!=_joint_finder->getDistantNodeCell()[isource][itarget].end(); it++)
01354             {
01355               int globalNode=_topology->convertNodeToGlobal(idomain,(*it).first);
01356               int globalCell=(*it).second;
01357               node2cell.insert(make_pair(globalNode, globalCell));
01358               cell2node.insert(make_pair(globalCell, globalNode));
01359             }
01360         }
01361 #endif
01362     }  //endfor idomain
01363 
01364   //creating graph arcs (cell to cell relations)
01365   //arcs are stored in terms of (index,value) notation
01366   // 0 3 5 6 6
01367   // 1 2 3 2 3 3 
01368   // means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
01369   // in present version arcs are not doubled but reflexive (1,1) arcs are present for each cell
01370  
01371   //warning here one node have less than or equal effective number of cell with it
01372   //but cell could have more than effective nodes
01373   //because other equals nodes in other domain (with other global inode)
01374   if (MyGlobals::_Verbose>100)
01375     std::cout<< "proc " << MyGlobals::_Rank << " : creating graph arcs on nbNodes " << _topology->nbNodes() << std::endl;
01376   for (int inode=0; inode<_topology->nbNodes(); inode++)  //on all nodes
01377     {
01378       typedef multimap<int,int>::const_iterator MI;
01379       pair <MI,MI> myRange=node2cell.equal_range(inode);
01380       for (MI cell1=myRange.first; cell1!=myRange.second; cell1++)  //on cells with inode
01381         {
01382           pair <MI,MI> myNodes1=cell2node.equal_range(cell1->second);  //nodes of one cell
01383           for (MI cell2=myRange.first; cell2!=myRange.second; cell2++)  //on one of these cell
01384             {
01385               if (cell1->second!=cell2->second)  //in cells of my domain
01386                 {
01387                   pair <MI,MI> myNodes2=cell2node.equal_range(cell2->second); //on nodes of this cell
01388                   int nbcn=0; //number of common nodes between cells: at least 3 for cells
01389                   for (MI it1=myNodes1.first; it1!=myNodes1.second; ++it1)
01390                     {
01391                       for (MI it2=myNodes2.first; it2!=myNodes2.second; ++it2)
01392                         {
01393                           if ((*it1).second==(*it2).second)
01394                             {
01395                               nbcn=nbcn+1 ; break;
01396                             }
01397                         }
01398                     }
01399                   if (nbcn>=3)  //cvw TODO if 2d cells set at 2
01400                     cell2cell.insert(make_pair(cell1->second,cell2->second));
01401                   //note here there is some global numerotation of cell which come from other domain (not mydomain)
01402                 }
01403             }
01404         }
01405     }
01406   
01407   if (MyGlobals::_Verbose>100)
01408     std::cout << "proc " << MyGlobals::_Rank << " : create skylinearray" << std::endl;
01409   //filling up index and value to create skylinearray structure
01410   std::vector <int> index,value;
01411   index.push_back(0);
01412   int idep=0;
01413   
01414   for (int idomain=0; idomain<nbdomain; idomain++)
01415     {
01416       if (isParallelMode() && !_domain_selector->isMyDomain(idomain)) continue;
01417       int nbCells=_mesh[idomain]->getNumberOfCells();
01418       for (int icell=0; icell<nbCells; icell++)
01419         {
01420           int size=0;
01421           int globalCell=_topology->convertCellToGlobal(idomain,icell);
01422           multimap<int,int>::iterator it;
01423           pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
01424           ret=cell2cell.equal_range(globalCell);
01425           for (it=ret.first; it!=ret.second; ++it)
01426             {
01427               int ival=(*it).second; //no adding one existing yet
01428               for (int i=idep ; i<idep+size ; i++)
01429                 {
01430                   if (value[i]==ival)
01431                     {
01432                       ival= -1; break;
01433                     }
01434                 }
01435               if (ival!= -1)
01436                 {
01437                   value.push_back(ival);
01438                   size++;
01439                 }
01440             }
01441           idep=index[index.size()-1]+size;
01442           index.push_back(idep);
01443         }
01444     }
01445   
01446   array=new MEDPARTITIONER::SkyLineArray(index,value);
01447 
01448   if (MyGlobals::_Verbose>100)
01449     {
01450       std::cout << "\nproc " << _domain_selector->rank() << " : end MeshCollection::buildCellGraph " <<
01451         index.size()-1 << " " << value.size() << std::endl;
01452       if (index.size()>1)
01453         {
01454           for (int i=0; i<10; ++i)
01455             std::cout<<index[i]<<" ";
01456           std::cout << "... " << index[index.size()-1] << std::endl;
01457           for (int i=0; i<15; ++i)
01458             std::cout<< value[i] << " ";
01459           int ll=index[index.size()-1]-1;
01460           std::cout << "... (" << ll << ") " << value[ll-1] << " " << value[ll] << std::endl;
01461         }
01462     }
01463   
01464 }
01465 
01466 
01473 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(int nbdomain,
01474                                                                           Graph::splitter_type split, 
01475                                                                           const std::string& options_string,
01476                                                                           int *user_edge_weights,
01477                                                                           int *user_vertices_weights)
01478 {
01479   if (MyGlobals::_Verbose>10)
01480     std::cout << "proc " << MyGlobals::_Rank << " : MeshCollection::createPartition : Building cell graph" << std::endl;
01481   
01482   if (nbdomain <1)
01483     throw INTERP_KERNEL::Exception("Number of subdomains must be > 0");
01484   MEDPARTITIONER::SkyLineArray* array=0;
01485   int* edgeweights=0;
01486   buildCellGraph(array,edgeweights);
01487   
01488   Graph* cellGraph;
01489   switch (split)
01490     {
01491     case Graph::METIS:
01492 #if defined(MED_ENABLE_PARMETIS) || defined(MED_ENABLE_METIS)
01493       if (MyGlobals::_Verbose>10)
01494         std::cout << "METISGraph" << std::endl;
01495       cellGraph=new METISGraph(array,edgeweights);
01496 #else
01497       throw INTERP_KERNEL::Exception("MeshCollection::createPartition : PARMETIS/METIS is not available. Check your products, please.");
01498 #endif
01499       break;
01500     case Graph::SCOTCH:
01501 #ifdef MED_ENABLE_SCOTCH
01502       if (MyGlobals::_Verbose>10)
01503         std::cout << "SCOTCHGraph" << std::endl;
01504       cellGraph=new SCOTCHGraph(array,edgeweights);
01505 #else
01506       throw INTERP_KERNEL::Exception("MeshCollection::createPartition : SCOTCH is not available. Check your products, please.");
01507 #endif
01508       break;
01509     }
01510 
01512   if (user_edge_weights!=0) 
01513     cellGraph->setEdgesWeights(user_edge_weights);
01514   if (user_vertices_weights!=0)
01515     cellGraph->setVerticesWeights(user_vertices_weights);
01516 
01517   if (MyGlobals::_Is0verbose>10)
01518     std::cout << "partitioning graph on " << nbdomain << " domains" << std::endl;
01519   cellGraph->partGraph(nbdomain, options_string, _domain_selector);
01520 
01521   if (MyGlobals::_Is0verbose>10)
01522     std::cout << "building new topology" << std::endl;
01523   //cellGraph is a shared pointer 
01524   Topology *topology=0;
01525   topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
01526   //cleaning
01527   delete [] edgeweights;
01528   delete cellGraph;
01529   if (MyGlobals::_Verbose>11)
01530     std::cout << "proc " << MyGlobals::_Rank << " : end MeshCollection::createPartition" << std::endl;
01531   return topology;
01532 }
01533 
01540 MEDPARTITIONER::Topology* MEDPARTITIONER::MeshCollection::createPartition(const int* partition)
01541 {
01542   MEDPARTITIONER::SkyLineArray* array=0;
01543   int* edgeweights=0;
01544 
01545   buildCellGraph(array,edgeweights);
01546   Graph* cellGraph;
01547   std::set<int> domains;
01548   for (int i=0; i<_topology->nbCells(); i++)
01549     {
01550       domains.insert(partition[i]);
01551     }
01552   cellGraph=new UserGraph(array, partition, _topology->nbCells());
01553   
01554   //cellGraph is a shared pointer 
01555   Topology *topology=0;
01556   int nbdomain=domains.size();
01557   topology=new ParallelTopology (cellGraph, getTopology(), nbdomain, getMeshDimension());
01558   //  if (array!=0) delete array;
01559   delete cellGraph;
01560   return topology;
01561 }
01562 
01563 void MEDPARTITIONER::MeshCollection::setDomainNames(const std::string& name)
01564 {
01565   for (int i=0; i<_topology->nbDomain(); i++)
01566     {
01567       std::ostringstream oss;
01568       oss<<name<<"_"<<i;
01569       if (!isParallelMode() || _domain_selector->isMyDomain(i))
01570         _mesh[i]->setName(oss.str().c_str());
01571     }
01572 }
01573 
01574 ParaMEDMEM::DataArrayDouble *MEDPARTITIONER::MeshCollection::getField(std::string descriptionField, int iold)
01575 //getField look for and read it if not done, and assume decrRef() in ~MeshCollection;
01576 //something like MEDCouplingFieldDouble *f2=MEDLoader::ReadFieldCell(name.c_str(),f1->getMesh()->getName(),0,f1->getName(),0,1);
01577 {
01578   int rank=MyGlobals::_Rank;
01579   std::string tag="ioldFieldDouble="+IntToStr(iold);
01580   std::string descriptionIold=descriptionField+SerializeFromString(tag);
01581   if (_map_dataarray_double.find(descriptionIold)!=_map_dataarray_double.end())
01582     {
01583       if (MyGlobals::_Verbose>300)
01584         std::cout << "proc " << rank << " : YET READ getField : " << descriptionIold << std::endl;
01585       ParaMEDMEM::DataArrayDouble* res=_map_dataarray_double[descriptionIold];
01586       return res;
01587     }
01588   if (MyGlobals::_Verbose>200)
01589     std::cout << "proc " << rank << " : TO BE READ getField : " << descriptionIold << std::endl;
01590   std::string description, fileName, meshName, fieldName;
01591   int typeField, DT, IT, entity;
01592   fileName=MyGlobals::_File_Names[iold];
01593   if (MyGlobals::_Verbose>10) 
01594     std::cout << "proc " << MyGlobals::_Rank << " : in " << fileName << " " << iold << " " << descriptionIold << std::endl;
01595   FieldShortDescriptionToData(descriptionIold, fieldName, typeField, entity, DT, IT);
01596   meshName=MyGlobals::_Mesh_Names[iold];
01597   
01598   ParaMEDMEM::MEDCouplingFieldDouble* f2=MEDLoader::ReadField((ParaMEDMEM::TypeOfField) typeField,
01599                                                               fileName.c_str(), meshName.c_str(), 0, fieldName.c_str(), DT, IT);
01600   
01601   ParaMEDMEM::DataArrayDouble* res=f2->getArray();
01602   //to know names of components
01603   std::vector<std::string> browse=BrowseFieldDouble(f2);
01604   std::string localFieldInformation=descriptionIold+SerializeFromVectorOfString(browse);
01605   if (MyGlobals::_Verbose>10)
01606     std::cout << "proc " << MyGlobals::_Rank << " : localFieldInformation : " << localFieldInformation << std::endl;
01607   MyGlobals::_General_Informations.push_back(localFieldInformation);
01608   res->incrRef();
01609   f2->decrRef();
01610   _map_dataarray_double[descriptionIold]=res; 
01611   return res;
01612 }
01613 
01614 void MEDPARTITIONER::MeshCollection::prepareFieldDescriptions()
01615 //to have unique valid fields names/pointers/descriptions for partitionning
01616 //filter _field_descriptions to be in all procs compliant and equal
01617 {
01618   int nbfiles=MyGlobals::_File_Names.size(); //nb domains
01619   std::vector<std::string> r2;
01620   //from allgatherv then vector(procs) of serialised vector(fields) of vector(description) data
01621   for (int i=0; i<(int)_field_descriptions.size(); i++)
01622     {
01623       std::vector<std::string> r1=DeserializeToVectorOfString(_field_descriptions[i]);
01624       for (int ii=0; ii<(int)r1.size(); ii++)
01625         r2.push_back(r1[ii]);
01626     }
01627   //here vector(procs*fields) of serialised vector(description) data
01628   _field_descriptions=r2;
01629   int nbfields=_field_descriptions.size(); //on all domains
01630   if ((nbfields%nbfiles)!=0)
01631     {
01632       if (MyGlobals::_Rank==0)
01633         {
01634           std::cerr<< "\nERROR : incoherent number of fields references in all files .med\n" << std::endl
01635               << "fileMedNames :" << std::endl
01636               << ReprVectorOfString(MyGlobals::_File_Names)
01637               << "field_descriptions :" << std::endl
01638               << ReprVectorOfString(MyGlobals::_Field_Descriptions);
01639         }
01640       throw INTERP_KERNEL::Exception("incoherent number of fields references in all files .med\n");
01641     }
01642   _field_descriptions.resize(nbfields/nbfiles);
01643   for (int i=0; i<(int)_field_descriptions.size(); i++)
01644     {
01645       std::string str=_field_descriptions[i];
01646       str=EraseTagSerialized(str,"idomain=");
01647       str=EraseTagSerialized(str,"fileName=");
01648       _field_descriptions[i]=str;
01649     }
01650 }
01651 
01652 //returns true if inodes of a face are in inodes of a cell
01653 bool isFaceOncell(std::vector< int >& inodesFace, std::vector< int >&  inodesCell)
01654 {
01655   int ires=0;
01656   int nbok=inodesFace.size();
01657   for (int i=0; i<nbok; i++)
01658     {
01659       int ii=inodesFace[i];
01660       if (ii<0)
01661         std::cout << "isFaceOncell problem inodeface<0" << std::endl;
01662       for (int j=0; j<(int)inodesCell.size(); j++)
01663         {
01664           if (ii==inodesCell[j])
01665             {
01666               ires=ires+1;
01667               break; //inode of face found
01668             }
01669         }
01670       if (ires<i+1)
01671         break; //inode of face not found do not continue...
01672     }
01673   return (ires==nbok);
01674 }
01675 
01676 void MEDPARTITIONER::MeshCollection::filterFaceOnCell()
01677 {
01678   for (int inew=0; inew<_topology->nbDomain(); inew++)
01679     {
01680       if (isParallelMode() && _domain_selector->isMyDomain(inew))
01681         {
01682           if (MyGlobals::_Verbose>200) 
01683             std::cout << "proc " << MyGlobals::_Rank << " : filterFaceOnCell on inewDomain " << inew << " nbOfFaces " << _face_mesh[inew]->getNumberOfCells() << std::endl;
01684           ParaMEDMEM::MEDCouplingUMesh* mcel=_mesh[inew];
01685           ParaMEDMEM::MEDCouplingUMesh* mfac=_face_mesh[inew];
01686       
01687           //to have cellnode=f(facenode)... inodeCell=nodeIds[inodeFace]
01688           std::vector<int> nodeIds;
01689           getNodeIds(*mcel, *mfac, nodeIds);
01690           if (nodeIds.size()==0)
01691             continue;  //one empty mesh nothing to do
01692       
01693           ParaMEDMEM::DataArrayInt *revNodalCel=ParaMEDMEM::DataArrayInt::New();
01694           ParaMEDMEM::DataArrayInt *revNodalIndxCel=ParaMEDMEM::DataArrayInt::New();
01695           mcel->getReverseNodalConnectivity(revNodalCel,revNodalIndxCel);
01696           int *revC=revNodalCel->getPointer();
01697           int *revIndxC=revNodalIndxCel->getPointer();
01698       
01699           std::vector< int > faceOnCell;
01700           std::vector< int > faceNotOnCell;
01701           int nbface=mfac->getNumberOfCells();
01702           for (int iface=0; iface<nbface; iface++)
01703             {
01704               bool ok;
01705               std::vector< int > inodesFace;
01706               mfac->getNodeIdsOfCell(iface, inodesFace);
01707               int nbnodFace=inodesFace.size();
01708               //set inodesFace in mcel
01709               for (int i=0; i<nbnodFace; i++) inodesFace[i]=nodeIds[inodesFace[i]];
01710               int inod=inodesFace[0];
01711               if (inod<0)
01712                 std::cout << "filterFaceOnCell problem 1" << std::endl;
01713               int nbcell=revIndxC[inod+1]-revIndxC[inod];
01714               for (int j=0; j<nbcell; j++) //look for each cell with inod
01715                 {
01716                   int icel=revC[revIndxC[inod]+j];
01717                   std::vector< int > inodesCell;
01718                   mcel->getNodeIdsOfCell(icel, inodesCell);
01719                   ok=isFaceOncell(inodesFace, inodesCell);
01720                   if (ok) break;
01721                 }
01722               if (ok)
01723                 {
01724                   faceOnCell.push_back(iface);
01725                 }
01726               else
01727                 {
01728                   faceNotOnCell.push_back(iface);
01729                   if (MyGlobals::_Is0verbose>300)
01730                     std::cout << "face NOT on cell " << iface << " " << faceOnCell.size()-1 << std::endl;
01731                 }
01732             }
01733       
01734           revNodalCel->decrRef();
01735           revNodalIndxCel->decrRef();
01736       
01737           std::string keyy;
01738           keyy=Cle1ToStr("filterFaceOnCell",inew);
01739           _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceOnCell);
01740           keyy=Cle1ToStr("filterNotFaceOnCell",inew);
01741           _map_dataarray_int[keyy]=CreateDataArrayIntFromVector(faceNotOnCell);
01742         }
01743     }
01744 }