1 // Copyright (C) 2007-2017 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (EDF R&D)
21 #include "MEDCouplingVoronoi.hxx"
22 #include "MEDCoupling1GTUMesh.hxx"
23 #include "MEDCouplingCMesh.hxx"
26 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
27 #include "Interpolation2D.txx"
28 #include "Interpolation3DSurf.hxx"
30 using namespace MEDCoupling;
32 Voronizer::~Voronizer()
36 int Voronizer1D::getDimension() const
41 int Voronizer2D::getDimension() const
46 int Voronizer3D::getDimension() const
51 MCAuto<MEDCouplingUMesh> ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector<double>& bbox, double eps)
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.};
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];
67 double nor(sqrt(v[0]*v[0]+v[1]*v[1]));
70 MCAuto<MEDCouplingUMesh> line(MEDCouplingUMesh::New("line",1));
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]);
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;
81 DataArrayInt *cellNb1(0),*cellNb2(0);
82 MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0);
83 MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2);
85 MCAuto<DataArrayInt> cellNb10(cellNb1),cellNb20(cellNb2);
88 sp2->getCellsContainingPoint(pt1,eps,ccp);
90 throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !");
91 MCAuto<MEDCouplingUMesh> ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true));
97 MCAuto<MEDCouplingUMesh> MergeVorCells2D(MEDCouplingUMesh *p, double eps, bool isZip)
99 MCAuto<DataArrayInt> edgeToKeep;
100 MCAuto<MEDCouplingUMesh> p0;
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);
107 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
110 skinOfRes->zipCoords();
111 if(skinOfRes->getNumberOfCells()!=skinOfRes->getNumberOfNodes())
112 throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !");
114 MCAuto<DataArrayInt> d(skinOfRes->orderConsecutiveCells1D());
115 MCAuto<MEDCoupling1SGTUMesh> skinOfRes2;
117 MCAuto<MEDCouplingUMesh> part(skinOfRes->buildPartOfMySelf(d->begin(),d->end()));
118 skinOfRes2=MEDCoupling1SGTUMesh::New(part);
120 MCAuto<DataArrayInt> c(skinOfRes2->getNodalConnectivity()->deepCopy());
121 c->circularPermutation(1);
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());
129 m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin());
133 MCAuto<MEDCouplingUMesh> MergeVorCells(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
135 std::size_t sz(vcs.size());
137 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
140 MCAuto<MEDCouplingUMesh> p;
142 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
143 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
147 bool dummy; int dummy2;
148 MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
150 return MergeVorCells2D(p,eps,true);
154 * suppress additional sub points on edges
156 MCAuto<MEDCouplingUMesh> SimplifyPolygon(const MEDCouplingUMesh *m, double eps)
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++)
166 int prev(conn[(i+nbPtsInPolygon-1)%nbPtsInPolygon+1]),current(conn[i%nbPtsInPolygon+1]),zeNext(conn[(i+1)%nbPtsInPolygon+1]);
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],
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],
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);
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]);
187 MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
189 std::size_t sz(vcs.size());
191 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
194 MCAuto<MEDCouplingUMesh> p;
196 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
197 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
201 bool dummy; int dummy2;
202 MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
204 MCAuto<DataArrayInt> edgeToKeep;
205 MCAuto<MEDCouplingUMesh> p0;
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);
212 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
213 MCAuto<DataArrayDouble> eqn(skinOfRes->computePlaneEquationOf3DFaces());
214 MCAuto<DataArrayInt> comm,commI;
216 DataArrayInt *a(0),*b(0);
217 eqn->findCommonTuples(eps,0,a,b);
219 //comm=DataArrayInt::New(); comm->alloc(0,1); commI=DataArrayInt::New(); commI->alloc(1,1); commI->setIJ(0,0,0);
221 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",3));
222 ret->setCoords(skinOfRes->getCoords());
223 ret->allocateCells();
224 std::vector<int> conn;
226 for(int i=0;i<commI->getNumberOfTuples()-1;i++,jj++)
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)
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]);
240 MCAuto<DataArrayInt> remain(comm->buildComplement(skinOfRes->getNumberOfCells()));
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++)
248 conn.insert(conn.end(),cPtr+ciPtr[i]+1,cPtr+ciPtr[i+1]);
251 ret->insertNextCell(INTERP_KERNEL::NORM_POLYHED,conn.size(),&conn[0]);
255 MCAuto<MEDCouplingUMesh> MergeVorCells1D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
257 static const int CONN_SEG2_DFT[2]={0,1};
259 throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 1 !");
263 throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 2 !");
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]];
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]];
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);
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); }
284 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer1D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
286 static const int CONN_SEG2_DFT[2]={0,1};
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());
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++)
304 MCAuto<MEDCouplingUMesh> vorTess;
306 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
307 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
312 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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++)
326 double seed(pts[poly]),zept(*pt);
327 double mid((seed+zept)/2.);
329 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
333 const int *connPtr(tile->getNodalConnectivity()->begin());
334 const double *coordPtr(tile->getCoords()->begin());
335 a=coordPtr[connPtr[1]]; b=coordPtr[connPtr[2]];
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);
340 if(t0->isMonotonic(true,eps))
341 { pol0[0]=a; pol0[1]=mid; pol1[0]=mid; pol1[1]=b; }
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);
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);
352 l0[poly]=modifiedCell;
353 newVorCells.push_back(newVorCell);
355 l0.push_back(MergeVorCells1D(newVorCells,eps));
357 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
358 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
360 bool dummy; int dummy2;
361 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
366 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer2D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
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());
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++)
385 MCAuto<MEDCouplingUMesh> vorTess;
387 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
388 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
393 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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())
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]);
409 MCAuto<MEDCouplingUMesh> a;
410 MCAuto<DataArrayInt> b,c;
412 DataArrayInt *bPtr(0),*cPtr(0);
413 a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
416 MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
417 if(part->getNumberOfTuples()!=1)
418 throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
419 MCAuto<MEDCouplingUMesh> newVorCell;
421 MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
422 newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
424 newVorCell->zipCoords();
425 MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
426 modifiedCell->zipCoords();
427 l0[poly]=modifiedCell;
429 MCAuto<DataArrayInt> ids;
431 DataArrayInt *tmp(0);
432 bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
435 throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
437 MCAuto<DataArrayDouble> newCoords;
439 MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
440 newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
442 const double *cPtr(newCoords->begin());
443 for(int j=0;j<newCoords->getNumberOfTuples();j++,cPtr+=2)
445 std::set<int> zeCandidates;
447 std::vector<int> zeCandidatesTmp;
448 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
449 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
451 std::set<int> tmp,newElementsToDo;
452 std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
453 std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
454 elemsToDo=newElementsToDo;
456 newVorCells.push_back(newVorCell);
458 l0.push_back(MergeVorCells(newVorCells,eps));
460 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
461 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
463 bool dummy; int dummy2;
464 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
469 MCAuto<MEDCouplingUMesh> Split3DCellInParts(const MEDCouplingUMesh *m, const double pt[3], const double seed[3], double eps, int tmp[2])
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));
479 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
481 double eps2(1.-sqrt(eps));// 2nd eps for interpolation. Here the eps is computed to feet cos(eps) ~ 1-eps^2
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());
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++)
497 MCAuto<MEDCouplingUMesh> vorTess;
499 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
500 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
505 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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::set<int> elemsToDo(polygsToIterOn.begin(),polygsToIterOn.end()),elemsDone;
514 std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
515 MCAuto<DataArrayInt> d(DataArrayInt::New()),dI(DataArrayInt::New()),rd(DataArrayInt::New()),rdI(DataArrayInt::New());
516 MCAuto<MEDCouplingUMesh> faces(vorTess->buildDescendingConnectivity(d,dI,rd,rdI));
518 while(!elemsToDo.empty())
520 int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
521 const double *seed(pts+3*poly);
522 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
525 MCAuto<MEDCouplingUMesh> cells;
528 cells=Split3DCellInParts(tile,pt,seed,eps,tmp);
530 catch(INTERP_KERNEL::Exception& e)
534 MCAuto<MEDCouplingUMesh> newVorCell(cells->buildPartOfMySelfSlice(1,2,1,true));
535 newVorCell->zipCoords();
536 MCAuto<MEDCouplingUMesh> modifiedCell(cells->buildPartOfMySelfSlice(0,1,1,true));
537 modifiedCell->zipCoords();
538 l0[poly]=modifiedCell;
539 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
541 MCAuto<MEDCouplingUMesh> faces2;
543 MCAuto<DataArrayInt> d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New());
544 faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I);
546 MCAuto<MEDCouplingUMesh> faces3(faces2->buildPartOfMySelfSlice(1,faces2->getNumberOfCells(),1,true));// suppress internal face
547 MCAuto<MEDCouplingUMesh> facesOfCurSplitPol(faces->buildPartOfMySelf(d->begin()+dI->getIJ(poly,0),d->begin()+dI->getIJ(poly+1,0),true));
548 // intersection between the out faces of newVorCell and the neighbor faces of poly polyhedron -> candidates
549 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(facesOfCurSplitPol);
550 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(faces3);
551 INTERP_KERNEL::Interpolation3DSurf interpolation;
552 interpolation.setMinDotBtwPlane3DSurfIntersect(eps2);
553 interpolation.setMaxDistance3DSurfIntersect(eps);
554 interpolation.setPrecision(1e-12);
555 std::vector<std::map<int,double> > matrix;
556 interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix,"P0P0");
557 std::set<int> zeCandidates;
558 for(std::vector<std::map<int,double> >::const_iterator it2=matrix.begin();it2!=matrix.end();it2++)
559 for(std::map<int,double>::const_iterator it3=(*it2).begin();it3!=(*it2).end();it3++)
561 int faceIdInVorTess(d->getIJ(dI->getIJ(poly,0)+(*it3).first,0));
562 for(const int *it4=rd->begin()+rdI->getIJ(faceIdInVorTess,0);it4!=rd->begin()+rdI->getIJ(faceIdInVorTess+1,0);it4++)
565 zeCandidates.insert(*it4);
568 std::set<int> tmp2,newElementsToDo;
569 std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp2,tmp2.begin()));
570 std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp2.begin(),tmp2.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
571 elemsToDo=newElementsToDo;
574 newVorCells.push_back(newVorCell);
575 l0[poly]=modifiedCell;
577 l0.push_back(MergeVorCells3D(newVorCells,eps));
579 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
580 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
582 bool dummy; int dummy2;
583 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));