1 // Copyright (C) 2007-2023 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"
24 #include "MEDCouplingFieldDouble.hxx"
27 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
28 #include "Interpolation2D.txx"
29 #include "Interpolation3DSurf.hxx"
31 using namespace MEDCoupling;
33 Voronizer::~Voronizer()
37 int Voronizer1D::getDimension() const
42 int Voronizer2D::getDimension() const
47 int Voronizer3D::getDimension() const
52 MCAuto<MEDCouplingUMesh> ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector<double>& bbox, double eps)
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.};
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];
68 double nor(sqrt(v[0]*v[0]+v[1]*v[1]));
71 MCAuto<MEDCouplingUMesh> line(MEDCouplingUMesh::New("line",1));
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]);
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;
82 DataArrayIdType *cellNb1(0),*cellNb2(0);
83 MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0);
84 MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2);
86 MCAuto<DataArrayIdType> cellNb10(cellNb1),cellNb20(cellNb2);
88 std::vector<mcIdType> ccp;
89 sp2->getCellsContainingPoint(pt1,eps,ccp);
91 throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !");
92 MCAuto<MEDCouplingUMesh> ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true));
95 MCAuto<MEDCouplingFieldDouble> tmp(ret->getMeasureField(false));
96 if(tmp->getArray()->getIJ(0,0)<0)
97 ret->invertOrientationOfAllCells();
103 MCAuto<MEDCouplingUMesh> MergeVorCells2D(MEDCouplingUMesh *p, double eps, bool isZip)
105 MCAuto<DataArrayIdType> edgeToKeep;
106 MCAuto<MEDCouplingUMesh> p0;
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);
113 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
116 skinOfRes->zipCoords();
117 if(ToIdType(skinOfRes->getNumberOfCells())!=skinOfRes->getNumberOfNodes())
118 throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !");
120 MCAuto<DataArrayIdType> d(skinOfRes->orderConsecutiveCells1D());
121 MCAuto<MEDCoupling1SGTUMesh> skinOfRes2;
123 MCAuto<MEDCouplingUMesh> part(skinOfRes->buildPartOfMySelf(d->begin(),d->end()));
124 skinOfRes2=MEDCoupling1SGTUMesh::New(part);
126 MCAuto<DataArrayIdType> c(skinOfRes2->getNodalConnectivity()->deepCopy());
127 c->circularPermutation(1);
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());
135 m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin());
139 MCAuto<MEDCouplingUMesh> MergeVorCells(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
141 std::size_t sz(vcs.size());
143 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
146 MCAuto<MEDCouplingUMesh> p;
148 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
149 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
153 bool dummy; mcIdType dummy2;
154 MCAuto<DataArrayIdType> dummy3(p->mergeNodes(eps,dummy,dummy2));
156 return MergeVorCells2D(p,eps,true);
160 * suppress additional sub points on edges
162 MCAuto<MEDCouplingUMesh> SimplifyPolygon(const MEDCouplingUMesh *m, double eps)
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++)
172 mcIdType prev(conn[(i+nbPtsInPolygon-1)%nbPtsInPolygon+1]),current(conn[i%nbPtsInPolygon+1]),zeNext(conn[(i+1)%nbPtsInPolygon+1]);
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],
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],
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);
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]);
193 MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
195 std::size_t sz(vcs.size());
197 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
200 MCAuto<MEDCouplingUMesh> p;
202 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
203 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
207 bool dummy; mcIdType dummy2;
208 MCAuto<DataArrayIdType> dummy3(p->mergeNodes(eps,dummy,dummy2));
210 MCAuto<DataArrayIdType> edgeToKeep;
211 MCAuto<MEDCouplingUMesh> p0;
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);
218 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
219 MCAuto<DataArrayDouble> eqn(skinOfRes->computePlaneEquationOf3DFaces());
220 MCAuto<DataArrayIdType> comm,commI;
222 DataArrayIdType *a(0),*b(0);
223 eqn->findCommonTuples(eps,0,a,b);
225 //comm=DataArrayIdType::New(); comm->alloc(0,1); commI=DataArrayIdType::New(); commI->alloc(1,1); commI->setIJ(0,0,0);
227 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",3));
228 ret->setCoords(skinOfRes->getCoords());
229 ret->allocateCells();
230 std::vector<mcIdType> conn;
232 for(mcIdType i=0;i<commI->getNumberOfTuples()-1;i++,jj++)
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)
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]);
246 MCAuto<DataArrayIdType> remain(comm->buildComplement(ToIdType(skinOfRes->getNumberOfCells())));
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++)
254 conn.insert(conn.end(),cPtr+ciPtr[i]+1,cPtr+ciPtr[i+1]);
257 ret->insertNextCell(INTERP_KERNEL::NORM_POLYHED,ToIdType(conn.size()),&conn[0]);
261 MCAuto<MEDCouplingUMesh> MergeVorCells1D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
263 static const mcIdType CONN_SEG2_DFT[2]={0,1};
265 throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 1 !");
269 throw INTERP_KERNEL::Exception("MergeVorCells1D : internal error 2 !");
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]];
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]];
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);
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); }
290 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer1D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
292 static const mcIdType CONN_SEG2_DFT[2]={0,1};
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());
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++)
310 MCAuto<MEDCouplingUMesh> vorTess;
312 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
313 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
318 MCAuto<DataArrayIdType> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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++)
332 double seed(pts[poly]),zept(*pt);
333 double mid((seed+zept)/2.);
335 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
339 const mcIdType *connPtr(tile->getNodalConnectivity()->begin());
340 const double *coordPtr(tile->getCoords()->begin());
341 a=coordPtr[connPtr[1]]; b=coordPtr[connPtr[2]];
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);
346 if(t0->isMonotonic(true,eps))
347 { pol0[0]=a; pol0[1]=mid; pol1[0]=mid; pol1[1]=b; }
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);
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);
358 l0[poly]=modifiedCell;
359 newVorCells.push_back(newVorCell);
361 l0.push_back(MergeVorCells1D(newVorCells,eps));
363 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
364 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
366 bool dummy; mcIdType dummy2;
367 MCAuto<DataArrayIdType> dummy3(ret->mergeNodes(eps,dummy,dummy2));
372 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer2D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
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());
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++)
391 MCAuto<MEDCouplingUMesh> vorTess;
393 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
394 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
399 MCAuto<DataArrayIdType> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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())
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]);
415 MCAuto<MEDCouplingUMesh> a;
416 MCAuto<DataArrayIdType> b,c;
418 DataArrayIdType *bPtr(0),*cPtr(0);
419 a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
422 MCAuto<DataArrayIdType> part(c->findIdsEqual(-1));
423 if(part->getNumberOfTuples()!=1)
424 throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
425 MCAuto<MEDCouplingUMesh> newVorCell;
427 MCAuto<DataArrayIdType> tmp(part->buildComplement(ToIdType(a->getNumberOfCells())));
428 newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
430 newVorCell->zipCoords();
431 MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
432 modifiedCell->zipCoords();
434 MCAuto<MEDCouplingFieldDouble> tmp(modifiedCell->getMeasureField(false));
435 if(tmp->getArray()->getIJ(0,0)<0)
436 modifiedCell->invertOrientationOfAllCells();
438 l0[poly]=modifiedCell;
440 MCAuto<DataArrayIdType> ids;
442 DataArrayIdType *tmp(0);
443 bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
446 throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
448 MCAuto<DataArrayDouble> newCoords;
450 MCAuto<DataArrayIdType> tmp(ids->buildComplement(a->getNumberOfNodes()));
451 newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
453 const double *cPtr(newCoords->begin());
454 for(mcIdType j=0;j<newCoords->getNumberOfTuples();j++,cPtr+=2)
456 std::set<mcIdType> zeCandidates;
458 std::vector<mcIdType> zeCandidatesTmp;
459 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
460 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
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;
467 newVorCells.push_back(newVorCell);
469 MCAuto<MEDCouplingUMesh> mergedVorCell(MergeVorCells(newVorCells,eps));
471 MCAuto<MEDCouplingFieldDouble> tmp(mergedVorCell->getMeasureField(false));
472 if(tmp->getArray()->getIJ(0,0)<0)
473 mergedVorCell->invertOrientationOfAllCells();
475 l0.push_back(mergedVorCell);
477 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
478 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
480 bool dummy; mcIdType dummy2;
481 MCAuto<DataArrayIdType> dummy3(ret->mergeNodes(eps,dummy,dummy2));
486 MCAuto<MEDCouplingUMesh> Split3DCellInParts(const MEDCouplingUMesh *m, const double pt[3], const double seed[3], double eps, mcIdType tmp[2])
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));
496 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
498 //double eps2(1.-sqrt(eps));// 2nd eps for interpolation. Here the eps is computed to feet cos(eps) ~ 1-eps^2
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());
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++)
514 MCAuto<MEDCouplingUMesh> vorTess;
516 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
517 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
522 MCAuto<DataArrayIdType> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
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++)
532 const double *seed(pts+3*poly);
533 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
536 MCAuto<MEDCouplingUMesh> cells;
539 cells=Split3DCellInParts(tile,pt,seed,eps,tmp);
541 catch(INTERP_KERNEL::Exception&)
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;
552 l0.push_back(MergeVorCells3D(newVorCells,eps));
554 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
555 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
557 bool dummy; mcIdType dummy2;
558 MCAuto<DataArrayIdType> dummy3(ret->mergeNodes(eps,dummy,dummy2));