src/MEDCoupling/MEDCouplingVoronoi.cxx

   1 // Copyright (C) 2007-2017  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (EDF R&D)
  20
  21 #include "MEDCouplingVoronoi.hxx"
  22 #include "MEDCoupling1GTUMesh.hxx"
  23 #include "MEDCouplingCMesh.hxx"
  24 #include "MCAuto.txx"
  25
  26 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
  27 #include "Interpolation2D.txx"
  28 #include "Interpolation3DSurf.hxx"
  29
  30 using namespace MEDCoupling;
  31
  32 Voronizer::~Voronizer()
  33 {
  34 }
  35
  36 int Voronizer2D::getDimension() const
  37 {
  38   return 2;
  39 }
  40
  41 int Voronizer3D::getDimension() const
  42 {
  43   return 3;
  44 }
  45
  46 MCAuto<MEDCouplingUMesh> ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector<double>& bbox, double eps)
  47 {
  48   static const double FACT=1.2;
  49   MCAuto<MEDCouplingCMesh> m(MEDCouplingCMesh::New());
  50   MCAuto<DataArrayDouble> arr1(DataArrayDouble::New()); arr1->alloc(2,1) ; arr1->setIJ(0,0,bbox[0]); arr1->setIJ(1,0,bbox[1]);
  51   MCAuto<DataArrayDouble> arr2(DataArrayDouble::New()); arr2->alloc(2,1) ; arr2->setIJ(0,0,bbox[2]); arr2->setIJ(1,0,bbox[3]);
  52   m->setCoords(arr1,arr2);
  53   static const double PT[2]={0.,0.};
  54   m->scale(PT,FACT);
  55   MCAuto<MEDCouplingUMesh> mu(m->buildUnstructured());
  56   double l(std::max(bbox[1]-bbox[0],bbox[3]-bbox[2]));
  57   double middle[2]={(pt1[0]+pt2[0])/2.,(pt1[1]+pt2[1])/2.};
  58   double v[2]={pt1[0],pt1[1]};
  59   DataArrayDouble::Rotate2DAlg(middle,M_PI/2,1,v,v);
  60   v[0]=middle[0]-v[0]; v[1]=middle[1]-v[1];
  61   {
  62     double nor(sqrt(v[0]*v[0]+v[1]*v[1]));
  63     v[0]/=nor; v[1]/=nor;
  64   }
  65   MCAuto<MEDCouplingUMesh> line(MEDCouplingUMesh::New("line",1));
  66   {
  67     MCAuto<DataArrayDouble> coo(DataArrayDouble::New()); coo->alloc(2,2);
  68     coo->setIJ(0,0,middle[0]-2.*l*v[0]); coo->setIJ(0,1,middle[1]-2.*l*v[1]); coo->setIJ(1,0,middle[0]+2.*l*v[0]); coo->setIJ(1,1,middle[1]+2.*l*v[1]);
  69     line->setCoords(coo);
  70   }
  71   line->allocateCells();
  72   static const int CONN[2]={0,1};
  73   line->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN);
  74   MCAuto<MEDCouplingUMesh> sp2,sp1;
  75   {
  76     DataArrayInt *cellNb1(0),*cellNb2(0);
  77     MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0);
  78     MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2);
  79     sp1=sp1Pt; sp2=sp2Pt;
  80     MCAuto<DataArrayInt> cellNb10(cellNb1),cellNb20(cellNb2);
  81   }
  82   std::vector<int> ccp;
  83   sp2->getCellsContainingPoint(pt1,eps,ccp);
  84   if(ccp.size()!=1)
  85     throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !");
  86   MCAuto<MEDCouplingUMesh> ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true));
  87   ret->zipCoords();
  88   return ret;
  89 }
  90
  91
  92 MCAuto<MEDCouplingUMesh> MergeVorCells2D(MEDCouplingUMesh *p, double eps, bool isZip)
  93 {
  94   MCAuto<DataArrayInt> edgeToKeep;
  95   MCAuto<MEDCouplingUMesh> p0;
  96   {
  97     MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
  98     p0=p->buildDescendingConnectivity(d,di,rd,rdi);
  99     MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
 100     edgeToKeep=dsi->findIdsEqual(1);
 101   }
 102   MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
 103   if(isZip)
 104     {
 105       skinOfRes->zipCoords();
 106       if(skinOfRes->getNumberOfCells()!=skinOfRes->getNumberOfNodes())
 107         throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !");
 108     }
 109   MCAuto<DataArrayInt> d(skinOfRes->orderConsecutiveCells1D());
 110   MCAuto<MEDCoupling1SGTUMesh> skinOfRes2;
 111   {
 112     MCAuto<MEDCouplingUMesh> part(skinOfRes->buildPartOfMySelf(d->begin(),d->end()));
 113     skinOfRes2=MEDCoupling1SGTUMesh::New(part);
 114   }
 115   MCAuto<DataArrayInt> c(skinOfRes2->getNodalConnectivity()->deepCopy());
 116   c->circularPermutation(1);
 117   c->rearrange(2);
 118   std::vector< MCAuto<DataArrayInt> > vdi(c->explodeComponents());
 119   if(!vdi[0]->isEqual(*vdi[1]))
 120     throw INTERP_KERNEL::Exception("MergeVorCells : internal error !");
 121   MCAuto<MEDCouplingUMesh> m(MEDCouplingUMesh::New("",2));
 122   m->setCoords(skinOfRes2->getCoords());
 123   m->allocateCells();
 124   m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin());
 125   return m;
 126 }
 127
 128 MCAuto<MEDCouplingUMesh> MergeVorCells(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
 129 {
 130   std::size_t sz(vcs.size());
 131   if(sz<1)
 132     throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
 133   if(sz==1)
 134     return vcs[0];
 135   MCAuto<MEDCouplingUMesh> p;
 136   {
 137     std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
 138     p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
 139   }
 140   p->zipCoords();
 141   {
 142     bool dummy; int dummy2;
 143     MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
 144   }
 145   return MergeVorCells2D(p,eps,true);
 146 }
 147
 148 MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
 149 {
 150   std::size_t sz(vcs.size());
 151   if(sz<1)
 152     throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
 153   if(sz==1)
 154     return vcs[0];
 155   MCAuto<MEDCouplingUMesh> p;
 156   {
 157     std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
 158     p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
 159   }
 160   p->zipCoords();
 161   {
 162     bool dummy; int dummy2;
 163     MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
 164   }
 165   MCAuto<DataArrayInt> edgeToKeep;
 166   MCAuto<MEDCouplingUMesh> p0;
 167   {
 168     MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
 169     p0=p->buildDescendingConnectivity(d,di,rd,rdi);
 170     MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
 171     edgeToKeep=dsi->findIdsEqual(1);
 172   }
 173   MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
 174   MCAuto<DataArrayDouble> eqn(skinOfRes->computePlaneEquationOf3DFaces());
 175   MCAuto<DataArrayInt> comm,commI;
 176   {
 177     DataArrayInt *a(0),*b(0);
 178     eqn->findCommonTuples(eps,0,a,b);
 179     comm=a; commI=b;
 180     //comm=DataArrayInt::New(); comm->alloc(0,1); commI=DataArrayInt::New(); commI->alloc(1,1); commI->setIJ(0,0,0);
 181   }
 182   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",3));
 183   ret->setCoords(skinOfRes->getCoords());
 184   ret->allocateCells();
 185   std::vector<int> conn;
 186   int jj(0);
 187   for(int i=0;i<commI->getNumberOfTuples()-1;i++,jj++)
 188     {
 189       if(jj!=0)
 190         conn.push_back(-1);
 191       MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(comm->begin()+commI->getIJ(i,0),comm->begin()+commI->getIJ(i+1,0),true));
 192       MCAuto<MEDCouplingUMesh> tmp2;
 193       if(commI->getIJ(i+1,0)-commI->getIJ(i,0)==1)
 194         tmp2=tmp;
 195       else
 196         tmp2=MergeVorCells2D(tmp,eps,false);
 197       const int *cPtr(tmp2->getNodalConnectivity()->begin()),*ciPtr(tmp2->getNodalConnectivityIndex()->begin());
 198       conn.insert(conn.end(),cPtr+1,cPtr+ciPtr[1]);
 199     }
 200   MCAuto<DataArrayInt> remain(comm->buildComplement(skinOfRes->getNumberOfCells()));
 201   {
 202     MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(remain->begin(),remain->end(),true));
 203     const int *cPtr(tmp->getNodalConnectivity()->begin()),*ciPtr(tmp->getNodalConnectivityIndex()->begin());
 204     for(int i=0;i<remain->getNumberOfTuples();i++,jj++)
 205       {
 206         if(jj!=0)
 207           conn.push_back(-1);
 208         conn.insert(conn.end(),cPtr+ciPtr[i]+1,cPtr+ciPtr[i+1]);
 209       }
 210   }
 211   ret->insertNextCell(INTERP_KERNEL::NORM_POLYHED,conn.size(),&conn[0]);
 212   return ret;
 213 }
 214
 215 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer2D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 216 {
 217   if(!m || !points)
 218     throw INTERP_KERNEL::Exception("Voronoize2D : null pointer !");
 219   m->checkConsistencyLight();
 220   points->checkAllocated();
 221   if(m->getMeshDimension()!=2 || m->getSpaceDimension()!=2 || points->getNumberOfComponents()!=2)
 222     throw INTERP_KERNEL::Exception("Voronoize2D : spacedim must be equal to 2 and meshdim also equal to 2 !");
 223   if(m->getNumberOfCells()!=1)
 224     throw INTERP_KERNEL::Exception("Voronoize2D : mesh is expected to have only one cell !");
 225   int nbPts(points->getNumberOfTuples());
 226   if(nbPts<1)
 227     throw INTERP_KERNEL::Exception("Voronoize2D : at least one point expected !");
 228   std::vector<double> bbox(4);
 229   m->getBoundingBox(&bbox[0]);
 230   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 231   const double *pts(points->begin());
 232   for(int i=1;i<nbPts;i++)
 233     {
 234       MCAuto<MEDCouplingUMesh> vorTess;
 235       {
 236         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 237         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 238       }
 239       {
 240         bool dummy;
 241         int newNbNodes;
 242         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 243       }
 244       std::vector<int> polygsToIterOn;
 245       const double *pt(pts+i*2);
 246       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 247       if(polygsToIterOn.size()<1)
 248         throw INTERP_KERNEL::Exception("Voronoize2D : presence of a point outside the given cell !");
 249       std::set<int> elemsToDo,elemsDone; elemsToDo.insert(polygsToIterOn[0]);
 250       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 251       while(!elemsToDo.empty())
 252         {
 253           int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
 254           const double *seed(pts+2*poly);
 255           MCAuto<MEDCouplingUMesh> cell(ComputeBigCellFrom(pt,seed,bbox,eps));
 256           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 257           tile->zipCoords();
 258           MCAuto<MEDCouplingUMesh> a;
 259           MCAuto<DataArrayInt> b,c;
 260           {
 261             DataArrayInt *bPtr(0),*cPtr(0);
 262             a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
 263             b=bPtr; c=cPtr;
 264           }
 265           MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
 266           if(part->getNumberOfTuples()!=1)
 267             throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
 268           MCAuto<MEDCouplingUMesh> newVorCell;
 269           {
 270             MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
 271             newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
 272           }
 273           newVorCell->zipCoords();
 274           MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
 275           modifiedCell->zipCoords();
 276           l0[poly]=modifiedCell;
 277           //
 278           MCAuto<DataArrayInt> ids;
 279           {
 280             DataArrayInt *tmp(0);
 281             bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
 282             ids=tmp;
 283             if(!sta)
 284               throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
 285           }
 286           MCAuto<DataArrayDouble> newCoords;
 287           {
 288             MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
 289             newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
 290           }
 291           const double *cPtr(newCoords->begin());
 292           for(int i=0;i<newCoords->getNumberOfTuples();i++,cPtr+=2)
 293             {
 294               std::set<int> zeCandidates;
 295               {
 296                 std::vector<int> zeCandidatesTmp;
 297                 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
 298                 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
 299               }
 300               std::set<int> tmp,newElementsToDo;
 301               std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
 302               std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
 303               elemsToDo=newElementsToDo;
 304             }
 305           newVorCells.push_back(newVorCell);
 306         }
 307       l0.push_back(MergeVorCells(newVorCells,eps));
 308     }
 309   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 310   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 311   {
 312     bool dummy; int dummy2;
 313     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 314   }
 315   return ret;
 316 }
 317
 318 MCAuto<MEDCouplingUMesh> Split3DCellInParts(const MEDCouplingUMesh *m, const double pt[3], const double seed[3], double eps, int tmp[2])
 319 {
 320   if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || m->getNumberOfCells()!=1)
 321     throw INTERP_KERNEL::Exception("Split3DCellInParts : expecting a 3D with exactly one cell !");
 322   double middle[3]={(pt[0]+seed[0])/2.,(pt[1]+seed[1])/2.,(pt[2]+seed[2])/2.};
 323   double vec[3]={pt[0]-seed[0],pt[1]-seed[1],pt[2]-seed[2]};
 324   MCAuto<MEDCouplingUMesh> res(m->clipSingle3DCellByPlane(middle,vec,eps));
 325   return res;
 326 }
 327
 328 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 329 {
 330   double eps2(1.-sqrt(eps));// 2nd eps for interpolation. Here the eps is computed to feet cos(eps) ~ 1-eps^2
 331   if(!m || !points)
 332     throw INTERP_KERNEL::Exception("Voronoize3D : null pointer !");
 333   m->checkConsistencyLight();
 334   points->checkAllocated();
 335   if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || points->getNumberOfComponents()!=3)
 336     throw INTERP_KERNEL::Exception("Voronoize3D : spacedim must be equal to 3 and meshdim also equal to 3 !");
 337   if(m->getNumberOfCells()!=1)
 338     throw INTERP_KERNEL::Exception("Voronoize3D : mesh is expected to have only one cell !");
 339   int nbPts(points->getNumberOfTuples());
 340   if(nbPts<1)
 341     throw INTERP_KERNEL::Exception("Voronoize3D : at least one point expected !");
 342   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 343   const double *pts(points->begin());
 344   for(int i=1;i<nbPts;i++)
 345     {
 346       MCAuto<MEDCouplingUMesh> vorTess;
 347       {
 348         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 349         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 350       }
 351       {
 352         bool dummy;
 353         int newNbNodes;
 354         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 355       }
 356       std::vector<int> polygsToIterOn;
 357       const double *pt(pts+i*3);
 358       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 359       if(polygsToIterOn.size()<1)
 360         throw INTERP_KERNEL::Exception("Voronoize3D : presence of a point outside the given cell !");
 361       std::set<int> elemsToDo(polygsToIterOn.begin(),polygsToIterOn.end()),elemsDone;
 362       std::size_t ii(0);
 363       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 364       MCAuto<DataArrayInt> d(DataArrayInt::New()),dI(DataArrayInt::New()),rd(DataArrayInt::New()),rdI(DataArrayInt::New());
 365       MCAuto<MEDCouplingUMesh> faces(vorTess->buildDescendingConnectivity(d,dI,rd,rdI));
 366       //
 367       while(!elemsToDo.empty())
 368         {
 369           int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
 370           const double *seed(pts+3*poly);
 371           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 372           tile->zipCoords();
 373           int tmp[2];
 374           MCAuto<MEDCouplingUMesh> cells;
 375           try
 376             {
 377               cells=Split3DCellInParts(tile,pt,seed,eps,tmp);
 378             }
 379           catch(INTERP_KERNEL::Exception& e)
 380             {
 381               continue;
 382             }
 383           MCAuto<MEDCouplingUMesh> newVorCell(cells->buildPartOfMySelfSlice(1,2,1,true));
 384           newVorCell->zipCoords();
 385           MCAuto<MEDCouplingUMesh> modifiedCell(cells->buildPartOfMySelfSlice(0,1,1,true));
 386           modifiedCell->zipCoords();
 387           l0[poly]=modifiedCell;
 388           if(std::find(polygsToIterOn.begin(),polygsToIterOn.end(),poly)!=polygsToIterOn.end())// we iterate on a polyhedron containg the point to add pt -> add cells sharing faces with just computed newVorCell
 389             {
 390               MCAuto<MEDCouplingUMesh> faces2;
 391               {
 392                 MCAuto<DataArrayInt> d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New());
 393                 faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I);
 394               }
 395               MCAuto<MEDCouplingUMesh> faces3(faces2->buildPartOfMySelfSlice(1,faces2->getNumberOfCells(),1,true));// suppress internal face
 396               MCAuto<MEDCouplingUMesh> facesOfCurSplitPol(faces->buildPartOfMySelf(d->begin()+dI->getIJ(poly,0),d->begin()+dI->getIJ(poly+1,0),true));
 397               // intersection between the out faces of newVorCell and the neighbor faces of poly polyhedron -> candidates
 398               MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(facesOfCurSplitPol);
 399               MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(faces3);
 400               INTERP_KERNEL::Interpolation3DSurf interpolation;
 401               interpolation.setMinDotBtwPlane3DSurfIntersect(eps2);
 402               interpolation.setMaxDistance3DSurfIntersect(eps);
 403               interpolation.setPrecision(1e-12);
 404               std::vector<std::map<int,double> > matrix;
 405               interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix,"P0P0");
 406               std::set<int> zeCandidates;
 407               for(std::vector<std::map<int,double> >::const_iterator it2=matrix.begin();it2!=matrix.end();it2++)
 408                 for(std::map<int,double>::const_iterator it3=(*it2).begin();it3!=(*it2).end();it3++)
 409                   {
 410                     int faceIdInVorTess(d->getIJ(dI->getIJ(poly,0)+(*it3).first,0));
 411                     for(const int *it4=rd->begin()+rdI->getIJ(faceIdInVorTess,0);it4!=rd->begin()+rdI->getIJ(faceIdInVorTess+1,0);it4++)
 412                       {
 413                         if(*it4!=poly)
 414                           zeCandidates.insert(*it4);
 415                       }
 416                   }
 417               std::set<int> tmp,newElementsToDo;
 418               std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
 419               std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
 420               elemsToDo=newElementsToDo;
 421             }
 422           //
 423           newVorCells.push_back(newVorCell);
 424           ii++;
 425         }
 426       l0.push_back(MergeVorCells3D(newVorCells,eps));
 427     }
 428   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 429   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 430   {
 431     bool dummy; int dummy2;
 432     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 433   }
 434   return ret;
 435 }