X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FGeometric2D%2FInterpKernelGeo2DQuadraticPolygon.cxx;h=f6c1ee3164a5cbcf525d754ca07fad61c474a9b4;hb=aafcf704892f03308a84407e898d9e8b19496a1c;hp=26c7a72231b15a7ca231edaa76fad324fed389da;hpb=66423193c07c9c8f10816594d23edd7271c7a2a1;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 26c7a7223..f6c1ee316 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D +// Copyright (C) 2007-2022 CEA/DEN, EDF R&D // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -284,12 +284,12 @@ double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other) * the cell id in global other mesh. */ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, - const std::map& mapThis, const std::map& mapOther, - int offset1, int offset2 , - const std::vector& otherEdgeIds, - std::vector& edgesThis, int cellIdThis, - std::vector< std::vector >& edgesInOtherColinearWithThis, std::vector< std::vector >& subDivOther, - std::vector& addCoo, std::map& mergedNodes) + const std::map& mapThis, const std::map& mapOther, + mcIdType offset1, mcIdType offset2 , + const std::vector& otherEdgeIds, + std::vector& edgesThis, mcIdType cellIdThis, + std::vector< std::vector >& edgesInOtherColinearWithThis, std::vector< std::vector >& subDivOther, + std::vector& addCoo, std::map& mergedNodes) { double xBaryBB, yBaryBB; double fact=normalizeExt(&other, xBaryBB, yBaryBB); @@ -300,7 +300,7 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, ComposedEdge *cThis=new ComposedEdge; ComposedEdge *cOther=new ComposedEdge; int i=0; - std::map mapAddCoo; + std::map mapAddCoo; for(itOther.first();!itOther.finished();itOther.next(),i++) { // For each edge of 'other', proceed with intersections: the edge might split into sub-edges, 'otherTmp' will hold the final split result. @@ -322,9 +322,9 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, ElementaryEdge* curThis=itThis.current(); merge.clear(); // - std::map::const_iterator thisStart(mapThis.find(curThis->getStartNode())),thisEnd(mapThis.find(curThis->getEndNode())), + std::map::const_iterator thisStart(mapThis.find(curThis->getStartNode())),thisEnd(mapThis.find(curThis->getEndNode())), otherStart(mapOther.find(curOtherTmp->getStartNode())),otherEnd(mapOther.find(curOtherTmp->getEndNode())); - int thisStart2(thisStart==mapThis.end()?-1:(*thisStart).second), thisEnd2(thisEnd==mapThis.end()?-1:(*thisEnd).second), + mcIdType thisStart2(thisStart==mapThis.end()?-1:(*thisStart).second), thisEnd2(thisEnd==mapThis.end()?-1:(*thisEnd).second), otherStart2(otherStart==mapOther.end()?-1:(*otherStart).second+offset1),otherEnd2(otherEnd==mapOther.end()?-1:(*otherEnd).second+offset1); // if(curThis->getPtr()->intersectWith(curOtherTmp->getPtr(),merge,*cThis,*cOther)) @@ -362,8 +362,14 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, // Converting back to integer connectivity: if(otherTmp._sub_edges.size()>1) // only if a new point has been added (i.e. an actual intersection was done) { - for(std::list::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++) - (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo); + std::size_t jj = 0, sz(otherTmp._sub_edges.size()); + for(std::list::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++, jj++) + { + short skipStartOrEnd = jj == 0 ? -1 : (jj == sz-1 ? 1 : 0); // -1 means START, 1 means END, 0 other + (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2, + fact,xBaryBB,yBaryBB, skipStartOrEnd, + /*out*/ subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo); + } } } Delete(cThis); @@ -380,8 +386,8 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, * orientation of edge (see buildDescendingConnectivity2() method). * See appendEdgeFromCrudeDataArray() for params description. */ -void QuadraticPolygon::buildFromCrudeDataArray(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, - const int *descBg, const int *descEnd, const std::vector >& intersectEdges) +void QuadraticPolygon::buildFromCrudeDataArray(const std::map& mapp, bool isQuad, const mcIdType *nodalBg, const double *coords, + const mcIdType *descBg, const mcIdType *descEnd, const std::vector >& intersectEdges) { std::size_t nbOfSeg=std::distance(descBg,descEnd); for(std::size_t i=0;i& mapp, bool isQuad, - const int *nodalBg, const double *coords, - const int *descBg, const int *descEnd, const std::vector >& intersectEdges) +void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map& mapp, bool isQuad, + const mcIdType *nodalBg, const double *coords, + const mcIdType *descBg, const mcIdType *descEnd, const std::vector >& intersectEdges) { if(!isQuad) { bool direct=descBg[edgePos]>0; - int edgeId=abs(descBg[edgePos])-1; // back to C indexing mode - const std::vector& subEdge=intersectEdges[edgeId]; + mcIdType edgeId=std::abs(descBg[edgePos])-1; // back to C indexing mode + const std::vector& subEdge=intersectEdges[edgeId]; std::size_t nbOfSubEdges=subEdge.size()/2; for(std::size_t j=0;j0; - int edgeId=abs(descBg[edgePos])-1; - const std::vector& subEdge=intersectEdges[edgeId]; + mcIdType edgeId=std::abs(descBg[edgePos])-1; + const std::vector& subEdge=intersectEdges[edgeId]; std::size_t nbOfSubEdges=subEdge.size()/2; if(colinearity) { @@ -439,7 +445,7 @@ void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const s } } -void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, int edgeId, const std::vector& subEdge, const std::map& mapp) +void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, mcIdType edgeId, const std::vector& subEdge, const std::map& mapp) { std::size_t nbOfSubEdges=subEdge.size()/2; if(!baseEdge) @@ -463,17 +469,17 @@ void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size * 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 * orientation of edge. */ -void QuadraticPolygon::buildFromCrudeDataArray2(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector >& intersectEdges, - const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector >& intersectEdges1, - const std::vector< std::vector >& colinear1, - std::map >& alreadyExistingIn2) +void QuadraticPolygon::buildFromCrudeDataArray2(const std::map& mapp, bool isQuad, const mcIdType *nodalBg, const double *coords, const mcIdType *descBg, const mcIdType *descEnd, const std::vector >& intersectEdges2, + const INTERP_KERNEL::QuadraticPolygon& pol1, const mcIdType *descBg1, const mcIdType *descEnd1, const std::vector >& intersectEdges1, + const std::vector< std::vector >& colinear1, + std::map >& alreadyExistingIn2) { std::size_t nbOfSeg=std::distance(descBg,descEnd); for(std::size_t i=0;i0; - int edgeId=abs(descBg[i])-1;//current edge id of pol2 - std::map >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]); + mcIdType edgeId=std::abs(descBg[i])-1;//current edge id of pol2 + std::map >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]); if(it1!=alreadyExistingIn2.end() || it2!=alreadyExistingIn2.end()) { bool sameDir=(it1!=alreadyExistingIn2.end()); @@ -497,28 +503,28 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map > > idIns1; - int offset1=0; + std::vector > > idIns1; + mcIdType offset1=0; if(!directos) {// if the current edge of pol2 has one or more colinear edges part into pol1 - const std::vector& c=colinear1[edgeId]; + const std::vector& c=colinear1[edgeId]; std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1); for(std::size_t j=0;j >(edgeId1,std::pair(descBg1[j]>0,offset1)));// it exists an edge into pol1 given by tuple (idIn1,direct1) that is colinear at edge 'edgeId' in pol2 + idIns1.push_back(std::pair >(edgeId1,std::pair(descBg1[j]>0,offset1)));// it exists an edge into pol1 given by tuple (idIn1,direct1) that is colinear at edge 'edgeId' in pol2 //std::pair0; } - offset1+=intersectEdges1[edgeId1].size()/2;//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2 + 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 } directos=idIns1.empty(); } if(directos) {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method std::size_t oldSz=_sub_edges.size(); - appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges); + appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges2); std::size_t newSz=_sub_edges.size(); std::size_t zeSz=newSz-oldSz; alreadyExistingIn2[descBg[i]].resize(zeSz); @@ -528,25 +534,25 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map& subEdge=intersectEdges[edgeId]; + const std::vector& subEdge=intersectEdges2[edgeId]; std::size_t nbOfSubEdges=subEdge.size()/2; for(std::size_t j=0;j > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++) + mcIdType idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1]; + mcIdType idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2]; + bool direction11=false,found=false; + bool direct1=false;//store if needed the direction in 1 + mcIdType offset2=0; + mcIdType nbOfSubEdges1=0; + for(std::vector > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++) { - int idIn1=(*it).first;//store if needed the cell id in 1 + mcIdType idIn1=(*it).first;//store if needed the cell id in 1 direct1=(*it).second.first; offset1=(*it).second.second; - const std::vector& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1]; - nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2; + const std::vector& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1]; + nbOfSubEdges1=ToIdType(subEdge1PossiblyAlreadyIn1.size()/2); offset2=0; - for(std::size_t k=0;k reuse Edge instance of pol1 - ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)]; + ElementaryEdge *e=pol1[FromIdType(offset1+(direct1?offset2:nbOfSubEdges1-offset2-1))]; Edge *ee=e->getPtr(); ee->incrRef(); ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11)); @@ -583,39 +589,39 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map >& intersectEdges, - const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, - const std::vector >& intersectEdges1, const std::vector< std::vector >& colinear1) const +void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const mcIdType *descBg, const mcIdType *descEnd, const std::vector >& intersectEdges, + const INTERP_KERNEL::QuadraticPolygon& pol1, const mcIdType *descBg1, const mcIdType *descEnd1, + const std::vector >& intersectEdges1, const std::vector< std::vector >& colinear1) const { std::size_t nbOfSeg=std::distance(descBg,descEnd); for(std::size_t i=0;i0; - int edgeId=abs(descBg[i])-1;//current edge id of pol2 - const std::vector& c=colinear1[edgeId]; + mcIdType edgeId=std::abs(descBg[i])-1;//current edge id of pol2 + const std::vector& c=colinear1[edgeId]; if(c.empty()) continue; - const std::vector& subEdge=intersectEdges[edgeId]; + const std::vector& subEdge=intersectEdges[edgeId]; std::size_t nbOfSubEdges=subEdge.size()/2; // std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1); - int offset1=0; + mcIdType offset1=0; for(std::size_t j=0;j0; - const std::vector& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1]; - std::size_t nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2; - int offset2=0; + const std::vector& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1]; + mcIdType nbOfSubEdges1=ToIdType(subEdge1PossiblyAlreadyIn1.size()/2); + mcIdType offset2=0; bool found=false; - for(std::size_t kk=0;kkgetPtr()->declareOn(); } } } - offset1+=intersectEdges1[edgeId1].size()/2;//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2 + 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 } } } -void QuadraticPolygon::appendCrudeData(const std::map& mapp, double xBary, double yBary, double fact, int offset, std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI) const +void QuadraticPolygon::appendCrudeData(const std::map& mapp, double xBary, double yBary, double fact, mcIdType offset, std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI) const { int nbOfNodesInPg=0; bool presenceOfQuadratic=presenceOfQuadraticEdge(); @@ -642,14 +648,14 @@ void QuadraticPolygon::appendCrudeData(const std::map { Node *tmp=0; tmp=(*it)->getStartNode(); - std::map::const_iterator it1=mapp.find(tmp); + std::map::const_iterator it1=mapp.find(tmp); conn.push_back((*it1).second); nbOfNodesInPg++; } if(presenceOfQuadratic) { int j=0; - int off=offset+((int)addCoordsQuadratic.size())/2; + mcIdType off=offset+ToIdType(addCoordsQuadratic.size())/2; for(std::list::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++) { INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf(); @@ -666,13 +672,13 @@ void QuadraticPolygon::appendCrudeData(const std::map /*! * 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. * 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. - * @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 - * @param [in,out] edgesBoundaryOther, parameter that stores all edges in result of intersection that are not + * @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 + * @param [in,out] edgesBoundaryOther parameter that stores all edges in result of intersection that are not */ void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set& edgesThis, std::set& edgesBoundaryOther, - const std::map& mapp, int idThis, int idOther, int offset, - std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI, - std::vector& nbThis, std::vector& nbOther) + const std::map& mapp, mcIdType idThis, mcIdType idOther, mcIdType offset, + std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI, + std::vector& nbThis, std::vector& nbOther) { double xBaryBB, yBaryBB; double fact=normalizeExt(&other, xBaryBB, yBaryBB); @@ -703,6 +709,28 @@ void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set 2) + for(it.first();!it.finished();it.next()) + { + ElementaryEdge * cur = it.current(); + if (prevEdge && prevEdge->hasSameExtremities(*cur)) + { + it.eraseCurrent(); + it.eraseCurrent(); + prevEdge = it.current(); + } + else + prevEdge = cur; + } +} + /*! * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method. * 'other' is a QuadraticPolygon of \b non closed edges. @@ -1052,7 +1080,7 @@ std::list QuadraticPolygon::zipConsecutiveInSegments() const } /*! - * @param [in] pol1zip is a list of set of edges (=an opened polygon) coming from split polygon 1. + * @param [in] pol1Zip is a list of set of edges (=an opened polygon) coming from split polygon 1. * @param [in] pol1 should be considered as pol1Simplified. * @param [in] pol2 is split pol2. * @param [out] results the resulting \b CLOSED polygons. @@ -1061,7 +1089,8 @@ void QuadraticPolygon::ClosePolygons(std::list& pol1Zip, con std::vector& results) { bool directionKnownInPol2=false; - bool directionInPol2; + bool directionInPol2=false; + bool needCleaning = false; for(std::list::iterator iter=pol1Zip.begin();iter!=pol1Zip.end();) { // Build incrementally the full closed cells from the consecutive line parts already built in pol1Zip. @@ -1069,14 +1098,17 @@ void QuadraticPolygon::ClosePolygons(std::list& pol1Zip, con // This process can produce several cells. if((*iter)->completed()) { + if (needCleaning) + (*iter)->cleanDegeneratedConsecutiveEdges(); results.push_back(*iter); directionKnownInPol2=false; + needCleaning=false; iter=pol1Zip.erase(iter); continue; } if(!directionKnownInPol2) { - if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol2)) + if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol2, needCleaning)) { delete *iter; iter=pol1Zip.erase(iter); continue; } else directionKnownInPol2=true; @@ -1097,12 +1129,15 @@ void QuadraticPolygon::ClosePolygons(std::list& pol1Zip, con /*! * 'this' is expected to be set of edges (not closed) of pol1 split. */ -bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted, bool& direction) const +bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted, + bool& direction, bool& needCleaning) const { + needCleaning = false; IteratorOnComposedEdge it2(const_cast(&pol2Splitted)); bool found=false; Node *n=getEndNode(); ElementaryEdge *cur=it2.current(); + // Find edge in pol2 whose start node is the end node of the current piece in pol1Zip (*this) for(it2.first();!it2.finished() && !found;) { cur=it2.current(); @@ -1118,20 +1153,26 @@ bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1N { if(e->getPtr()==cur->getPtr()) { - direction=false; - it2.previousLoop(); + // if we have the same edge, several possibilities: + // - either this means that pol1 and pol2 have opposite orientation (since we matched end node with start node before) + // - or (more tricky, see testIntersect2DMeshes11()) that pol1 and pol2 have same orientation but 'this' turns in such + // a way that it attaches to pol2 on an edge in opposite orientation. + // To sort this out, inspect localisation of next edge in pol2 wrt pol1NotSplitted. + it2.nextLoop(); cur=it2.current(); Node *repr=cur->getPtr()->buildRepresentantOfMySelf(); bool ret=pol1NotSplitted.isInOrOut(repr); repr->decrRef(); + direction = ret; + needCleaning = ret; // if true we are in tricky case 2 above, we know that we will produce two consecutive overlapping edges in result return ret; } - else + else // here we don't need to go prev or next: { - direction=true; Node *repr=cur->getPtr()->buildRepresentantOfMySelf(); bool ret=pol1NotSplitted.isInOrOut(repr); repr->decrRef(); + direction = ret; return ret; } } @@ -1203,9 +1244,9 @@ std::list::iterator QuadraticPolygon::CheckInList(Node *n, s * intersecting cells */ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::set& notUsedInPol1, const std::set& edgesInPol2OnBoundary, - const std::map& mapp, int offset, int idThis, - std::vector& addCoordsQuadratic, std::vector& conn, - std::vector& connI, std::vector& nb1, std::vector& nb2) + const std::map& mapp, mcIdType offset, mcIdType idThis, + std::vector& addCoordsQuadratic, std::vector& conn, + std::vector& connI, std::vector& nb1, std::vector& nb2) { // Initialise locations on pol1. Remember that edges found in 'notUsedInPol1' are also part of the edges forming pol1. pol1.initLocations(); @@ -1263,24 +1304,55 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: // retPolsUnderConstruction initially empty -> see if(!pol1Zip.empty()) below ... for(std::list::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();) { - if((*itConstr)->getStartNode()==(*itConstr)->getEndNode()) // reconstruction of a cell is finished - { - itConstr++; - continue; - } - Node *curN=(*itConstr)->getEndNode(); - bool smthHappened=false; - // Complete a partially reconstructed polygon with boundary edges by matching nodes: + Node *startN = (*itConstr)->getStartNode(); + Node *curN = (*itConstr)->getEndNode(); + if(startN == curN) // reconstruction of a cell is finished + { itConstr++; continue; } + + bool smthHappened=false, doneEarly=false; + // Complete a partially reconstructed polygon with boundary edges of pol2 by matching nodes: for(std::list::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();) { - if(curN==(*it2)->getStartNode()) - { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); } - else if(curN==(*it2)->getEndNode()) - { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); } + if(curN==(*it2)->getEndNode()) // only end node should be considered if orientation is correct for input meshes + // in the funny case of cells exactly included (see test case testIntersect2DMeshesTmp13) this is mandatory to take edges from pol2 in the right order. + { + (*it2)->incrRef(); + (*itConstr)->pushBack(new ElementaryEdge(*it2,false)); + curN=(*it2)->getStartNode(); + smthHappened=true; + it2=edgesInPol2OnBoundaryL.erase(it2); + } else it2++; + if (curN == startN) // we might end here + { doneEarly = true; break; } + } + + // It might be the case that the lookup on start nodes made above failed because pol2 is wrongly oriented. + // Be somewhat flexible and keep on supporting this case here (useful for voronisation notably): + if(!smthHappened) + { + for(std::list::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();) + { + if(curN==(*it2)->getStartNode()) + { + (*it2)->incrRef(); + (*itConstr)->pushBack(new ElementaryEdge(*it2,true)); + curN=(*it2)->getEndNode(); + smthHappened=true; + it2=edgesInPol2OnBoundaryL.erase(it2); + } + else + it2++; + if (curN == startN) // we might end here + { doneEarly = true; break; } + } } - if(smthHappened) + + if (doneEarly) + { itConstr++; continue; } + + if(smthHappened) // Now continue the construction by finding the next bit in pol1Zip. Not too sure what are the cases where the boolean if False ... { for(std::list::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();) { @@ -1297,7 +1369,7 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: itZip++; } } - else + else // Nothing happened. { for(std::list::const_iterator it5=(*itConstr)->_sub_edges.begin();it5!=(*itConstr)->_sub_edges.end();it5++) { @@ -1335,8 +1407,9 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: // Convert to integer connectivity: for(std::list::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();itConstr++) { - if((*itConstr)->getStartNode()==(*itConstr)->getEndNode()) // take only fully closed reconstructed polygon (?? might there be others??) + if((*itConstr)->getStartNode()==(*itConstr)->getEndNode()) // take only fully closed reconstructed polygon { + (*itConstr)->cleanDegeneratedConsecutiveEdges(); (*itConstr)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1); for(std::list::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++) delete *it6;