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 /*!
 154  * suppress additional sub points on edges
 155  */
 156 MCAuto<MEDCouplingUMesh> SimplifyPolygon(const MEDCouplingUMesh *m, double eps)
 157 {
 158   if(m->getNumberOfCells()!=1)
 159     throw INTERP_KERNEL::Exception("SimplifyPolygon : internal error !");
 160   const int *conn(m->getNodalConnectivity()->begin()),*conni(m->getNodalConnectivityIndex()->begin());
 161   int nbPtsInPolygon(conni[1]-conni[0]-1);
 162   const double *coo(m->getCoords()->begin());
 163   std::vector<int> resConn;
 164   for(int i=0;i<nbPtsInPolygon;i++)
 165     {
 166       int prev(conn[(i+nbPtsInPolygon-1)%nbPtsInPolygon+1]),current(conn[i%nbPtsInPolygon+1]),zeNext(conn[(i+1)%nbPtsInPolygon+1]);
 167       double a[3]={
 168         coo[3*prev+0]-coo[3*current+0],
 169         coo[3*prev+1]-coo[3*current+1],
 170         coo[3*prev+2]-coo[3*current+2],
 171       },b[3]={
 172         coo[3*current+0]-coo[3*zeNext+0],
 173         coo[3*current+1]-coo[3*zeNext+1],
 174         coo[3*current+2]-coo[3*zeNext+2],
 175       };
 176       double c[3]={a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]};
 177       if(sqrt(c[0]*c[0]+c[1]*c[1]+c[2]*c[2])>eps)
 178         resConn.push_back(current);
 179     }
 180   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
 181   ret->setCoords(m->getCoords());
 182   ret->allocateCells();
 183   ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,resConn.size(),&resConn[0]);
 184   return ret;
 185 }
 186
 187 MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
 188 {
 189   std::size_t sz(vcs.size());
 190   if(sz<1)
 191     throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
 192   if(sz==1)
 193     return vcs[0];
 194   MCAuto<MEDCouplingUMesh> p;
 195   {
 196     std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
 197     p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
 198   }
 199   p->zipCoords();
 200   {
 201     bool dummy; int dummy2;
 202     MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
 203   }
 204   MCAuto<DataArrayInt> edgeToKeep;
 205   MCAuto<MEDCouplingUMesh> p0;
 206   {
 207     MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
 208     p0=p->buildDescendingConnectivity(d,di,rd,rdi);
 209     MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
 210     edgeToKeep=dsi->findIdsEqual(1);
 211   }
 212   MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
 213   MCAuto<DataArrayDouble> eqn(skinOfRes->computePlaneEquationOf3DFaces());
 214   MCAuto<DataArrayInt> comm,commI;
 215   {
 216     DataArrayInt *a(0),*b(0);
 217     eqn->findCommonTuples(eps,0,a,b);
 218     comm=a; commI=b;
 219     //comm=DataArrayInt::New(); comm->alloc(0,1); commI=DataArrayInt::New(); commI->alloc(1,1); commI->setIJ(0,0,0);
 220   }
 221   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",3));
 222   ret->setCoords(skinOfRes->getCoords());
 223   ret->allocateCells();
 224   std::vector<int> conn;
 225   int jj(0);
 226   for(int i=0;i<commI->getNumberOfTuples()-1;i++,jj++)
 227     {
 228       if(jj!=0)
 229         conn.push_back(-1);
 230       MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(comm->begin()+commI->getIJ(i,0),comm->begin()+commI->getIJ(i+1,0),true));
 231       MCAuto<MEDCouplingUMesh> tmp2;
 232       if(commI->getIJ(i+1,0)-commI->getIJ(i,0)==1)
 233         tmp2=tmp;
 234       else
 235         tmp2=MergeVorCells2D(tmp,eps,false);
 236       tmp2=SimplifyPolygon(tmp2,eps);
 237       const int *cPtr(tmp2->getNodalConnectivity()->begin()),*ciPtr(tmp2->getNodalConnectivityIndex()->begin());
 238       conn.insert(conn.end(),cPtr+1,cPtr+ciPtr[1]);
 239     }
 240   MCAuto<DataArrayInt> remain(comm->buildComplement(skinOfRes->getNumberOfCells()));
 241   {
 242     MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(remain->begin(),remain->end(),true));
 243     const int *cPtr(tmp->getNodalConnectivity()->begin()),*ciPtr(tmp->getNodalConnectivityIndex()->begin());
 244     for(int i=0;i<remain->getNumberOfTuples();i++,jj++)
 245       {
 246         if(jj!=0)
 247           conn.push_back(-1);
 248         conn.insert(conn.end(),cPtr+ciPtr[i]+1,cPtr+ciPtr[i+1]);
 249       }
 250   }
 251   ret->insertNextCell(INTERP_KERNEL::NORM_POLYHED,conn.size(),&conn[0]);
 252   return ret;
 253 }
 254
 255 MCAuto<MEDCouplingUMesh> MergeVorCells1D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
 256 {
 257   static const int CONN_SEG2_DFT[2]={0,1};
 258   if(vcs.empty())
 259     throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 1 !");
 260   if(vcs.size()==1)
 261     return vcs[0];
 262   if(vcs.size()>2)
 263     throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 2 !");
 264   double a0,b0,a1,b1;
 265   {
 266     const int *connPtr(vcs[0]->getNodalConnectivity()->begin());
 267     const double *coordPtr(vcs[0]->getCoords()->begin());
 268     a0=coordPtr[connPtr[1]]; b0=coordPtr[connPtr[2]];
 269   }
 270   {
 271     const int *connPtr(vcs[1]->getNodalConnectivity()->begin());
 272     const double *coordPtr(vcs[1]->getCoords()->begin());
 273     a1=coordPtr[connPtr[1]]; b1=coordPtr[connPtr[2]];
 274   }
 275   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",1)); ret->allocateCells(); ret->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN_SEG2_DFT);
 276   MCAuto<DataArrayDouble> coo(DataArrayDouble::New()); coo->alloc(2,1); ret->setCoords(coo);
 277   if(fabs(b0-a1)<eps)
 278     { coo->setIJ(0,0,a0); coo->setIJ(1,0,b1); }
 279   else if(fabs(b1-a0)<eps)
 280     { coo->setIJ(0,0,b0); coo->setIJ(1,0,a1); }
 281   return ret;
 282 }
 283
 284 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer1D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 285 {
 286   static const int CONN_SEG2_DFT[2]={0,1};
 287   if(!m || !points)
 288     throw INTERP_KERNEL::Exception("Voronoize1D : null pointer !");
 289   m->checkConsistencyLight();
 290   points->checkAllocated();
 291   if(m->getMeshDimension()!=1 || m->getSpaceDimension()!=1 || points->getNumberOfComponents()!=1)
 292     throw INTERP_KERNEL::Exception("Voronoize1D : spacedim must be equal to 1 and meshdim also equal to 1 !");
 293   if(m->getNumberOfCells()!=1)
 294     throw INTERP_KERNEL::Exception("Voronoize1D : mesh is expected to have only one cell !");
 295   int nbPts(points->getNumberOfTuples());
 296   if(nbPts<1)
 297     throw INTERP_KERNEL::Exception("Voronoize1D : at least one point expected !");
 298   std::vector<double> bbox(4);
 299   m->getBoundingBox(&bbox[0]);
 300   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 301   const double *pts(points->begin());
 302   for(int i=1;i<nbPts;i++)
 303     {
 304       MCAuto<MEDCouplingUMesh> vorTess;
 305       {
 306         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 307         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 308       }
 309       {
 310         bool dummy;
 311         int newNbNodes;
 312         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 313       }
 314       std::vector<int> polygsToIterOn;
 315       const double *pt(pts+i);
 316       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 317       if(polygsToIterOn.empty())
 318         throw INTERP_KERNEL::Exception("Voronoize1D : a point is outside domain !");
 319       if(polygsToIterOn.size()>2)
 320         throw INTERP_KERNEL::Exception("Voronoize1D : overlap of points !");
 321       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 322       for(std::vector<int>::const_iterator it=polygsToIterOn.begin();it!=polygsToIterOn.end();it++)
 323         {
 324           int poly(*it);
 325           //
 326           double seed(pts[poly]),zept(*pt);
 327           double mid((seed+zept)/2.);
 328           //
 329           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 330           tile->zipCoords();
 331           double a,b;
 332           {
 333             const int *connPtr(tile->getNodalConnectivity()->begin());
 334             const double *coordPtr(tile->getCoords()->begin());
 335             a=coordPtr[connPtr[1]]; b=coordPtr[connPtr[2]];
 336           }
 337           double pol0[2],pol1[2];
 338           MCAuto<DataArrayDouble> t0(DataArrayDouble::New()); t0->alloc(3,1); t0->setIJ(0,0,zept); t0->setIJ(1,0,mid); t0->setIJ(2,0,seed);
 339           t0->applyLin(1.,-a);
 340           if(t0->isMonotonic(true,eps))
 341             { pol0[0]=a; pol0[1]=mid; pol1[0]=mid; pol1[1]=b; }
 342           else
 343             { pol1[0]=a; pol1[1]=mid; pol0[0]=mid; pol0[1]=b; }
 344           MCAuto<MEDCouplingUMesh> modifiedCell(MEDCouplingUMesh::New("",1)); modifiedCell->allocateCells();
 345           MCAuto<DataArrayDouble> coo1(DataArrayDouble::New()); coo1->alloc(2,1); coo1->setIJ(0,0,pol1[0]); coo1->setIJ(1,0,pol1[1]);
 346           modifiedCell->setCoords(coo1); modifiedCell->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN_SEG2_DFT);
 347           //
 348           MCAuto<MEDCouplingUMesh> newVorCell(MEDCouplingUMesh::New("",1)); newVorCell->allocateCells();
 349           MCAuto<DataArrayDouble> coo2(DataArrayDouble::New()); coo2->alloc(2,1); coo2->setIJ(0,0,pol0[0]); coo2->setIJ(1,0,pol0[1]);
 350           newVorCell->setCoords(coo2); newVorCell->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN_SEG2_DFT);
 351           //
 352           l0[poly]=modifiedCell;
 353           newVorCells.push_back(newVorCell);
 354         }
 355       l0.push_back(MergeVorCells1D(newVorCells,eps));
 356     }
 357   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 358   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 359   {
 360     bool dummy; int dummy2;
 361     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 362   }
 363   return ret;
 364 }
 365
 366 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer2D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 367 {
 368   if(!m || !points)
 369     throw INTERP_KERNEL::Exception("Voronoize2D : null pointer !");
 370   m->checkConsistencyLight();
 371   points->checkAllocated();
 372   if(m->getMeshDimension()!=2 || m->getSpaceDimension()!=2 || points->getNumberOfComponents()!=2)
 373     throw INTERP_KERNEL::Exception("Voronoize2D : spacedim must be equal to 2 and meshdim also equal to 2 !");
 374   if(m->getNumberOfCells()!=1)
 375     throw INTERP_KERNEL::Exception("Voronoize2D : mesh is expected to have only one cell !");
 376   int nbPts(points->getNumberOfTuples());
 377   if(nbPts<1)
 378     throw INTERP_KERNEL::Exception("Voronoize2D : at least one point expected !");
 379   std::vector<double> bbox(4);
 380   m->getBoundingBox(&bbox[0]);
 381   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 382   const double *pts(points->begin());
 383   for(int i=1;i<nbPts;i++)
 384     {
 385       MCAuto<MEDCouplingUMesh> vorTess;
 386       {
 387         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 388         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 389       }
 390       {
 391         bool dummy;
 392         int newNbNodes;
 393         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 394       }
 395       std::vector<int> polygsToIterOn;
 396       const double *pt(pts+i*2);
 397       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 398       if(polygsToIterOn.size()<1)
 399         throw INTERP_KERNEL::Exception("Voronoize2D : presence of a point outside the given cell !");
 400       std::set<int> elemsToDo,elemsDone; elemsToDo.insert(polygsToIterOn[0]);
 401       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 402       while(!elemsToDo.empty())
 403         {
 404           int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
 405           const double *seed(pts+2*poly);
 406           MCAuto<MEDCouplingUMesh> cell(ComputeBigCellFrom(pt,seed,bbox,eps));
 407           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 408           tile->zipCoords();
 409           MCAuto<MEDCouplingUMesh> a;
 410           MCAuto<DataArrayInt> b,c;
 411           {
 412             DataArrayInt *bPtr(0),*cPtr(0);
 413             a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
 414             b=bPtr; c=cPtr;
 415           }
 416           MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
 417           if(part->getNumberOfTuples()!=1)
 418             throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
 419           MCAuto<MEDCouplingUMesh> newVorCell;
 420           {
 421             MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
 422             newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
 423           }
 424           newVorCell->zipCoords();
 425           MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
 426           modifiedCell->zipCoords();
 427           l0[poly]=modifiedCell;
 428           //
 429           MCAuto<DataArrayInt> ids;
 430           {
 431             DataArrayInt *tmp(0);
 432             bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
 433             ids=tmp;
 434             if(!sta)
 435               throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
 436           }
 437           MCAuto<DataArrayDouble> newCoords;
 438           {
 439             MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
 440             newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
 441           }
 442           const double *cPtr(newCoords->begin());
 443           for(int j=0;j<newCoords->getNumberOfTuples();j++,cPtr+=2)
 444             {
 445               std::set<int> zeCandidates;
 446               {
 447                 std::vector<int> zeCandidatesTmp;
 448                 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
 449                 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
 450               }
 451               std::set<int> tmp2,newElementsToDo;
 452               std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp2,tmp2.begin()));
 453               std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp2.begin(),tmp2.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
 454               elemsToDo=newElementsToDo;
 455             }
 456           newVorCells.push_back(newVorCell);
 457         }
 458       l0.push_back(MergeVorCells(newVorCells,eps));
 459     }
 460   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 461   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 462   {
 463     bool dummy; int dummy2;
 464     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 465   }
 466   return ret;
 467 }
 468
 469 MCAuto<MEDCouplingUMesh> Split3DCellInParts(const MEDCouplingUMesh *m, const double pt[3], const double seed[3], double eps, int tmp[2])
 470 {
 471   if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || m->getNumberOfCells()!=1)
 472     throw INTERP_KERNEL::Exception("Split3DCellInParts : expecting a 3D with exactly one cell !");
 473   double middle[3]={(pt[0]+seed[0])/2.,(pt[1]+seed[1])/2.,(pt[2]+seed[2])/2.};
 474   double vec[3]={pt[0]-seed[0],pt[1]-seed[1],pt[2]-seed[2]};
 475   MCAuto<MEDCouplingUMesh> res(m->clipSingle3DCellByPlane(middle,vec,eps));
 476   return res;
 477 }
 478
 479 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
 480 {
 481   double eps2(1.-sqrt(eps));// 2nd eps for interpolation. Here the eps is computed to feet cos(eps) ~ 1-eps^2
 482   if(!m || !points)
 483     throw INTERP_KERNEL::Exception("Voronoize3D : null pointer !");
 484   m->checkConsistencyLight();
 485   points->checkAllocated();
 486   if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || points->getNumberOfComponents()!=3)
 487     throw INTERP_KERNEL::Exception("Voronoize3D : spacedim must be equal to 3 and meshdim also equal to 3 !");
 488   if(m->getNumberOfCells()!=1)
 489     throw INTERP_KERNEL::Exception("Voronoize3D : mesh is expected to have only one cell !");
 490   int nbPts(points->getNumberOfTuples());
 491   if(nbPts<1)
 492     throw INTERP_KERNEL::Exception("Voronoize3D : at least one point expected !");
 493   std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
 494   const double *pts(points->begin());
 495   for(int i=1;i<nbPts;i++)
 496     {
 497       MCAuto<MEDCouplingUMesh> vorTess;
 498       {
 499         std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 500         vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
 501       }
 502       {
 503         bool dummy;
 504         int newNbNodes;
 505         MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
 506       }
 507       std::vector<int> polygsToIterOn;
 508       const double *pt(pts+i*3);
 509       vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
 510       if(polygsToIterOn.size()<1)
 511         throw INTERP_KERNEL::Exception("Voronoize3D : presence of a point outside the given cell !");
 512       std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
 513       for(int poly=0;poly<vorTess->getNumberOfCells();poly++)
 514         {
 515           const double *seed(pts+3*poly);
 516           MCAuto<MEDCouplingUMesh> tile(l0[poly]);
 517           tile->zipCoords();
 518           int tmp[2];
 519           MCAuto<MEDCouplingUMesh> cells;
 520           try
 521             {
 522               cells=Split3DCellInParts(tile,pt,seed,eps,tmp);
 523             }
 524           catch(INTERP_KERNEL::Exception& e)
 525             {
 526               continue;
 527             }
 528           MCAuto<MEDCouplingUMesh> newVorCell(cells->buildPartOfMySelfSlice(1,2,1,true));
 529           newVorCell->zipCoords();
 530           MCAuto<MEDCouplingUMesh> modifiedCell(cells->buildPartOfMySelfSlice(0,1,1,true));
 531           modifiedCell->zipCoords();
 532           newVorCells.push_back(newVorCell);
 533           l0[poly]=modifiedCell;
 534         }
 535       l0.push_back(MergeVorCells3D(newVorCells,eps));
 536     }
 537   std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
 538   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
 539   {
 540     bool dummy; int dummy2;
 541     MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
 542   }
 543   return ret;
 544 }