1 // Copyright (C) 2007-2020 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 (CEA/DEN)
21 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
22 #include "InterpKernelGeo2DElementaryEdge.hxx"
23 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
24 #include "InterpKernelGeo2DAbstractEdge.hxx"
25 #include "InterpKernelGeo2DEdgeLin.hxx"
26 #include "InterpKernelGeo2DBounds.hxx"
27 #include "InterpKernelGeo2DEdge.txx"
29 #include "NormalizedUnstructuredMesh.hxx"
37 using namespace INTERP_KERNEL;
39 namespace INTERP_KERNEL
41 const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
44 QuadraticPolygon::QuadraticPolygon(const char *file)
46 char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
47 std::ifstream stream(file);
48 stream.exceptions(std::ios_base::eofbit);
52 stream.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
53 while(strcmp(currentLine,"1200 2")!=0);
56 Edge *newEdge=Edge::BuildFromXfigLine(stream);
58 newEdge->changeStartNodeWith(back()->getEndNode());
63 catch(const std::ifstream::failure&)
66 catch(const std::exception & ex)
68 // Some code before this catch throws the C++98 version of the exception (mangled
69 // name is " NSt8ios_base7failureE"), but FED24 compilation of the current version of the code
70 // tries to catch the C++11 version of it (mangled name "NSt8ios_base7failureB5cxx11E").
71 // So we have this nasty hack to catch both versions ...
73 // TODO: the below should be replaced by a better handling avoiding exception throwing.
74 if (std::string(ex.what()) == "basic_ios::clear")
76 //std::cout << "std::ios_base::failure C++11\n";
81 front()->changeStartNodeWith(back()->getEndNode());
84 QuadraticPolygon::~QuadraticPolygon()
88 QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes)
90 QuadraticPolygon *ret(new QuadraticPolygon);
91 std::size_t size=nodes.size();
92 for(std::size_t i=0;i<size;i++)
94 ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%size]));
100 QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes)
102 QuadraticPolygon *ret(new QuadraticPolygon);
103 std::size_t size=nodes.size();
104 for(std::size_t i=0;i<size/2;i++)
108 e1=new EdgeLin(nodes[i],nodes[i+size/2]);
109 e2=new EdgeLin(nodes[i+size/2],nodes[(i+1)%(size/2)]);
110 SegSegIntersector inters(*e1,*e2);
111 bool colinearity=inters.areColinears();
112 delete e1; delete e2;
114 ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%(size/2)]));
116 ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+size/2],nodes[(i+1)%(size/2)]));
117 nodes[i]->decrRef(); nodes[i+size/2]->decrRef();
122 Edge *QuadraticPolygon::BuildLinearEdge(std::vector<Node *>& nodes)
125 throw INTERP_KERNEL::Exception("QuadraticPolygon::BuildLinearEdge : input vector is expected to be of size 2 !");
126 Edge *ret(new EdgeLin(nodes[0],nodes[1]));
127 nodes[0]->decrRef(); nodes[1]->decrRef();
131 Edge *QuadraticPolygon::BuildArcCircleEdge(std::vector<Node *>& nodes)
134 throw INTERP_KERNEL::Exception("QuadraticPolygon::BuildArcCircleEdge : input vector is expected to be of size 3 !");
135 EdgeLin *e1(new EdgeLin(nodes[0],nodes[2])),*e2(new EdgeLin(nodes[2],nodes[1]));
136 SegSegIntersector inters(*e1,*e2);
137 bool colinearity=inters.areColinears();
138 delete e1; delete e2;
141 ret=new EdgeLin(nodes[0],nodes[1]);
143 ret=new EdgeArcCircle(nodes[0],nodes[2],nodes[1]);
144 nodes[0]->decrRef(); nodes[1]->decrRef(); nodes[2]->decrRef();
148 void QuadraticPolygon::BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName)
150 std::ofstream file(fileName);
151 file << std::setprecision(16);
152 file << " double coords[]=" << std::endl << " { ";
153 for(std::vector<Node *>::const_iterator iter=nodes.begin();iter!=nodes.end();iter++)
155 if(iter!=nodes.begin())
156 file << "," << std::endl << " ";
157 file << (*(*iter))[0] << ", " << (*(*iter))[1];
159 file << "};" << std::endl;
162 void QuadraticPolygon::closeMe() const
164 if(!front()->changeStartNodeWith(back()->getEndNode()))
165 throw(Exception("big error: not closed polygon..."));
168 void QuadraticPolygon::circularPermute()
170 if(_sub_edges.size()>1)
172 ElementaryEdge *first=_sub_edges.front();
173 _sub_edges.pop_front();
174 _sub_edges.push_back(first);
178 bool QuadraticPolygon::isButterflyAbs()
180 INTERP_KERNEL::Bounds b;
182 b.prepareForAggregation();
184 double dimChar=b.getCaracteristicDim();
185 b.getBarycenter(xBary,yBary);
186 applyGlobalSimilarity(xBary,yBary,dimChar);
188 return isButterfly();
191 bool QuadraticPolygon::isButterfly() const
193 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
195 Edge *e1=(*it)->getPtr();
196 std::list<ElementaryEdge *>::const_iterator it2=it;
198 for(;it2!=_sub_edges.end();it2++)
200 MergePoints commonNode;
201 ComposedEdge *outVal1=new ComposedEdge;
202 ComposedEdge *outVal2=new ComposedEdge;
203 Edge *e2=(*it2)->getPtr();
204 if(e1->intersectWith(e2,commonNode,*outVal1,*outVal2))
217 void QuadraticPolygon::dumpInXfigFileWithOther(const ComposedEdge& other, const char *fileName) const
219 std::ofstream file(fileName);
220 const int resolution=1200;
222 box.prepareForAggregation();
224 other.fillBounds(box);
225 dumpInXfigFile(file,resolution,box);
226 other.ComposedEdge::dumpInXfigFile(file,resolution,box);
229 void QuadraticPolygon::dumpInXfigFile(const char *fileName) const
231 std::ofstream file(fileName);
232 const int resolution=1200;
234 box.prepareForAggregation();
236 dumpInXfigFile(file,resolution,box);
239 void QuadraticPolygon::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const
241 stream << "#FIG 3.2 Produced by xfig version 3.2.5-alpha5" << std::endl;
242 stream << "Landscape" << std::endl;
243 stream << "Center" << std::endl;
244 stream << "Metric" << std::endl;
245 stream << "Letter" << std::endl;
246 stream << "100.00" << std::endl;
247 stream << "Single" << std::endl;
248 stream << "-2" << std::endl;
249 stream << resolution << " 2" << std::endl;
250 ComposedEdge::dumpInXfigFile(stream,resolution,box);
254 * Warning contrary to intersectWith method this method is \b NOT const. 'this' and 'other' are modified after call of this method.
256 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other)
258 double ret=0.,xBaryBB,yBaryBB;
259 double fact=normalize(&other,xBaryBB,yBaryBB);
260 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
261 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
263 ret+=fabs((*iter)->getArea());
266 return ret*fact*fact;
270 * This method splits 'this' with 'other' into smaller pieces localizable. 'mapThis' is a map that gives the correspondence
271 * between nodes contained in 'this' and node ids in a global mesh.
272 * In the same way, 'mapOther' gives the correspondence between nodes contained in 'other' and node ids in a
273 * global mesh from which 'other' is extracted.
274 * This method has 1 out parameter : 'edgesThis', After the call of this method, it contains the nodal connectivity (including type)
275 * of 'this' into globlal "this mesh".
276 * This method has 2 in/out parameters : 'subDivOther' and 'addCoo'.'otherEdgeIds' is useful to put values in
277 * 'edgesThis', 'subDivOther' and 'addCoo'.
278 * Size of 'otherEdgeIds' has to be equal to number of ElementaryEdges in 'other'. No check of that will be done.
279 * The term 'abs' in the name recalls that we normalize the mesh (spatially) so that node coordinates fit into [0;1].
280 * @param offset1 is the number of nodes contained in global mesh from which 'this' is extracted.
281 * @param offset2 is the sum of nodes contained in global mesh from which 'this' is extracted and 'other' is extracted.
282 * @param edgesInOtherColinearWithThis will be appended at the end of the vector with colinear edge ids of other (if any)
283 * @param otherEdgeIds is a vector with the same size than other before calling this method. It gives in the same order
284 * the cell id in global other mesh.
286 void QuadraticPolygon::splitAbs(QuadraticPolygon& other,
287 const std::map<INTERP_KERNEL::Node *,mcIdType>& mapThis, const std::map<INTERP_KERNEL::Node *,mcIdType>& mapOther,
288 mcIdType offset1, mcIdType offset2 ,
289 const std::vector<mcIdType>& otherEdgeIds,
290 std::vector<mcIdType>& edgesThis, mcIdType cellIdThis,
291 std::vector< std::vector<mcIdType> >& edgesInOtherColinearWithThis, std::vector< std::vector<mcIdType> >& subDivOther,
292 std::vector<double>& addCoo, std::map<mcIdType,mcIdType>& mergedNodes)
294 double xBaryBB, yBaryBB;
295 double fact=normalizeExt(&other, xBaryBB, yBaryBB);
298 IteratorOnComposedEdge itThis(this),itOther(&other); // other is (part of) the tool mesh
300 ComposedEdge *cThis=new ComposedEdge;
301 ComposedEdge *cOther=new ComposedEdge;
303 std::map<INTERP_KERNEL::Node *,mcIdType> mapAddCoo;
304 for(itOther.first();!itOther.finished();itOther.next(),i++)
306 // For each edge of 'other', proceed with intersections: the edge might split into sub-edges, 'otherTmp' will hold the final split result.
307 // In the process of going through all edges of 'other', 'this' (which contains initially only one edge)
308 // is sub-divised into several edges : each of them has to be tested when intersecting the next candidate stored in 'other'.
309 QuadraticPolygon otherTmp;
310 ElementaryEdge* curOther=itOther.current();
311 otherTmp.pushBack(new ElementaryEdge(curOther->getPtr(),curOther->getDirection())); curOther->getPtr()->incrRef();
312 IteratorOnComposedEdge itOtherTmp(&otherTmp);
313 for(itOtherTmp.first();!itOtherTmp.finished();itOtherTmp.next())
315 ElementaryEdge* curOtherTmp=itOtherTmp.current();
316 if(!curOtherTmp->isThereStartPoint())
317 itThis.first(); // reset iterator on 'this'
319 itThis=curOtherTmp->getIterator();
320 for(;!itThis.finished();)
322 ElementaryEdge* curThis=itThis.current();
325 std::map<INTERP_KERNEL::Node *,mcIdType>::const_iterator thisStart(mapThis.find(curThis->getStartNode())),thisEnd(mapThis.find(curThis->getEndNode())),
326 otherStart(mapOther.find(curOtherTmp->getStartNode())),otherEnd(mapOther.find(curOtherTmp->getEndNode()));
327 mcIdType thisStart2(thisStart==mapThis.end()?-1:(*thisStart).second), thisEnd2(thisEnd==mapThis.end()?-1:(*thisEnd).second),
328 otherStart2(otherStart==mapOther.end()?-1:(*otherStart).second+offset1),otherEnd2(otherEnd==mapOther.end()?-1:(*otherEnd).second+offset1);
330 if(curThis->getPtr()->intersectWith(curOtherTmp->getPtr(),merge,*cThis,*cOther))
332 if(!curThis->getDirection()) cThis->reverse();
333 if(!curOtherTmp->getDirection()) cOther->reverse();
334 // Substitution of a single simple edge by two sub-edges resulting from the intersection
335 // First modify the edges currently pointed by itThis and itOtherTmp so that the newly created node
336 // becomes the end of the previous sub-edge and the beginning of the next one.
337 UpdateNeighbours(merge,itThis,itOtherTmp,cThis,cOther);
338 delete curThis; // <-- destroying simple edge coming from pol1
339 delete curOtherTmp; // <-- destroying simple edge coming from pol2
340 // Then insert second part of the intersection.
341 itThis.insertElemEdges(cThis,true); // <-- 2nd param is true to go next.
342 itOtherTmp.insertElemEdges(cOther,false); // <-- 2nd param is false to avoid to go next.
343 curOtherTmp=itOtherTmp.current();
345 itThis.assignMySelfToAllElems(cOther);
348 cThis=new ComposedEdge;
349 cOther=new ComposedEdge;
353 UpdateNeighbours(merge,itThis,itOtherTmp,curThis,curOtherTmp);
356 merge.updateMergedNodes(thisStart2,thisEnd2,otherStart2,otherEnd2,mergedNodes);
359 // If one sub-edge of otherTmp is "ON" an edge of this, then we have colinearity (all edges in otherTmp are //)
360 if(otherTmp.presenceOfOn())
361 edgesInOtherColinearWithThis[otherEdgeIds[i]].push_back(cellIdThis);
362 // Converting back to integer connectivity:
363 if(otherTmp._sub_edges.size()>1) // only if a new point has been added (i.e. an actual intersection was done)
365 std::size_t jj = 0, sz(otherTmp._sub_edges.size());
366 for(std::list<ElementaryEdge *>::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++, jj++)
368 short skipStartOrEnd = jj == 0 ? -1 : (jj == sz-1 ? 1 : 0); // -1 means START, 1 means END, 0 other
369 (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,
370 fact,xBaryBB,yBaryBB, skipStartOrEnd,
371 /*out*/ subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo);
378 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
379 (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/edgesThis,addCoo,mapAddCoo);
384 * This method builds 'this' from its descending conn stored in crude mode (MEDCoupling).
385 * Descending conn is in FORTRAN relative mode in order to give the
386 * orientation of edge (see buildDescendingConnectivity2() method).
387 * See appendEdgeFromCrudeDataArray() for params description.
389 void QuadraticPolygon::buildFromCrudeDataArray(const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp, bool isQuad, const mcIdType *nodalBg, const double *coords,
390 const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges)
392 std::size_t nbOfSeg=std::distance(descBg,descEnd);
393 for(std::size_t i=0;i<nbOfSeg;i++)
395 appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
399 void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp, bool isQuad,
400 const mcIdType *nodalBg, const double *coords,
401 const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges)
405 bool direct=descBg[edgePos]>0;
406 mcIdType edgeId=std::abs(descBg[edgePos])-1; // back to C indexing mode
407 const std::vector<mcIdType>& subEdge=intersectEdges[edgeId];
408 std::size_t nbOfSubEdges=subEdge.size()/2;
409 for(std::size_t j=0;j<nbOfSubEdges;j++)
410 appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
414 std::size_t nbOfSeg=std::distance(descBg,descEnd);
415 const double *st=coords+2*(nodalBg[edgePos]);
416 INTERP_KERNEL::Node *st0=new INTERP_KERNEL::Node(st[0],st[1]);
417 const double *endd=coords+2*(nodalBg[(edgePos+1)%nbOfSeg]);
418 INTERP_KERNEL::Node *endd0=new INTERP_KERNEL::Node(endd[0],endd[1]);
419 const double *middle=coords+2*(nodalBg[edgePos+nbOfSeg]);
420 INTERP_KERNEL::Node *middle0=new INTERP_KERNEL::Node(middle[0],middle[1]);
422 e1=new EdgeLin(st0,middle0);
423 e2=new EdgeLin(middle0,endd0);
424 SegSegIntersector inters(*e1,*e2);
425 bool colinearity=inters.areColinears();
426 delete e1; delete e2;
428 bool direct=descBg[edgePos]>0;
429 mcIdType edgeId=std::abs(descBg[edgePos])-1;
430 const std::vector<mcIdType>& subEdge=intersectEdges[edgeId];
431 std::size_t nbOfSubEdges=subEdge.size()/2;
434 for(std::size_t j=0;j<nbOfSubEdges;j++)
435 appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
439 Edge *e=new EdgeArcCircle(st0,middle0,endd0,true);
440 for(std::size_t j=0;j<nbOfSubEdges;j++)
441 appendSubEdgeFromCrudeDataArray(e,j,direct,edgeId,subEdge,mapp);
444 st0->decrRef(); endd0->decrRef(); middle0->decrRef();
448 void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, mcIdType edgeId, const std::vector<mcIdType>& subEdge, const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp)
450 std::size_t nbOfSubEdges=subEdge.size()/2;
452 {//it is not a quadratic subedge
453 Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
454 Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
455 ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
459 {//it is a quadratic subedge
460 Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
461 Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
462 Edge *ee=baseEdge->buildEdgeLyingOnMe(start,end);
463 ElementaryEdge *eee=new ElementaryEdge(ee,true);
469 * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the
470 * orientation of edge.
472 void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp, bool isQuad, const mcIdType *nodalBg, const double *coords, const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges2,
473 const INTERP_KERNEL::QuadraticPolygon& pol1, const mcIdType *descBg1, const mcIdType *descEnd1, const std::vector<std::vector<mcIdType> >& intersectEdges1,
474 const std::vector< std::vector<mcIdType> >& colinear1,
475 std::map<mcIdType,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
477 std::size_t nbOfSeg=std::distance(descBg,descEnd);
478 for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
480 bool direct=descBg[i]>0;
481 mcIdType edgeId=std::abs(descBg[i])-1;//current edge id of pol2
482 std::map<mcIdType,std::vector<INTERP_KERNEL::ElementaryEdge *> >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]);
483 if(it1!=alreadyExistingIn2.end() || it2!=alreadyExistingIn2.end())
485 bool sameDir=(it1!=alreadyExistingIn2.end());
486 const std::vector<INTERP_KERNEL::ElementaryEdge *>& edgesAlreadyBuilt=sameDir?(*it1).second:(*it2).second;
489 for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_iterator it3=edgesAlreadyBuilt.begin();it3!=edgesAlreadyBuilt.end();it3++)
491 Edge *ee=(*it3)->getPtr(); ee->incrRef();
492 pushBack(new ElementaryEdge(ee,(*it3)->getDirection()));
497 for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_reverse_iterator it4=edgesAlreadyBuilt.rbegin();it4!=edgesAlreadyBuilt.rend();it4++)
499 Edge *ee=(*it4)->getPtr(); ee->incrRef();
500 pushBack(new ElementaryEdge(ee,!(*it4)->getDirection()));
505 bool directos=colinear1[edgeId].empty();
506 std::vector<std::pair<mcIdType,std::pair<bool,mcIdType> > > idIns1;
509 {// if the current edge of pol2 has one or more colinear edges part into pol1
510 const std::vector<mcIdType>& c=colinear1[edgeId];
511 std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
512 for(std::size_t j=0;j<nbOfEdgesIn1;j++)
514 mcIdType edgeId1=std::abs(descBg1[j])-1;
515 if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
517 idIns1.push_back(std::pair<mcIdType,std::pair<bool,mcIdType> >(edgeId1,std::pair<bool,mcIdType>(descBg1[j]>0,offset1)));// it exists an edge into pol1 given by tuple (idIn1,direct1) that is colinear at edge 'edgeId' in pol2
518 //std::pair<edgeId1); direct1=descBg1[j]>0;
520 offset1+=ToIdType(intersectEdges1[edgeId1].size()/2);//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
522 directos=idIns1.empty();
525 {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
526 std::size_t oldSz=_sub_edges.size();
527 appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges2);
528 std::size_t newSz=_sub_edges.size();
529 std::size_t zeSz=newSz-oldSz;
530 alreadyExistingIn2[descBg[i]].resize(zeSz);
531 std::list<ElementaryEdge *>::const_reverse_iterator it5=_sub_edges.rbegin();
532 for(std::size_t p=0;p<zeSz;p++,it5++)
533 alreadyExistingIn2[descBg[i]][zeSz-p-1]=*it5;
536 {//there is subpart of edge 'edgeId' of pol2 inside pol1
537 const std::vector<mcIdType>& subEdge=intersectEdges2[edgeId];
538 std::size_t nbOfSubEdges=subEdge.size()/2;
539 for(std::size_t j=0;j<nbOfSubEdges;j++)
541 mcIdType idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
542 mcIdType idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
543 bool direction11,found=false;
544 bool direct1;//store if needed the direction in 1
546 mcIdType nbOfSubEdges1;
547 for(std::vector<std::pair<mcIdType,std::pair<bool,mcIdType> > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++)
549 mcIdType idIn1=(*it).first;//store if needed the cell id in 1
550 direct1=(*it).second.first;
551 offset1=(*it).second.second;
552 const std::vector<mcIdType>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
553 nbOfSubEdges1=ToIdType(subEdge1PossiblyAlreadyIn1.size()/2);
555 for(mcIdType k=0;k<nbOfSubEdges1 && !found;k++)
556 {//perform a loop on all subedges of pol1 that includes edge 'edgeId' of pol2. For the moment we iterate only on subedges of ['idIn1']... To improve
557 if(subEdge1PossiblyAlreadyIn1[2*k]==idBg && subEdge1PossiblyAlreadyIn1[2*k+1]==idEnd)
558 { direction11=true; found=true; }
559 else if(subEdge1PossiblyAlreadyIn1[2*k]==idEnd && subEdge1PossiblyAlreadyIn1[2*k+1]==idBg)
560 { direction11=false; found=true; }
566 {//the current subedge of edge 'edgeId' of pol2 is not a part of the colinear edge 'idIn1' of pol1 -> build new Edge instance
567 //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
568 Node *start=(*mapp.find(idBg)).second;
569 Node *end=(*mapp.find(idEnd)).second;
570 ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
572 alreadyExistingIn2[descBg[i]].push_back(e);
575 {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1
576 ElementaryEdge *e=pol1[FromIdType<int>(offset1+(direct1?offset2:nbOfSubEdges1-offset2-1))];
577 Edge *ee=e->getPtr();
579 ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11));
581 alreadyExistingIn2[descBg[i]].push_back(e2);
589 * Method expected to be called on pol2. Every params not suffixed by numbered are supposed to refer to pol2 (this).
590 * Method to find edges that are ON.
592 void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges,
593 const INTERP_KERNEL::QuadraticPolygon& pol1, const mcIdType *descBg1, const mcIdType *descEnd1,
594 const std::vector<std::vector<mcIdType> >& intersectEdges1, const std::vector< std::vector<mcIdType> >& colinear1) const
596 std::size_t nbOfSeg=std::distance(descBg,descEnd);
597 for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
599 bool direct=descBg[i]>0;
600 mcIdType edgeId=std::abs(descBg[i])-1;//current edge id of pol2
601 const std::vector<mcIdType>& c=colinear1[edgeId];
604 const std::vector<mcIdType>& subEdge=intersectEdges[edgeId];
605 std::size_t nbOfSubEdges=subEdge.size()/2;
607 std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
609 for(std::size_t j=0;j<nbOfEdgesIn1;j++)
611 mcIdType edgeId1=std::abs(descBg1[j])-1;
612 if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
614 for(std::size_t k=0;k<nbOfSubEdges;k++)
616 mcIdType idBg=direct?subEdge[2*k]:subEdge[2*nbOfSubEdges-2*k-1];
617 mcIdType idEnd=direct?subEdge[2*k+1]:subEdge[2*nbOfSubEdges-2*k-2];
618 mcIdType idIn1=edgeId1;
619 bool direct1=descBg1[j]>0;
620 const std::vector<mcIdType>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
621 mcIdType nbOfSubEdges1=ToIdType(subEdge1PossiblyAlreadyIn1.size()/2);
624 for(mcIdType kk=0;kk<nbOfSubEdges1 && !found;kk++)
626 found=(subEdge1PossiblyAlreadyIn1[2*kk]==idBg && subEdge1PossiblyAlreadyIn1[2*kk+1]==idEnd) || (subEdge1PossiblyAlreadyIn1[2*kk]==idEnd && subEdge1PossiblyAlreadyIn1[2*kk+1]==idBg);
632 ElementaryEdge *e=pol1[(int)(offset1+(direct1?offset2:nbOfSubEdges1-offset2-1))];
633 e->getPtr()->declareOn();
637 offset1+=ToIdType(intersectEdges1[edgeId1].size()/2);//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
642 void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,mcIdType>& mapp, double xBary, double yBary, double fact, mcIdType offset, std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& conn, std::vector<mcIdType>& connI) const
645 bool presenceOfQuadratic=presenceOfQuadraticEdge();
646 conn.push_back(presenceOfQuadratic?NORM_QPOLYG:NORM_POLYGON);
647 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
650 tmp=(*it)->getStartNode();
651 std::map<INTERP_KERNEL::Node *,mcIdType>::const_iterator it1=mapp.find(tmp);
652 conn.push_back((*it1).second);
655 if(presenceOfQuadratic)
658 mcIdType off=offset+ToIdType(addCoordsQuadratic.size())/2;
659 for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++)
661 INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf();
662 node->unApplySimilarity(xBary,yBary,fact);
663 addCoordsQuadratic.push_back((*node)[0]);
664 addCoordsQuadratic.push_back((*node)[1]);
665 conn.push_back(off+j);
669 connI.push_back(connI.back()+nbOfNodesInPg+1);
673 * This method make the hypothesis that \a this and \a other are split at the minimum into edges that are fully IN, OUT or ON.
674 * This method returns newly created polygons in \a conn and \a connI and the corresponding ids ( \a idThis, \a idOther) are stored respectively into \a nbThis and \a nbOther.
675 * @param [in,out] edgesThis, parameter that keep informed the caller about the edges in this not shared by the result of intersection of \a this with \a other
676 * @param [in,out] edgesBoundaryOther, parameter that stores all edges in result of intersection that are not
678 void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTERP_KERNEL::Edge *>& edgesThis, std::set<INTERP_KERNEL::Edge *>& edgesBoundaryOther,
679 const std::map<INTERP_KERNEL::Node *,mcIdType>& mapp, mcIdType idThis, mcIdType idOther, mcIdType offset,
680 std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& conn, std::vector<mcIdType>& connI,
681 std::vector<mcIdType>& nbThis, std::vector<mcIdType>& nbOther)
683 double xBaryBB, yBaryBB;
684 double fact=normalizeExt(&other, xBaryBB, yBaryBB);
685 //Locate \a this relative to \a other (edges of \a this, aka \a pol1 are marked as IN or OUT)
686 other.performLocatingOperationSlow(*this); // without any assumption
687 std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(*this,other);
688 for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
690 (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI);
691 INTERP_KERNEL::IteratorOnComposedEdge it1(*it);
692 for(it1.first();!it1.finished();it1.next())
694 Edge *e=it1.current()->getPtr();
695 if(edgesThis.find(e)!=edgesThis.end())
699 if(edgesBoundaryOther.find(e)!=edgesBoundaryOther.end())
700 edgesBoundaryOther.erase(e);
702 edgesBoundaryOther.insert(e);
705 nbThis.push_back(idThis);
706 nbOther.push_back(idOther);
709 unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
713 * Remove the two 'identical' edges from the list, and drecRef the objects.
715 void QuadraticPolygon::cleanDegeneratedConsecutiveEdges()
717 IteratorOnComposedEdge it(this);
718 ElementaryEdge * prevEdge = 0;
719 if (recursiveSize() > 2)
720 for(it.first();!it.finished();it.next())
722 ElementaryEdge * cur = it.current();
723 if (prevEdge && prevEdge->hasSameExtremities(*cur))
727 prevEdge = it.current();
735 * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method.
736 * 'other' is a QuadraticPolygon of \b non closed edges.
738 double QuadraticPolygon::intersectWithAbs1D(QuadraticPolygon& other, bool& isColinear)
740 double ret = 0., xBaryBB, yBaryBB;
741 double fact = normalize(&other, xBaryBB, yBaryBB);
743 QuadraticPolygon cpyOfThis(*this);
744 QuadraticPolygon cpyOfOther(other);
746 SplitPolygonsEachOther(cpyOfThis, cpyOfOther, nbOfSplits);
747 //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
748 performLocatingOperation(cpyOfOther);
750 for(std::list<ElementaryEdge *>::const_iterator it=cpyOfOther._sub_edges.begin();it!=cpyOfOther._sub_edges.end();it++)
752 switch((*it)->getLoc())
756 ret += fabs((*it)->getPtr()->getCurveLength());
762 ret += fabs((*it)->getPtr()->getCurveLength());
774 * Warning contrary to intersectWith method this method is \b NOT const. 'this' and 'other' are modified after call of this method.
776 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other, double* barycenter)
778 double ret=0.,bary[2],area,xBaryBB,yBaryBB;
779 barycenter[0] = barycenter[1] = 0.;
780 double fact=normalize(&other,xBaryBB,yBaryBB);
781 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
782 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
784 area=fabs((*iter)->getArea());
785 (*iter)->getBarycenter(bary);
788 barycenter[0] += bary[0]*area;
789 barycenter[1] += bary[1]*area;
791 if ( ret > std::numeric_limits<double>::min() )
793 barycenter[0]=barycenter[0]/ret*fact+xBaryBB;
794 barycenter[1]=barycenter[1]/ret*fact+yBaryBB;
797 return ret*fact*fact;
801 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
802 * This is possible because loc attribute in Edge class is mutable.
803 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
805 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other) const
808 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
809 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
811 ret+=fabs((*iter)->getArea());
818 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
819 * This is possible because loc attribute in Edge class is mutable.
820 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
822 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other, double* barycenter) const
824 double ret=0., bary[2];
825 barycenter[0] = barycenter[1] = 0.;
826 std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
827 for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
829 double area = fabs((*iter)->getArea());
830 (*iter)->getBarycenter(bary);
833 barycenter[0] += bary[0]*area;
834 barycenter[1] += bary[1]*area;
836 if ( ret > std::numeric_limits<double>::min() )
838 barycenter[0] /= ret;
839 barycenter[1] /= ret;
845 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
846 * This is possible because loc attribute in Edge class is mutable.
847 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
849 void QuadraticPolygon::intersectForPerimeter(const QuadraticPolygon& other, double& perimeterThisPart, double& perimeterOtherPart, double& perimeterCommonPart) const
851 perimeterThisPart=0.; perimeterOtherPart=0.; perimeterCommonPart=0.;
852 QuadraticPolygon cpyOfThis(*this);
853 QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
854 SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
855 performLocatingOperation(cpyOfOther);
856 other.performLocatingOperation(cpyOfThis);
857 cpyOfThis.dispatchPerimeterExcl(perimeterThisPart,perimeterCommonPart);
858 cpyOfOther.dispatchPerimeterExcl(perimeterOtherPart,perimeterCommonPart);
859 perimeterCommonPart/=2.;
863 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
864 * This is possible because loc attribute in Edge class is mutable.
865 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
867 * polThis.size()==this->size() and polOther.size()==other.size().
868 * For each ElementaryEdge of 'this', the corresponding contribution in resulting polygon is in 'polThis'.
869 * For each ElementaryEdge of 'other', the corresponding contribution in resulting polygon is in 'polOther'.
870 * As consequence common part are counted twice (in polThis \b and in polOther).
872 void QuadraticPolygon::intersectForPerimeterAdvanced(const QuadraticPolygon& other, std::vector< double >& polThis, std::vector< double >& polOther) const
874 polThis.resize(size());
875 polOther.resize(other.size());
876 IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
878 for(it1.first();!it1.finished();it1.next(),edgeId++)
880 ElementaryEdge* curE1=it1.current();
881 QuadraticPolygon cpyOfOther(other);
882 QuadraticPolygon tmp;
883 tmp.pushBack(curE1->clone());
885 SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
886 other.performLocatingOperation(tmp);
887 tmp.dispatchPerimeter(polThis[edgeId]);
890 IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&other));
892 for(it2.first();!it2.finished();it2.next(),edgeId++)
894 ElementaryEdge* curE2=it2.current();
895 QuadraticPolygon cpyOfThis(*this);
896 QuadraticPolygon tmp;
897 tmp.pushBack(curE2->clone());
899 SplitPolygonsEachOther(tmp,cpyOfThis,tmp2);
900 performLocatingOperation(tmp);
901 tmp.dispatchPerimeter(polOther[edgeId]);
907 * numberOfCreatedPointsPerEdge is resized to the number of edges of 'this'.
908 * This method returns in ordered maner the number of newly created points per edge.
909 * This method performs a split process between 'this' and 'other' that gives the result PThis.
910 * Then for each edges of 'this' this method counts how many edges in Pthis have the same id.
912 void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vector< int >& numberOfCreatedPointsPerEdge) const
914 numberOfCreatedPointsPerEdge.resize(size());
915 IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
917 for(it1.first();!it1.finished();it1.next(),edgeId++)
919 ElementaryEdge* curE1=it1.current();
920 QuadraticPolygon cpyOfOther(other);
921 QuadraticPolygon tmp;
922 tmp.pushBack(curE1->clone());
924 SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
925 numberOfCreatedPointsPerEdge[edgeId]=tmp.recursiveSize()-1;
930 * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
931 * This is possible because loc attribute in Edge class is mutable.
932 * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
933 * This method is currently not used by any high level functionality.
935 std::vector<QuadraticPolygon *> QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const
937 QuadraticPolygon cpyOfThis(*this);
938 QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
939 SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
940 //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
941 performLocatingOperation(cpyOfOther);
942 return other.buildIntersectionPolygons(cpyOfOther, cpyOfThis);
946 * This method is typically the first step of boolean operations between pol1 and pol2.
947 * This method perform the minimal splitting so that at the end each edges constituting pol1 are fully either IN or OUT or ON.
948 * @param pol1 IN/OUT param that is equal to 'this' when called.
950 void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits)
952 IteratorOnComposedEdge it1(&pol1),it2(&pol2);
954 ComposedEdge *c1=new ComposedEdge;
955 ComposedEdge *c2=new ComposedEdge;
956 for(it2.first();!it2.finished();it2.next())
958 ElementaryEdge* curE2=it2.current();
959 if(!curE2->isThereStartPoint())
962 it1=curE2->getIterator();
963 for(;!it1.finished();)
966 ElementaryEdge* curE1=it1.current();
967 merge.clear(); nbOfSplits++;
968 if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
970 if(!curE1->getDirection()) c1->reverse();
971 if(!curE2->getDirection()) c2->reverse();
972 UpdateNeighbours(merge,it1,it2,c1,c2);
973 //Substitution of simple edge by sub-edges.
974 delete curE1; // <-- destroying simple edge coming from pol1
975 delete curE2; // <-- destroying simple edge coming from pol2
976 it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
977 it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
980 it1.assignMySelfToAllElems(c2);//To avoid that others
988 UpdateNeighbours(merge,it1,it2,curE1,curE2);
997 void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol1) const
999 IteratorOnComposedEdge it(&pol1);
1000 TypeOfEdgeLocInPolygon loc=FULL_ON_1;
1001 for(it.first();!it.finished();it.next())
1003 ElementaryEdge *cur=it.current();
1004 loc=cur->locateFullyMySelf(*this,loc);//*this=pol2=other
1008 void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) const
1010 IteratorOnComposedEdge it(&pol2);
1011 for(it.first();!it.finished();it.next())
1013 ElementaryEdge *cur=it.current();
1014 cur->locateFullyMySelfAbsolute(*this);
1019 * Given 2 polygons \a pol1 and \a pol2 (localized) the resulting polygons are returned.
1021 * this : pol1 simplified.
1022 * @param [in] pol1 pol1 split.
1023 * @param [in] pol2 pol2 split.
1025 std::vector<QuadraticPolygon *> QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const
1027 std::vector<QuadraticPolygon *> ret;
1028 // Extract from pol1, and pol1 only, all consecutive edges.
1029 // pol1Zip contains concatenated pieces of pol1 which are part of the resulting intersecting cell being built.
1030 std::list<QuadraticPolygon *> pol1Zip=pol1.zipConsecutiveInSegments();
1031 if(!pol1Zip.empty())
1032 ClosePolygons(pol1Zip,*this,pol2,ret);
1034 {//borders of pol1 do not cross pol2,and pol1 borders are outside of pol2. That is to say, either pol1 and pol2
1035 //do not overlap or pol2 is fully inside pol1. So in the first case no intersection, in the other case
1036 //the intersection is pol2.
1037 ElementaryEdge *e1FromPol2=pol2[0];
1038 TypeOfEdgeLocInPolygon loc=FULL_ON_1;
1039 loc=e1FromPol2->locateFullyMySelf(*this,loc);
1041 ret.push_back(new QuadraticPolygon(pol2));
1047 * Returns parts of potentially non closed-polygons. Each returned polygons are not mergeable.
1048 * this : pol1 split and localized.
1050 std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
1052 std::list<QuadraticPolygon *> ret;
1053 IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(this));
1054 int nbOfTurns=recursiveSize();
1056 if(!it.goToNextInOn(false,i,nbOfTurns))
1062 QuadraticPolygon *tmp1=new QuadraticPolygon;
1063 TypeOfEdgeLocInPolygon loc=it.current()->getLoc();
1064 while(loc!=FULL_OUT_1 && i<nbOfTurns)
1066 ElementaryEdge *tmp3=it.current()->clone();
1067 tmp1->pushBack(tmp3);
1069 loc=it.current()->getLoc();
1076 ret.push_back(tmp1);
1077 it.goToNextInOn(true,i,nbOfTurns);
1083 * @param [in] pol1zip is a list of set of edges (=an opened polygon) coming from split polygon 1.
1084 * @param [in] pol1 should be considered as pol1Simplified.
1085 * @param [in] pol2 is split pol2.
1086 * @param [out] results the resulting \b CLOSED polygons.
1088 void QuadraticPolygon::ClosePolygons(std::list<QuadraticPolygon *>& pol1Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2,
1089 std::vector<QuadraticPolygon *>& results)
1091 bool directionKnownInPol2=false;
1092 bool directionInPol2;
1093 bool needCleaning = false;
1094 for(std::list<QuadraticPolygon *>::iterator iter=pol1Zip.begin();iter!=pol1Zip.end();)
1096 // Build incrementally the full closed cells from the consecutive line parts already built in pol1Zip.
1097 // At the end of the process the item currently iterated has been totally completed (start_node=end_node)
1098 // This process can produce several cells.
1099 if((*iter)->completed())
1102 (*iter)->cleanDegeneratedConsecutiveEdges();
1103 results.push_back(*iter);
1104 directionKnownInPol2=false;
1106 iter=pol1Zip.erase(iter);
1109 if(!directionKnownInPol2)
1111 if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol2, needCleaning))
1112 { delete *iter; iter=pol1Zip.erase(iter); continue; }
1114 directionKnownInPol2=true;
1116 std::list<QuadraticPolygon *>::iterator iter2=iter; iter2++;
1117 // Fill as much as possible the current iterate (=a part of pol1) with consecutive pieces from pol2:
1118 std::list<QuadraticPolygon *>::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol2,iter2,pol1Zip.end(),directionInPol2);
1119 // and now add a full connected piece from pol1Zip:
1120 if(iter3!=pol1Zip.end())
1122 (*iter)->pushBack(*iter3);
1124 pol1Zip.erase(iter3);
1130 * 'this' is expected to be set of edges (not closed) of pol1 split.
1132 bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted,
1133 bool& direction, bool& needCleaning) const
1135 needCleaning = false;
1136 IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&pol2Splitted));
1138 Node *n=getEndNode();
1139 ElementaryEdge *cur=it2.current();
1140 // Find edge in pol2 whose start node is the end node of the current piece in pol1Zip (*this)
1141 for(it2.first();!it2.finished() && !found;)
1144 found=(cur->getStartNode()==n);
1149 throw Exception("Internal error: polygons incompatible with each others. Should never happen!");
1150 //Ok we found correspondence between this and pol2. Searching for right direction to close polygon.
1151 ElementaryEdge *e=_sub_edges.back();
1152 if(e->getLoc()==FULL_ON_1)
1154 if(e->getPtr()==cur->getPtr())
1156 // if we have the same edge, several possibilities:
1157 // - either this means that pol1 and pol2 have opposite orientation (since we matched end node with start node before)
1158 // - or (more tricky, see testIntersect2DMeshes11()) that pol1 and pol2 have same orientation but 'this' turns in such
1159 // a way that it attaches to pol2 on an edge in opposite orientation.
1160 // To sort this out, inspect localisation of next edge in pol2 wrt pol1NotSplitted.
1163 Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
1164 bool ret=pol1NotSplitted.isInOrOut(repr);
1167 needCleaning = ret; // if true we are in tricky case 2 above, we know that we will produce two consecutive overlapping edges in result
1170 else // here we don't need to go prev or next:
1172 Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
1173 bool ret=pol1NotSplitted.isInOrOut(repr);
1180 direction=cur->locateFullyMySelfAbsolute(pol1NotSplitted)==FULL_IN_1;
1185 * This method fills as much as possible \a this (a sub-part of pol1 split) with edges of \a pol2Splitted.
1187 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol2Splitted,
1188 std::list<QuadraticPolygon *>::iterator iStart,
1189 std::list<QuadraticPolygon *>::iterator iEnd,
1192 IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(&pol2Splitted));
1194 Node *n=getEndNode();
1195 ElementaryEdge *cur1;
1196 for(it1.first();!it1.finished() && !found;)
1199 found=(cur1->getStartNode()==n);
1206 int szMax(pol2Splitted.size()+1),ii(0); // protection against aggressive users of IntersectMeshes using invalid input meshes ...
1207 std::list<QuadraticPolygon *>::iterator ret;
1209 { // Stack (consecutive) edges of pol1 into the result (no need to care about ordering, edges from pol1 are already consecutive)
1211 ElementaryEdge *tmp=cur1->clone();
1215 nodeToTest=tmp->getEndNode();
1216 direction?it1.nextLoop():it1.previousLoop();
1217 ret=CheckInList(nodeToTest,iStart,iEnd);
1222 while(ret==iEnd && ii<szMax);
1223 if(ii==szMax)// here a protection against aggressive users of IntersectMeshes of invalid input meshes
1224 throw INTERP_KERNEL::Exception("QuadraticPolygon::fillAsMuchAsPossibleWith : Something is invalid with input polygons !");
1228 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
1229 std::list<QuadraticPolygon *>::iterator iEnd)
1231 for(std::list<QuadraticPolygon *>::iterator iter=iStart;iter!=iEnd;iter++)
1232 if((*iter)->isNodeIn(n))
1238 * Compute the remaining parts of the intersection of mesh1 by mesh2.
1239 * The general algorithm is to :
1240 * - either return full cells from pol1 that were simply not touched by mesh2
1242 * - concatenate pieces from pol1 into consecutive pieces (a bit like zipConsecutiveSegments())
1243 * - complete those pieces by edges found in edgesInPol2OnBoundary, which are edges from pol2 located on the boundary of the previously built
1244 * intersecting cells
1246 void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::set<Edge *>& notUsedInPol1, const std::set<Edge *>& edgesInPol2OnBoundary,
1247 const std::map<INTERP_KERNEL::Node *,mcIdType>& mapp, mcIdType offset, mcIdType idThis,
1248 std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& conn,
1249 std::vector<mcIdType>& connI, std::vector<mcIdType>& nb1, std::vector<mcIdType>& nb2)
1251 // Initialise locations on pol1. Remember that edges found in 'notUsedInPol1' are also part of the edges forming pol1.
1252 pol1.initLocations();
1253 for(std::set<Edge *>::const_iterator it1=notUsedInPol1.begin();it1!=notUsedInPol1.end();it1++)
1254 { (*it1)->initLocs(); (*it1)->declareOn(); }
1255 for(std::set<Edge *>::const_iterator it2=edgesInPol2OnBoundary.begin();it2!=edgesInPol2OnBoundary.end();it2++)
1256 { (*it2)->initLocs(); (*it2)->declareIn(); }
1258 std::set<Edge *> notUsedInPol1L(notUsedInPol1);
1259 IteratorOnComposedEdge itPol1(const_cast<QuadraticPolygon *>(&pol1));
1261 std::list<QuadraticPolygon *> pol1Zip;
1262 // If none of the edges of pol1 was consumed by the rebuilding process, we can directly take pol1 as it is to form a cell:
1263 if(pol1.size()==(int)notUsedInPol1.size() && edgesInPol2OnBoundary.empty())
1265 pol1.appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
1268 // Zip consecutive edges found in notUsedInPol1L and which are not overlapping with boundary edge from edgesInPol2OnBoundary:
1269 while(!notUsedInPol1L.empty())
1271 // If all nodes are IN or edge is ON (i.e. as initialised at the begining of the method) then error
1272 for(int i=0;i<sz && (itPol1.current()->getStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1);i++)
1274 if(itPol1.current()->getStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1)
1275 throw INTERP_KERNEL::Exception("Presence of a target polygon fully included in source polygon ! The partition of this leads to a non simply connex cell (with hole) ! Impossible ! Such resulting cell cannot be stored in MED cell format !");
1276 QuadraticPolygon *tmp1=new QuadraticPolygon;
1279 Edge *ee=itPol1.current()->getPtr();
1280 if(ee->getLoc()==FULL_ON_1)
1282 ee->incrRef(); notUsedInPol1L.erase(ee);
1283 tmp1->pushBack(new ElementaryEdge(ee,itPol1.current()->getDirection()));
1287 while(itPol1.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty());
1288 pol1Zip.push_back(tmp1);
1292 std::list<QuadraticPolygon *> retPolsUnderContruction;
1293 std::list<Edge *> edgesInPol2OnBoundaryL(edgesInPol2OnBoundary.begin(),edgesInPol2OnBoundary.end());
1294 std::map<QuadraticPolygon *, std::list<QuadraticPolygon *> > pol1ZipConsumed; // for memory management only.
1295 std::size_t maxNbOfTurn=edgesInPol2OnBoundaryL.size(),nbOfTurn=0,iiMNT=0;
1296 for(std::list<QuadraticPolygon *>::const_iterator itMNT=pol1Zip.begin();itMNT!=pol1Zip.end();itMNT++,iiMNT++)
1297 nbOfTurn+=(*itMNT)->size();
1298 maxNbOfTurn=maxNbOfTurn*nbOfTurn; maxNbOfTurn*=maxNbOfTurn;
1299 // [ABN] at least 3 turns for very small cases (typically one (quad) edge against one (quad or lin) edge forming a new cell)!
1300 maxNbOfTurn = maxNbOfTurn<3 ? 3 : maxNbOfTurn;
1302 while(nbOfTurn<maxNbOfTurn) // the 'normal' way out of this loop is the break towards the end when pol1Zip is empty.
1304 // retPolsUnderConstruction initially empty -> see if(!pol1Zip.empty()) below ...
1305 for(std::list<QuadraticPolygon *>::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();)
1307 if((*itConstr)->getStartNode()==(*itConstr)->getEndNode()) // reconstruction of a cell is finished
1312 Node *curN=(*itConstr)->getEndNode();
1313 bool smthHappened=false;
1314 // Complete a partially reconstructed polygon with boundary edges by matching nodes:
1315 for(std::list<Edge *>::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();)
1317 if(curN==(*it2)->getStartNode())
1318 { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
1319 else if(curN==(*it2)->getEndNode())
1320 { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
1326 for(std::list<QuadraticPolygon *>::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();)
1328 if(curN==(*itZip)->getStartNode()) // we found a matching piece to append in pol1Zip. Append all of it to the current polygon being reconstr
1330 for(std::list<ElementaryEdge *>::const_iterator it4=(*itZip)->_sub_edges.begin();it4!=(*itZip)->_sub_edges.end();it4++)
1331 { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*itConstr)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
1332 pol1ZipConsumed[*itConstr].push_back(*itZip);
1333 curN=(*itZip)->getEndNode();
1334 itZip=pol1Zip.erase(itZip); // one zipped piece has been consumed
1335 break; // we can stop here, pieces in pol1Zip are not connected, by definition.
1343 for(std::list<ElementaryEdge *>::const_iterator it5=(*itConstr)->_sub_edges.begin();it5!=(*itConstr)->_sub_edges.end();it5++)
1345 Edge *ee=(*it5)->getPtr();
1346 if(edgesInPol2OnBoundary.find(ee)!=edgesInPol2OnBoundary.end())
1347 edgesInPol2OnBoundaryL.push_back(ee);
1349 for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++)
1350 pol1Zip.push_front(*it6);
1351 pol1ZipConsumed.erase(*itConstr);
1353 itConstr=retPolsUnderContruction.erase(itConstr);
1356 if(!pol1Zip.empty()) // the filling process of retPolsUnderConstruction starts here
1358 QuadraticPolygon *tmp=new QuadraticPolygon;
1359 QuadraticPolygon *first=*(pol1Zip.begin());
1360 for(std::list<ElementaryEdge *>::const_iterator it4=first->_sub_edges.begin();it4!=first->_sub_edges.end();it4++)
1361 { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); tmp->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
1362 pol1ZipConsumed[tmp].push_back(first);
1363 retPolsUnderContruction.push_back(tmp);
1364 pol1Zip.erase(pol1Zip.begin());
1370 if(nbOfTurn==maxNbOfTurn)
1372 std::ostringstream oss; oss << "Error during reconstruction of residual of cell ! It appears that either source or/and target mesh is/are not conform !";
1373 oss << " Number of turns is = " << nbOfTurn << " !";
1374 throw INTERP_KERNEL::Exception(oss.str().c_str());
1376 // Convert to integer connectivity:
1377 for(std::list<QuadraticPolygon *>::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();itConstr++)
1379 if((*itConstr)->getStartNode()==(*itConstr)->getEndNode()) // take only fully closed reconstructed polygon
1381 (*itConstr)->cleanDegeneratedConsecutiveEdges();
1382 (*itConstr)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
1383 for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++)
1389 std::ostringstream oss; oss << "Internal error during reconstruction of residual of cell! Non fully closed polygon built!";
1390 throw INTERP_KERNEL::Exception(oss.str().c_str());