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 Voronizer2D::getDimension() const
41 int Voronizer3D::getDimension() const
46 MCAuto<MEDCouplingUMesh> ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector<double>& bbox, double eps)
48 static const double FACT=1.2;
49 MCAuto<MEDCouplingCMesh> m(MEDCouplingCMesh::New());
50 MCAuto<DataArrayDouble> arr1(DataArrayDouble::New()); arr1->alloc(2,1) ; arr1->setIJ(0,0,bbox[0]); arr1->setIJ(1,0,bbox[1]);
51 MCAuto<DataArrayDouble> arr2(DataArrayDouble::New()); arr2->alloc(2,1) ; arr2->setIJ(0,0,bbox[2]); arr2->setIJ(1,0,bbox[3]);
52 m->setCoords(arr1,arr2);
53 static const double PT[2]={0.,0.};
55 MCAuto<MEDCouplingUMesh> mu(m->buildUnstructured());
56 double l(std::max(bbox[1]-bbox[0],bbox[3]-bbox[2]));
57 double middle[2]={(pt1[0]+pt2[0])/2.,(pt1[1]+pt2[1])/2.};
58 double v[2]={pt1[0],pt1[1]};
59 DataArrayDouble::Rotate2DAlg(middle,M_PI/2,1,v,v);
60 v[0]=middle[0]-v[0]; v[1]=middle[1]-v[1];
62 double nor(sqrt(v[0]*v[0]+v[1]*v[1]));
65 MCAuto<MEDCouplingUMesh> line(MEDCouplingUMesh::New("line",1));
67 MCAuto<DataArrayDouble> coo(DataArrayDouble::New()); coo->alloc(2,2);
68 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]);
71 line->allocateCells();
72 static const int CONN[2]={0,1};
73 line->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN);
74 MCAuto<MEDCouplingUMesh> sp2,sp1;
76 DataArrayInt *cellNb1(0),*cellNb2(0);
77 MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0);
78 MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2);
80 MCAuto<DataArrayInt> cellNb10(cellNb1),cellNb20(cellNb2);
83 sp2->getCellsContainingPoint(pt1,eps,ccp);
85 throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !");
86 MCAuto<MEDCouplingUMesh> ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true));
92 MCAuto<MEDCouplingUMesh> MergeVorCells2D(MEDCouplingUMesh *p, double eps, bool isZip)
94 MCAuto<DataArrayInt> edgeToKeep;
95 MCAuto<MEDCouplingUMesh> p0;
97 MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
98 p0=p->buildDescendingConnectivity(d,di,rd,rdi);
99 MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
100 edgeToKeep=dsi->findIdsEqual(1);
102 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
105 skinOfRes->zipCoords();
106 if(skinOfRes->getNumberOfCells()!=skinOfRes->getNumberOfNodes())
107 throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !");
109 MCAuto<DataArrayInt> d(skinOfRes->orderConsecutiveCells1D());
110 MCAuto<MEDCoupling1SGTUMesh> skinOfRes2;
112 MCAuto<MEDCouplingUMesh> part(skinOfRes->buildPartOfMySelf(d->begin(),d->end()));
113 skinOfRes2=MEDCoupling1SGTUMesh::New(part);
115 MCAuto<DataArrayInt> c(skinOfRes2->getNodalConnectivity()->deepCopy());
116 c->circularPermutation(1);
118 std::vector< MCAuto<DataArrayInt> > vdi(c->explodeComponents());
119 if(!vdi[0]->isEqual(*vdi[1]))
120 throw INTERP_KERNEL::Exception("MergeVorCells : internal error !");
121 MCAuto<MEDCouplingUMesh> m(MEDCouplingUMesh::New("",2));
122 m->setCoords(skinOfRes2->getCoords());
124 m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin());
128 MCAuto<MEDCouplingUMesh> MergeVorCells(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
130 std::size_t sz(vcs.size());
132 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
135 MCAuto<MEDCouplingUMesh> p;
137 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
138 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
142 bool dummy; int dummy2;
143 MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
145 return MergeVorCells2D(p,eps,true);
148 MCAuto<MEDCouplingUMesh> MergeVorCells3D(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
150 std::size_t sz(vcs.size());
152 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
155 MCAuto<MEDCouplingUMesh> p;
157 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
158 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
162 bool dummy; int dummy2;
163 MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
165 MCAuto<DataArrayInt> edgeToKeep;
166 MCAuto<MEDCouplingUMesh> p0;
168 MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
169 p0=p->buildDescendingConnectivity(d,di,rd,rdi);
170 MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
171 edgeToKeep=dsi->findIdsEqual(1);
173 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
174 MCAuto<DataArrayDouble> eqn(skinOfRes->computePlaneEquationOf3DFaces());
175 MCAuto<DataArrayInt> comm,commI;
177 DataArrayInt *a(0),*b(0);
178 eqn->findCommonTuples(eps,0,a,b);
180 //comm=DataArrayInt::New(); comm->alloc(0,1); commI=DataArrayInt::New(); commI->alloc(1,1); commI->setIJ(0,0,0);
182 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",3));
183 ret->setCoords(skinOfRes->getCoords());
184 ret->allocateCells();
185 std::vector<int> conn;
187 for(int i=0;i<commI->getNumberOfTuples()-1;i++,jj++)
191 MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(comm->begin()+commI->getIJ(i,0),comm->begin()+commI->getIJ(i+1,0),true));
192 MCAuto<MEDCouplingUMesh> tmp2;
193 if(commI->getIJ(i+1,0)-commI->getIJ(i,0)==1)
196 tmp2=MergeVorCells2D(tmp,eps,false);
197 const int *cPtr(tmp2->getNodalConnectivity()->begin()),*ciPtr(tmp2->getNodalConnectivityIndex()->begin());
198 conn.insert(conn.end(),cPtr+1,cPtr+ciPtr[1]);
200 MCAuto<DataArrayInt> remain(comm->buildComplement(skinOfRes->getNumberOfCells()));
202 MCAuto<MEDCouplingUMesh> tmp(skinOfRes->buildPartOfMySelf(remain->begin(),remain->end(),true));
203 const int *cPtr(tmp->getNodalConnectivity()->begin()),*ciPtr(tmp->getNodalConnectivityIndex()->begin());
204 for(int i=0;i<remain->getNumberOfTuples();i++,jj++)
208 conn.insert(conn.end(),cPtr+ciPtr[i]+1,cPtr+ciPtr[i+1]);
211 ret->insertNextCell(INTERP_KERNEL::NORM_POLYHED,conn.size(),&conn[0]);
215 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer2D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
218 throw INTERP_KERNEL::Exception("Voronoize2D : null pointer !");
219 m->checkConsistencyLight();
220 points->checkAllocated();
221 if(m->getMeshDimension()!=2 || m->getSpaceDimension()!=2 || points->getNumberOfComponents()!=2)
222 throw INTERP_KERNEL::Exception("Voronoize2D : spacedim must be equal to 2 and meshdim also equal to 2 !");
223 if(m->getNumberOfCells()!=1)
224 throw INTERP_KERNEL::Exception("Voronoize2D : mesh is expected to have only one cell !");
225 int nbPts(points->getNumberOfTuples());
227 throw INTERP_KERNEL::Exception("Voronoize2D : at least one point expected !");
228 std::vector<double> bbox(4);
229 m->getBoundingBox(&bbox[0]);
230 std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
231 const double *pts(points->begin());
232 for(int i=1;i<nbPts;i++)
234 MCAuto<MEDCouplingUMesh> vorTess;
236 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
237 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
242 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
244 std::vector<int> polygsToIterOn;
245 const double *pt(pts+i*2);
246 vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
247 if(polygsToIterOn.size()<1)
248 throw INTERP_KERNEL::Exception("Voronoize2D : presence of a point outside the given cell !");
249 std::set<int> elemsToDo,elemsDone; elemsToDo.insert(polygsToIterOn[0]);
250 std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
251 while(!elemsToDo.empty())
253 int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
254 const double *seed(pts+2*poly);
255 MCAuto<MEDCouplingUMesh> cell(ComputeBigCellFrom(pt,seed,bbox,eps));
256 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
258 MCAuto<MEDCouplingUMesh> a;
259 MCAuto<DataArrayInt> b,c;
261 DataArrayInt *bPtr(0),*cPtr(0);
262 a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
265 MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
266 if(part->getNumberOfTuples()!=1)
267 throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
268 MCAuto<MEDCouplingUMesh> newVorCell;
270 MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
271 newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
273 newVorCell->zipCoords();
274 MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
275 modifiedCell->zipCoords();
276 l0[poly]=modifiedCell;
278 MCAuto<DataArrayInt> ids;
280 DataArrayInt *tmp(0);
281 bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
284 throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
286 MCAuto<DataArrayDouble> newCoords;
288 MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
289 newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
291 const double *cPtr(newCoords->begin());
292 for(int i=0;i<newCoords->getNumberOfTuples();i++,cPtr+=2)
294 std::set<int> zeCandidates;
296 std::vector<int> zeCandidatesTmp;
297 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
298 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
300 std::set<int> tmp,newElementsToDo;
301 std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
302 std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
303 elemsToDo=newElementsToDo;
305 newVorCells.push_back(newVorCell);
307 l0.push_back(MergeVorCells(newVorCells,eps));
309 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
310 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
312 bool dummy; int dummy2;
313 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));
318 MCAuto<MEDCouplingUMesh> Split3DCellInParts(const MEDCouplingUMesh *m, const double pt[3], const double seed[3], double eps, int tmp[2])
320 if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || m->getNumberOfCells()!=1)
321 throw INTERP_KERNEL::Exception("Split3DCellInParts : expecting a 3D with exactly one cell !");
322 double middle[3]={(pt[0]+seed[0])/2.,(pt[1]+seed[1])/2.,(pt[2]+seed[2])/2.};
323 double vec[3]={pt[0]-seed[0],pt[1]-seed[1],pt[2]-seed[2]};
324 MCAuto<MEDCouplingUMesh> res(m->clipSingle3DCellByPlane(middle,vec,eps));
328 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronizer3D::doIt(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps) const
330 double eps2(1.-sqrt(eps));// 2nd eps for interpolation. Here the eps is computed to feet cos(eps) ~ 1-eps^2
332 throw INTERP_KERNEL::Exception("Voronoize3D : null pointer !");
333 m->checkConsistencyLight();
334 points->checkAllocated();
335 if(m->getMeshDimension()!=3 || m->getSpaceDimension()!=3 || points->getNumberOfComponents()!=3)
336 throw INTERP_KERNEL::Exception("Voronoize3D : spacedim must be equal to 3 and meshdim also equal to 3 !");
337 if(m->getNumberOfCells()!=1)
338 throw INTERP_KERNEL::Exception("Voronoize3D : mesh is expected to have only one cell !");
339 int nbPts(points->getNumberOfTuples());
341 throw INTERP_KERNEL::Exception("Voronoize3D : at least one point expected !");
342 std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
343 const double *pts(points->begin());
344 for(int i=1;i<nbPts;i++)
346 MCAuto<MEDCouplingUMesh> vorTess;
348 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
349 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
354 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
356 std::vector<int> polygsToIterOn;
357 const double *pt(pts+i*3);
358 vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
359 if(polygsToIterOn.size()<1)
360 throw INTERP_KERNEL::Exception("Voronoize3D : presence of a point outside the given cell !");
361 std::set<int> elemsToDo(polygsToIterOn.begin(),polygsToIterOn.end()),elemsDone;
363 std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
364 MCAuto<DataArrayInt> d(DataArrayInt::New()),dI(DataArrayInt::New()),rd(DataArrayInt::New()),rdI(DataArrayInt::New());
365 MCAuto<MEDCouplingUMesh> faces(vorTess->buildDescendingConnectivity(d,dI,rd,rdI));
367 while(!elemsToDo.empty())
369 int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
370 const double *seed(pts+3*poly);
371 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
374 MCAuto<MEDCouplingUMesh> cells;
377 cells=Split3DCellInParts(tile,pt,seed,eps,tmp);
379 catch(INTERP_KERNEL::Exception& e)
383 MCAuto<MEDCouplingUMesh> newVorCell(cells->buildPartOfMySelfSlice(1,2,1,true));
384 newVorCell->zipCoords();
385 MCAuto<MEDCouplingUMesh> modifiedCell(cells->buildPartOfMySelfSlice(0,1,1,true));
386 modifiedCell->zipCoords();
387 l0[poly]=modifiedCell;
388 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
390 MCAuto<MEDCouplingUMesh> faces2;
392 MCAuto<DataArrayInt> d2(DataArrayInt::New()),d2I(DataArrayInt::New()),rd2(DataArrayInt::New()),rd2I(DataArrayInt::New());
393 faces2=newVorCell->buildDescendingConnectivity(d2,d2I,rd2,rd2I);
395 MCAuto<MEDCouplingUMesh> faces3(faces2->buildPartOfMySelfSlice(1,faces2->getNumberOfCells(),1,true));// suppress internal face
396 MCAuto<MEDCouplingUMesh> facesOfCurSplitPol(faces->buildPartOfMySelf(d->begin()+dI->getIJ(poly,0),d->begin()+dI->getIJ(poly+1,0),true));
397 // intersection between the out faces of newVorCell and the neighbor faces of poly polyhedron -> candidates
398 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(facesOfCurSplitPol);
399 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(faces3);
400 INTERP_KERNEL::Interpolation3DSurf interpolation;
401 interpolation.setMinDotBtwPlane3DSurfIntersect(eps2);
402 interpolation.setMaxDistance3DSurfIntersect(eps);
403 interpolation.setPrecision(1e-12);
404 std::vector<std::map<int,double> > matrix;
405 interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix,"P0P0");
406 std::set<int> zeCandidates;
407 for(std::vector<std::map<int,double> >::const_iterator it2=matrix.begin();it2!=matrix.end();it2++)
408 for(std::map<int,double>::const_iterator it3=(*it2).begin();it3!=(*it2).end();it3++)
410 int faceIdInVorTess(d->getIJ(dI->getIJ(poly,0)+(*it3).first,0));
411 for(const int *it4=rd->begin()+rdI->getIJ(faceIdInVorTess,0);it4!=rd->begin()+rdI->getIJ(faceIdInVorTess+1,0);it4++)
414 zeCandidates.insert(*it4);
417 std::set<int> tmp,newElementsToDo;
418 std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
419 std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
420 elemsToDo=newElementsToDo;
423 newVorCells.push_back(newVorCell);
426 l0.push_back(MergeVorCells3D(newVorCells,eps));
428 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
429 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));
431 bool dummy; int dummy2;
432 MCAuto<DataArrayInt> dummy3(ret->mergeNodes(eps,dummy,dummy2));