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);
153 MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
155 std::size_t sz(vcs.size());
157 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
160 MCAuto<MEDCouplingUMesh> p;
162 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
163 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
167 bool dummy; int dummy2;
168 MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
170 MCAuto<DataArrayInt> edgeToKeep;
171 MCAuto<MEDCouplingUMesh> p0;
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);
178 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
179 MCAuto<DataArrayDouble> eqn(skinOfRes->computePlaneEquationOf3DFaces());
180 MCAuto<DataArrayInt> comm,commI;
182 DataArrayInt *a(0),*b(0);
183 eqn->findCommonTuples(eps,0,a,b);
185 //comm=DataArrayInt::New(); comm->alloc(0,1); commI=DataArrayInt::New(); commI->alloc(1,1); commI->setIJ(0,0,0);
187 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",3));
188 ret->setCoords(skinOfRes->getCoords());
189 ret->allocateCells();
190 std::vector<int> conn;
192 for(int i=0;i<commI->getNumberOfTuples()-1;i++,jj++)
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)
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]);
205 MCAuto<DataArrayInt> remain(comm->buildComplement(skinOfRes->getNumberOfCells()));
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++)
213 conn.insert(conn.end(),cPtr+ciPtr[i]+1,cPtr+ciPtr[i+1]);
216 ret->insertNextCell(INTERP_KERNEL::NORM_POLYHED,conn.size(),&conn[0]);
220 MCAuto<MEDCouplingUMesh> MergeVorCells1D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
222 static const int CONN_SEG2_DFT[2]={0,1};
224 throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 1 !");
228 throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 2 !");
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]];
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]];
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);
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); }
249 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer1D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
251 static const int CONN_SEG2_DFT[2]={0,1};
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());
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++)
269 MCAuto<MEDCouplingUMesh> vorTess;
271 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
272 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
277 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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++)
291 double seed(pts[poly]),zept(*pt);
292 double mid((seed+zept)/2.);
294 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
298 const int *connPtr(tile->getNodalConnectivity()->begin());
299 const double *coordPtr(tile->getCoords()->begin());
300 a=coordPtr[connPtr[1]]; b=coordPtr[connPtr[2]];
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);
305 if(t0->isMonotonic(true,eps))
306 { pol0[0]=a; pol0[1]=mid; pol1[0]=mid; pol1[1]=b; }
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);
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);
317 l0[poly]=modifiedCell;
318 newVorCells.push_back(newVorCell);
320 l0.push_back(MergeVorCells1D(newVorCells,eps));
322 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
323 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
325 bool dummy; int dummy2;
326 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
331 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer2D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
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());
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++)
350 MCAuto<MEDCouplingUMesh> vorTess;
352 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
353 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
358 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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())
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]);
374 MCAuto<MEDCouplingUMesh> a;
375 MCAuto<DataArrayInt> b,c;
377 DataArrayInt *bPtr(0),*cPtr(0);
378 a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
381 MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
382 if(part->getNumberOfTuples()!=1)
383 throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
384 MCAuto<MEDCouplingUMesh> newVorCell;
386 MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
387 newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
389 newVorCell->zipCoords();
390 MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
391 modifiedCell->zipCoords();
392 l0[poly]=modifiedCell;
394 MCAuto<DataArrayInt> ids;
396 DataArrayInt *tmp(0);
397 bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
400 throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
402 MCAuto<DataArrayDouble> newCoords;
404 MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
405 newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
407 const double *cPtr(newCoords->begin());
408 for(int i=0;i<newCoords->getNumberOfTuples();i++,cPtr+=2)
410 std::set<int> zeCandidates;
412 std::vector<int> zeCandidatesTmp;
413 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
414 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
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;
421 newVorCells.push_back(newVorCell);
423 l0.push_back(MergeVorCells(newVorCells,eps));
425 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
426 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
428 bool dummy; int dummy2;
429 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
434 MCAuto<MEDCouplingUMesh> Split3DCellInParts(const MEDCouplingUMesh *m, const double pt[3], const double seed[3], double eps, int tmp[2])
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));
444 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
446 double eps2(1.-sqrt(eps));// 2nd eps for interpolation. Here the eps is computed to feet cos(eps) ~ 1-eps^2
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());
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++)
462 MCAuto<MEDCouplingUMesh> vorTess;
464 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
465 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
470 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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;
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));
483 while(!elemsToDo.empty())
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]);
490 MCAuto<MEDCouplingUMesh> cells;
493 cells=Split3DCellInParts(tile,pt,seed,eps,tmp);
495 catch(INTERP_KERNEL::Exception& e)
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
506 MCAuto<MEDCouplingUMesh> faces2;
508 MCAuto<DataArrayInt> d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New());
509 faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I);
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++)
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++)
530 zeCandidates.insert(*it4);
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;
539 newVorCells.push_back(newVorCell);
542 l0.push_back(MergeVorCells3D(newVorCells,eps));
544 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
545 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
547 bool dummy; int dummy2;
548 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));