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 Voronizer1D::getDimension() const
  37 {
  38   return 1;
  39 }
  40
  41 int Voronizer2D::getDimension() const
  42 {
  43   return 2;
  44 }
  45
  46 int Voronizer3D::getDimension() const
  47 {
  48   return 3;
  49 }
  50
  51 MCAuto<MEDCouplingUMesh> ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector<double>& bbox, double eps)
  52 {
  53   static const double FACT=1.2;
  54   MCAuto<MEDCouplingCMesh> m(MEDCouplingCMesh::New());
  55   MCAuto<DataArrayDouble> arr1(DataArrayDouble::New()); arr1->alloc(2,1) ; arr1->setIJ(0,0,bbox[0]); arr1->setIJ(1,0,bbox[1]);
  56   MCAuto<DataArrayDouble> arr2(DataArrayDouble::New()); arr2->alloc(2,1) ; arr2->setIJ(0,0,bbox[2]); arr2->setIJ(1,0,bbox[3]);
  57   m->setCoords(arr1,arr2);
  58   static const double PT[2]={0.,0.};
  59   m->scale(PT,FACT);
  60   MCAuto<MEDCouplingUMesh> mu(m->buildUnstructured());
  61   double l(std::max(bbox[1]-bbox[0],bbox[3]-bbox[2]));
  62   double middle[2]={(pt1[0]+pt2[0])/2.,(pt1[1]+pt2[1])/2.};
  63   double v[2]={pt1[0],pt1[1]};
  64   DataArrayDouble::Rotate2DAlg(middle,M_PI/2,1,v,v);
  65   v[0]=middle[0]-v[0]; v[1]=middle[1]-v[1];
  66   {
  67     double nor(sqrt(v[0]*v[0]+v[1]*v[1]));
  68     v[0]/=nor; v[1]/=nor;
  69   }
  70   MCAuto<MEDCouplingUMesh> line(MEDCouplingUMesh::New("line",1));
  71   {
  72     MCAuto<DataArrayDouble> coo(DataArrayDouble::New()); coo->alloc(2,2);
  73     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]);
  74     line->setCoords(coo);
  75   }
  76   line->allocateCells();
  77   static const int CONN[2]={0,1};
  78   line->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN);
  79   MCAuto<MEDCouplingUMesh> sp2,sp1;
  80   {
  81     DataArrayInt *cellNb1(0),*cellNb2(0);
  82     MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0);
  83     MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2);
  84     sp1=sp1Pt; sp2=sp2Pt;
  85     MCAuto<DataArrayInt> cellNb10(cellNb1),cellNb20(cellNb2);
  86   }
  87   std::vector<int> ccp;
  88   sp2->getCellsContainingPoint(pt1,eps,ccp);
  89   if(ccp.size()!=1)
  90     throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !");
  91   MCAuto<MEDCouplingUMesh> ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true));
  92   ret->zipCoords();
  93   return ret;
  94 }
  95
  96
  97 MCAuto<MEDCouplingUMesh> MergeVorCells2D(MEDCouplingUMesh *p, double eps, bool isZip)
  98 {
  99   MCAuto<DataArrayInt> edgeToKeep;
 100   MCAuto<MEDCouplingUMesh> p0;
 101   {
 102     MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
 103     p0=p->buildDescendingConnectivity(d,di,rd,rdi);
 104     MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
 105     edgeToKeep=dsi->findIdsEqual(1);
 106   }
 107   MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
 108   if(isZip)
 109     {
 110       skinOfRes->zipCoords();
 111       if(skinOfRes->getNumberOfCells()!=skinOfRes->getNumberOfNodes())
 112         throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !");
 113     }
 114   MCAuto<DataArrayInt> d(skinOfRes->orderConsecutiveCells1D());
 115   MCAuto<MEDCoupling1SGTUMesh> skinOfRes2;
 116   {
 117     MCAuto<MEDCouplingUMesh> part(skinOfRes->buildPartOfMySelf(d->begin(),d->end()));
 118     skinOfRes2=MEDCoupling1SGTUMesh::New(part);
 119   }
 120   MCAuto<DataArrayInt> c(skinOfRes2->getNodalConnectivity()->deepCopy());
 121   c->circularPermutation(1);
 122   c->rearrange(2);
 123   std::vector< MCAuto<DataArrayInt> > vdi(c->explodeComponents());
 124   if(!vdi[0]->isEqual(*vdi[1]))
 125     throw INTERP_KERNEL::Exception("MergeVorCells : internal error !");
 126   MCAuto<MEDCouplingUMesh> m(MEDCouplingUMesh::New("",2));
 127   m->setCoords(skinOfRes2->getCoords());
 128   m->allocateCells();
 129   m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin());
 130   return m;
 131 }
 132
 133 MCAuto<MEDCouplingUMesh> MergeVorCells(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
 134 {
 135   std::size_t sz(vcs.size());
 136   if(sz<1)
 137     throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
 138   if(sz==1)
 139     return vcs[0];
 140   MCAuto<MEDCouplingUMesh> p;
 141   {
 142     std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
 143     p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
 144   }
 145   p->zipCoords();
 146   {
 147     bool dummy; int dummy2;
 148     MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
 149   }
 150   return MergeVorCells2D(p,eps,true);
 151 }
 152
 153 MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
 154 {
 155   std::size_t sz(vcs.size());
 156   if(sz<1)
 157     throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
 158   if(sz==1)
 159     return vcs[0];
 160   MCAuto<MEDCouplingUMesh> p;
 161   {
 162     std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
 163     p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
 164   }
 165   p->zipCoords();
 166   {
 167     bool dummy; int dummy2;
 168     MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
 169   }
 170   MCAuto<DataArrayInt> edgeToKeep;
 171   MCAuto<MEDCouplingUMesh> p0;
 172   {
 173     MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
 174     p0=p->buildDescendingConnectivity(d,di,rd,rdi);
 175     MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
 176     edgeToKeep=dsi->findIdsEqual(1);
 177   }
 178   MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
 179   MCAuto<DataArrayDouble> eqn(skinOfRes->computePlaneEquationOf3DFaces());
 180   MCAuto<DataArrayInt> comm,commI;
 181   {
 182     DataArrayInt *a(0),*b(0);
 183     eqn->findCommonTuples(eps,0,a,b);
 184     comm=a; commI=b;
 185     //comm=DataArrayInt::New(); comm->alloc(0,1); commI=DataArrayInt::New(); commI->alloc(1,1); commI->setIJ(0,0,0);
 186   }
 187   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",3));
 188   ret->setCoords(skinOfRes->getCoords());
 189   ret->allocateCells();
 190   std::vector<int> conn;
 191   int jj(0);
 192   for(int i=0;i<commI->getNumberOfTuples()-1;i++,jj++)
 193     {
 194       if(jj!=0)
 195         conn.push_back(-1);
 196       MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(comm->begin()+commI->getIJ(i,0),comm->begin()+commI->getIJ(i+1,0),true));
 197       MCAuto<MEDCouplingUMesh> tmp2;
 198       if(commI->getIJ(i+1,0)-commI->getIJ(i,0)==1)
 199         tmp2=tmp;
 200       else
 201         tmp2=MergeVorCells2D(tmp,eps,false);
 202       const int *cPtr(tmp2->getNodalConnectivity()->begin()),*ciPtr(tmp2->getNodalConnectivityIndex()->begin());
 203       conn.insert(conn.end(),cPtr+1,cPtr+ciPtr[1]);
 204     }
 205   MCAuto<DataArrayInt> remain(comm->buildComplement(skinOfRes->getNumberOfCells()));
 206   {
 207     MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(remain->begin(),remain->end(),true));
 208     const int *cPtr(tmp->getNodalConnectivity()->begin()),*ciPtr(tmp->getNodalConnectivityIndex()->begin());
 209     for(int i=0;i<remain->getNumberOfTuples();i++,jj++)
 210       {
 211         if(jj!=0)
 212           conn.push_back(-1);
 213         conn.insert(conn.end(),cPtr+ciPtr[i]+1,cPtr+ciPtr[i+1]);
 214       }
 215   }
 216   ret->insertNextCell(INTERP_KERNEL::NORM_POLYHED,conn.size(),&conn[0]);
 217   return ret;
 218 }
 219
 220 MCAuto<MEDCouplingUMesh> MergeVorCells1D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
 221 {
 222   static const int CONN_SEG2_DFT[2]={0,1};
 223   if(vcs.empty())
 224     throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 1 !");
 225   if(vcs.size()==1)
 226     return vcs[0];
 227   if(vcs.size()>2)
 228     throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 2 !");
 229   double a0,b0,a1,b1;
 230   {
 231     const int *connPtr(vcs[0]->getNodalConnectivity()->begin());
 232     const double *coordPtr(vcs[0]->getCoords()->begin());
 233     a0=coordPtr[connPtr[1]]; b0=coordPtr[connPtr[2]];
 234   }
 235   {
 236     const int *connPtr(vcs[1]->getNodalConnectivity()->begin());
 237     const double *coordPtr(vcs[1]->getCoords()->begin());
 238     a1=coordPtr[connPtr[1]]; b1=coordPtr[connPtr[2]];
 239   }
 240   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",1)); ret->allocateCells(); ret->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN_SEG2_DFT);
 241   MCAuto<DataArrayDouble> coo(DataArrayDouble::New()); coo->alloc(2,1); ret->setCoords(coo);
 242   if(fabs(b0-a1)<eps)
 243     { coo->setIJ(0,0,a0); coo->setIJ(1,0,b1); }
 244   else if(fabs(b1-a0)<eps)
 245     { coo->setIJ(0,0,b0); coo->setIJ(1,0,a1); }
 246   return ret;
 247 }
 248
 249 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer1D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 250 {
 251   static const int CONN_SEG2_DFT[2]={0,1};
 252   if(!m || !points)
 253     throw INTERP_KERNEL::Exception("Voronoize1D : null pointer !");
 254   m->checkConsistencyLight();
 255   points->checkAllocated();
 256   if(m->getMeshDimension()!=1 || m->getSpaceDimension()!=1 || points->getNumberOfComponents()!=1)
 257     throw INTERP_KERNEL::Exception("Voronoize1D : spacedim must be equal to 1 and meshdim also equal to 1 !");
 258   if(m->getNumberOfCells()!=1)
 259     throw INTERP_KERNEL::Exception("Voronoize1D : mesh is expected to have only one cell !");
 260   int nbPts(points->getNumberOfTuples());
 261   if(nbPts<1)
 262     throw INTERP_KERNEL::Exception("Voronoize1D : at least one point expected !");
 263   std::vector<double> bbox(4);
 264   m->getBoundingBox(&bbox[0]);
 265   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 266   const double *pts(points->begin());
 267   for(int i=1;i<nbPts;i++)
 268     {
 269       MCAuto<MEDCouplingUMesh> vorTess;
 270       {
 271         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 272         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 273       }
 274       {
 275         bool dummy;
 276         int newNbNodes;
 277         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 278       }
 279       std::vector<int> polygsToIterOn;
 280       const double *pt(pts+i);
 281       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 282       if(polygsToIterOn.empty())
 283         throw INTERP_KERNEL::Exception("Voronoize1D : a point is outside domain !");
 284       if(polygsToIterOn.size()>2)
 285         throw INTERP_KERNEL::Exception("Voronoize1D : overlap of points !");
 286       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 287       for(std::vector<int>::const_iterator it=polygsToIterOn.begin();it!=polygsToIterOn.end();it++)
 288         {
 289           int poly(*it);
 290           //
 291           double seed(pts[poly]),zept(*pt);
 292           double mid((seed+zept)/2.);
 293           //
 294           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 295           tile->zipCoords();
 296           double a,b;
 297           {
 298             const int *connPtr(tile->getNodalConnectivity()->begin());
 299             const double *coordPtr(tile->getCoords()->begin());
 300             a=coordPtr[connPtr[1]]; b=coordPtr[connPtr[2]];
 301           }
 302           double pol0[2],pol1[2];
 303           MCAuto<DataArrayDouble> t0(DataArrayDouble::New()); t0->alloc(3,1); t0->setIJ(0,0,zept); t0->setIJ(1,0,mid); t0->setIJ(2,0,seed);
 304           t0->applyLin(1.,-a);
 305           if(t0->isMonotonic(true,eps))
 306             { pol0[0]=a; pol0[1]=mid; pol1[0]=mid; pol1[1]=b; }
 307           else
 308             { pol1[0]=a; pol1[1]=mid; pol0[0]=mid; pol0[1]=b; }
 309           MCAuto<MEDCouplingUMesh> modifiedCell(MEDCouplingUMesh::New("",1)); modifiedCell->allocateCells();
 310           MCAuto<DataArrayDouble> coo1(DataArrayDouble::New()); coo1->alloc(2,1); coo1->setIJ(0,0,pol1[0]); coo1->setIJ(1,0,pol1[1]);
 311           modifiedCell->setCoords(coo1); modifiedCell->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN_SEG2_DFT);
 312           //
 313           MCAuto<MEDCouplingUMesh> newVorCell(MEDCouplingUMesh::New("",1)); newVorCell->allocateCells();
 314           MCAuto<DataArrayDouble> coo2(DataArrayDouble::New()); coo2->alloc(2,1); coo2->setIJ(0,0,pol0[0]); coo2->setIJ(1,0,pol0[1]);
 315           newVorCell->setCoords(coo2); newVorCell->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN_SEG2_DFT);
 316           //
 317           l0[poly]=modifiedCell;
 318           newVorCells.push_back(newVorCell);
 319         }
 320       l0.push_back(MergeVorCells1D(newVorCells,eps));
 321     }
 322   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 323   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 324   {
 325     bool dummy; int dummy2;
 326     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 327   }
 328   return ret;
 329 }
 330
 331 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer2D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 332 {
 333   if(!m || !points)
 334     throw INTERP_KERNEL::Exception("Voronoize2D : null pointer !");
 335   m->checkConsistencyLight();
 336   points->checkAllocated();
 337   if(m->getMeshDimension()!=2 || m->getSpaceDimension()!=2 || points->getNumberOfComponents()!=2)
 338     throw INTERP_KERNEL::Exception("Voronoize2D : spacedim must be equal to 2 and meshdim also equal to 2 !");
 339   if(m->getNumberOfCells()!=1)
 340     throw INTERP_KERNEL::Exception("Voronoize2D : mesh is expected to have only one cell !");
 341   int nbPts(points->getNumberOfTuples());
 342   if(nbPts<1)
 343     throw INTERP_KERNEL::Exception("Voronoize2D : at least one point expected !");
 344   std::vector<double> bbox(4);
 345   m->getBoundingBox(&bbox[0]);
 346   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 347   const double *pts(points->begin());
 348   for(int i=1;i<nbPts;i++)
 349     {
 350       MCAuto<MEDCouplingUMesh> vorTess;
 351       {
 352         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 353         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 354       }
 355       {
 356         bool dummy;
 357         int newNbNodes;
 358         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 359       }
 360       std::vector<int> polygsToIterOn;
 361       const double *pt(pts+i*2);
 362       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 363       if(polygsToIterOn.size()<1)
 364         throw INTERP_KERNEL::Exception("Voronoize2D : presence of a point outside the given cell !");
 365       std::set<int> elemsToDo,elemsDone; elemsToDo.insert(polygsToIterOn[0]);
 366       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 367       while(!elemsToDo.empty())
 368         {
 369           int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
 370           const double *seed(pts+2*poly);
 371           MCAuto<MEDCouplingUMesh> cell(ComputeBigCellFrom(pt,seed,bbox,eps));
 372           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 373           tile->zipCoords();
 374           MCAuto<MEDCouplingUMesh> a;
 375           MCAuto<DataArrayInt> b,c;
 376           {
 377             DataArrayInt *bPtr(0),*cPtr(0);
 378             a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
 379             b=bPtr; c=cPtr;
 380           }
 381           MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
 382           if(part->getNumberOfTuples()!=1)
 383             throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
 384           MCAuto<MEDCouplingUMesh> newVorCell;
 385           {
 386             MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
 387             newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
 388           }
 389           newVorCell->zipCoords();
 390           MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
 391           modifiedCell->zipCoords();
 392           l0[poly]=modifiedCell;
 393           //
 394           MCAuto<DataArrayInt> ids;
 395           {
 396             DataArrayInt *tmp(0);
 397             bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
 398             ids=tmp;
 399             if(!sta)
 400               throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
 401           }
 402           MCAuto<DataArrayDouble> newCoords;
 403           {
 404             MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
 405             newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
 406           }
 407           const double *cPtr(newCoords->begin());
 408           for(int i=0;i<newCoords->getNumberOfTuples();i++,cPtr+=2)
 409             {
 410               std::set<int> zeCandidates;
 411               {
 412                 std::vector<int> zeCandidatesTmp;
 413                 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
 414                 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
 415               }
 416               std::set<int> tmp,newElementsToDo;
 417               std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
 418               std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
 419               elemsToDo=newElementsToDo;
 420             }
 421           newVorCells.push_back(newVorCell);
 422         }
 423       l0.push_back(MergeVorCells(newVorCells,eps));
 424     }
 425   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 426   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 427   {
 428     bool dummy; int dummy2;
 429     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 430   }
 431   return ret;
 432 }
 433
 434 MCAuto<MEDCouplingUMesh> Split3DCellInParts(const MEDCouplingUMesh *m, const double pt[3], const double seed[3], double eps, int tmp[2])
 435 {
 436   if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || m->getNumberOfCells()!=1)
 437     throw INTERP_KERNEL::Exception("Split3DCellInParts : expecting a 3D with exactly one cell !");
 438   double middle[3]={(pt[0]+seed[0])/2.,(pt[1]+seed[1])/2.,(pt[2]+seed[2])/2.};
 439   double vec[3]={pt[0]-seed[0],pt[1]-seed[1],pt[2]-seed[2]};
 440   MCAuto<MEDCouplingUMesh> res(m->clipSingle3DCellByPlane(middle,vec,eps));
 441   return res;
 442 }
 443
 444 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 445 {
 446   double eps2(1.-sqrt(eps));// 2nd eps for interpolation. Here the eps is computed to feet cos(eps) ~ 1-eps^2
 447   if(!m || !points)
 448     throw INTERP_KERNEL::Exception("Voronoize3D : null pointer !");
 449   m->checkConsistencyLight();
 450   points->checkAllocated();
 451   if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || points->getNumberOfComponents()!=3)
 452     throw INTERP_KERNEL::Exception("Voronoize3D : spacedim must be equal to 3 and meshdim also equal to 3 !");
 453   if(m->getNumberOfCells()!=1)
 454     throw INTERP_KERNEL::Exception("Voronoize3D : mesh is expected to have only one cell !");
 455   int nbPts(points->getNumberOfTuples());
 456   if(nbPts<1)
 457     throw INTERP_KERNEL::Exception("Voronoize3D : at least one point expected !");
 458   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 459   const double *pts(points->begin());
 460   for(int i=1;i<nbPts;i++)
 461     {
 462       MCAuto<MEDCouplingUMesh> vorTess;
 463       {
 464         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 465         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 466       }
 467       {
 468         bool dummy;
 469         int newNbNodes;
 470         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 471       }
 472       std::vector<int> polygsToIterOn;
 473       const double *pt(pts+i*3);
 474       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 475       if(polygsToIterOn.size()<1)
 476         throw INTERP_KERNEL::Exception("Voronoize3D : presence of a point outside the given cell !");
 477       std::set<int> elemsToDo(polygsToIterOn.begin(),polygsToIterOn.end()),elemsDone;
 478       std::size_t ii(0);
 479       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 480       MCAuto<DataArrayInt> d(DataArrayInt::New()),dI(DataArrayInt::New()),rd(DataArrayInt::New()),rdI(DataArrayInt::New());
 481       MCAuto<MEDCouplingUMesh> faces(vorTess->buildDescendingConnectivity(d,dI,rd,rdI));
 482       //
 483       while(!elemsToDo.empty())
 484         {
 485           int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
 486           const double *seed(pts+3*poly);
 487           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 488           tile->zipCoords();
 489           int tmp[2];
 490           MCAuto<MEDCouplingUMesh> cells;
 491           try
 492             {
 493               cells=Split3DCellInParts(tile,pt,seed,eps,tmp);
 494             }
 495           catch(INTERP_KERNEL::Exception& e)
 496             {
 497               continue;
 498             }
 499           MCAuto<MEDCouplingUMesh> newVorCell(cells->buildPartOfMySelfSlice(1,2,1,true));
 500           newVorCell->zipCoords();
 501           MCAuto<MEDCouplingUMesh> modifiedCell(cells->buildPartOfMySelfSlice(0,1,1,true));
 502           modifiedCell->zipCoords();
 503           l0[poly]=modifiedCell;
 504           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
 505             {
 506               MCAuto<MEDCouplingUMesh> faces2;
 507               {
 508                 MCAuto<DataArrayInt> d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New());
 509                 faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I);
 510               }
 511               MCAuto<MEDCouplingUMesh> faces3(faces2->buildPartOfMySelfSlice(1,faces2->getNumberOfCells(),1,true));// suppress internal face
 512               MCAuto<MEDCouplingUMesh> facesOfCurSplitPol(faces->buildPartOfMySelf(d->begin()+dI->getIJ(poly,0),d->begin()+dI->getIJ(poly+1,0),true));
 513               // intersection between the out faces of newVorCell and the neighbor faces of poly polyhedron -> candidates
 514               MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(facesOfCurSplitPol);
 515               MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(faces3);
 516               INTERP_KERNEL::Interpolation3DSurf interpolation;
 517               interpolation.setMinDotBtwPlane3DSurfIntersect(eps2);
 518               interpolation.setMaxDistance3DSurfIntersect(eps);
 519               interpolation.setPrecision(1e-12);
 520               std::vector<std::map<int,double> > matrix;
 521               interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix,"P0P0");
 522               std::set<int> zeCandidates;
 523               for(std::vector<std::map<int,double> >::const_iterator it2=matrix.begin();it2!=matrix.end();it2++)
 524                 for(std::map<int,double>::const_iterator it3=(*it2).begin();it3!=(*it2).end();it3++)
 525                   {
 526                     int faceIdInVorTess(d->getIJ(dI->getIJ(poly,0)+(*it3).first,0));
 527                     for(const int *it4=rd->begin()+rdI->getIJ(faceIdInVorTess,0);it4!=rd->begin()+rdI->getIJ(faceIdInVorTess+1,0);it4++)
 528                       {
 529                         if(*it4!=poly)
 530                           zeCandidates.insert(*it4);
 531                       }
 532                   }
 533               std::set<int> tmp,newElementsToDo;
 534               std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
 535               std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
 536               elemsToDo=newElementsToDo;
 537             }
 538           //
 539           newVorCells.push_back(newVorCell);
 540           ii++;
 541         }
 542       l0.push_back(MergeVorCells3D(newVorCells,eps));
 543     }
 544   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 545   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 546   {
 547     bool dummy; int dummy2;
 548     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 549   }
 550   return ret;
 551 }