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 using namespace MEDCoupling;
28 MCAuto<MEDCouplingUMesh> ComputeBigCellFrom(const double pt1[2], const double pt2[2], const std::vector<double>& bbox, double eps)
30 static const double FACT=1.2;
31 MCAuto<MEDCouplingCMesh> m(MEDCouplingCMesh::New());
32 MCAuto<DataArrayDouble> arr1(DataArrayDouble::New()); arr1->alloc(2,1) ; arr1->setIJ(0,0,bbox[0]); arr1->setIJ(1,0,bbox[1]);
33 MCAuto<DataArrayDouble> arr2(DataArrayDouble::New()); arr2->alloc(2,1) ; arr2->setIJ(0,0,bbox[2]); arr2->setIJ(1,0,bbox[3]);
34 m->setCoords(arr1,arr2);
35 static const double PT[2]={0.,0.};
37 MCAuto<MEDCouplingUMesh> mu(m->buildUnstructured());
38 double l(std::max(bbox[1]-bbox[0],bbox[3]-bbox[2]));
39 double middle[2]={(pt1[0]+pt2[0])/2.,(pt1[1]+pt2[1])/2.};
40 double v[2]={pt1[0],pt1[1]};
41 DataArrayDouble::Rotate2DAlg(middle,M_PI/2,1,v,v);
42 v[0]=middle[0]-v[0]; v[1]=middle[1]-v[1];
44 double nor(sqrt(v[0]*v[0]+v[1]*v[1]));
47 MCAuto<MEDCouplingUMesh> line(MEDCouplingUMesh::New("line",1));
49 MCAuto<DataArrayDouble> coo(DataArrayDouble::New()); coo->alloc(2,2);
50 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]);
53 line->allocateCells();
54 static const int CONN[2]={0,1};
55 line->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,CONN);
56 MCAuto<MEDCouplingUMesh> sp2,sp1;
58 DataArrayInt *cellNb1(0),*cellNb2(0);
59 MEDCouplingUMesh *sp2Pt(0),*sp1Pt(0);
60 MEDCouplingUMesh::Intersect2DMeshWith1DLine(mu,line,eps,sp2Pt,sp1Pt,cellNb1,cellNb2);
62 MCAuto<DataArrayInt> cellNb10(cellNb1),cellNb20(cellNb2);
65 sp2->getCellsContainingPoint(pt1,eps,ccp);
67 throw INTERP_KERNEL::Exception("ComputeBigCellFrom : expected single element !");
68 MCAuto<MEDCouplingUMesh> ret(sp2->buildPartOfMySelfSlice(ccp[0],ccp[0]+1,1,true));
73 MCAuto<MEDCouplingUMesh> MergeVorCells(const std::vector< MCAuto<MEDCouplingUMesh> >& vcs, double eps)
75 std::size_t sz(vcs.size());
77 throw INTERP_KERNEL::Exception("MergeVorCells : len of input vec expected to be >= 1 !");
80 MCAuto<MEDCouplingUMesh> p;
82 std::vector< const MEDCouplingUMesh * > vcsBis(VecAutoToVecOfCstPt(vcs));
83 p=MEDCouplingUMesh::MergeUMeshes(vcsBis);
87 bool dummy; int dummy2;
88 MCAuto<DataArrayInt> dummy3(p->mergeNodes(eps,dummy,dummy2));
90 MCAuto<DataArrayInt> edgeToKeep;
91 MCAuto<MEDCouplingUMesh> p0;
93 MCAuto<DataArrayInt> d(DataArrayInt::New()),di(DataArrayInt::New()),rd(DataArrayInt::New()),rdi(DataArrayInt::New());
94 p0=p->buildDescendingConnectivity(d,di,rd,rdi);
95 MCAuto<DataArrayInt> dsi(rdi->deltaShiftIndex());
96 edgeToKeep=dsi->findIdsEqual(1);
98 MCAuto<MEDCouplingUMesh> skinOfRes(p0->buildPartOfMySelf(edgeToKeep->begin(),edgeToKeep->end()));
99 skinOfRes->zipCoords();
100 if(skinOfRes->getNumberOfCells()!=skinOfRes->getNumberOfNodes())
101 throw INTERP_KERNEL::Exception("MergeVorCells : result of merge looks bad !");
102 MCAuto<DataArrayInt> d(skinOfRes->orderConsecutiveCells1D());
103 MCAuto<MEDCoupling1SGTUMesh> skinOfRes2;
105 MCAuto<MEDCouplingUMesh> part(skinOfRes->buildPartOfMySelf(d->begin(),d->end()));
106 skinOfRes2=MEDCoupling1SGTUMesh::New(part);
108 MCAuto<DataArrayInt> c(skinOfRes2->getNodalConnectivity()->deepCopy());
109 c->circularPermutation(1);
111 std::vector< MCAuto<DataArrayInt> > vdi(c->explodeComponents());
112 if(!vdi[0]->isEqual(*vdi[1]))
113 throw INTERP_KERNEL::Exception("MergeVorCells : internal error !");
114 MCAuto<MEDCouplingUMesh> m(MEDCouplingUMesh::New("",2));
115 m->setCoords(skinOfRes2->getCoords());
117 m->insertNextCell(INTERP_KERNEL::NORM_POLYGON,vdi[0]->getNumberOfTuples(),vdi[0]->begin());
121 MCAuto<MEDCouplingUMesh> MEDCoupling::Voronoize2D(const MEDCouplingUMesh *m, const DataArrayDouble *points, double eps)
124 throw INTERP_KERNEL::Exception("Voronoize2D : null pointer !");
125 m->checkConsistencyLight();
126 points->checkAllocated();
127 if(m->getMeshDimension()!=2 || m->getSpaceDimension()!=2 || points->getNumberOfComponents()!=2)
128 throw INTERP_KERNEL::Exception("Voronoize2D : spacedim must be equal to 2 and meshdim also equal to 2 !");
129 if(m->getNumberOfCells()!=1)
130 throw INTERP_KERNEL::Exception("Voronoize2D : mesh is expected to have only one cell !");
131 int nbPts(points->getNumberOfTuples());
133 throw INTERP_KERNEL::Exception("Voronoize2D : at least one point expected !");
134 std::vector<double> bbox(4);
135 m->getBoundingBox(&bbox[0]);
136 std::vector< MCAuto<MEDCouplingUMesh> > l0(1,MCAuto<MEDCouplingUMesh>(m->deepCopy()));
137 const double *pts(points->begin());
138 for(int i=1;i<nbPts;i++)
140 MCAuto<MEDCouplingUMesh> vorTess;
142 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
143 vorTess=MEDCouplingUMesh::MergeUMeshes(l0Bis);
148 MCAuto<DataArrayInt> dummy3(vorTess->mergeNodes(eps,dummy,newNbNodes));
150 std::vector<int> polygsToIterOn;
151 const double *pt(pts+i*2);
152 vorTess->getCellsContainingPoint(pt,eps,polygsToIterOn);
153 if(polygsToIterOn.size()<1)
154 throw INTERP_KERNEL::Exception("Voronoize2D : presence of a point outside the given cell !");
155 std::set<int> elemsToDo,elemsDone; elemsToDo.insert(polygsToIterOn[0]);
156 std::vector< MCAuto<MEDCouplingUMesh> > newVorCells;
157 while(!elemsToDo.empty())
159 int poly(*elemsToDo.begin()); elemsToDo.erase(elemsToDo.begin()); elemsDone.insert(poly);
160 const double *seed(pts+2*poly);
161 MCAuto<MEDCouplingUMesh> cell(ComputeBigCellFrom(pt,seed,bbox,eps));
162 MCAuto<MEDCouplingUMesh> tile(l0[poly]);
164 MCAuto<MEDCouplingUMesh> a;
165 MCAuto<DataArrayInt> b,c;
167 DataArrayInt *bPtr(0),*cPtr(0);
168 a=MEDCouplingUMesh::Intersect2DMeshes(tile,cell,eps,bPtr,cPtr);
171 MCAuto<DataArrayInt> part(c->findIdsEqual(-1));
172 if(part->getNumberOfTuples()!=1)
173 throw INTERP_KERNEL::Exception("Voronoize2D : internal error");
174 MCAuto<MEDCouplingUMesh> newVorCell;
176 MCAuto<DataArrayInt> tmp(part->buildComplement(a->getNumberOfCells()));
177 newVorCell=a->buildPartOfMySelf(tmp->begin(),tmp->end());
179 newVorCell->zipCoords();
180 MCAuto<MEDCouplingUMesh> modifiedCell(a->buildPartOfMySelf(part->begin(),part->end()));
181 modifiedCell->zipCoords();
182 l0[poly]=modifiedCell;
184 MCAuto<DataArrayInt> ids;
186 DataArrayInt *tmp(0);
187 bool sta(a->getCoords()->areIncludedInMe(cell->getCoords(),eps,tmp));
190 throw INTERP_KERNEL::Exception("Voronoize2D : internal error 2 !");
192 MCAuto<DataArrayDouble> newCoords;
194 MCAuto<DataArrayInt> tmp(ids->buildComplement(a->getNumberOfNodes()));
195 newCoords=a->getCoords()->selectByTupleId(tmp->begin(),tmp->end());
197 const double *cPtr(newCoords->begin());
198 for(int i=0;i<newCoords->getNumberOfTuples();i++,cPtr+=2)
200 std::set<int> zeCandidates;
202 std::vector<int> zeCandidatesTmp;
203 vorTess->getCellsContainingPoint(cPtr,eps,zeCandidatesTmp);
204 zeCandidates.insert(zeCandidatesTmp.begin(),zeCandidatesTmp.end());
206 std::set<int> tmp,newElementsToDo;
207 std::set_difference(zeCandidates.begin(),zeCandidates.end(),elemsDone.begin(),elemsDone.end(),std::inserter(tmp,tmp.begin()));
208 std::set_union(elemsToDo.begin(),elemsToDo.end(),tmp.begin(),tmp.end(),std::inserter(newElementsToDo,newElementsToDo.begin()));
209 elemsToDo=newElementsToDo;
211 newVorCells.push_back(newVorCell);
213 l0.push_back(MergeVorCells(newVorCells,eps));
215 std::vector< const MEDCouplingUMesh * > l0Bis(VecAutoToVecOfCstPt(l0));
216 MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::MergeUMeshes(l0Bis));