src/MEDCoupling/MEDCouplingVoronoi.cxx

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