X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FINTERP_KERNEL%2FGeometric2D%2FInterpKernelGeo2DQuadraticPolygon.cxx;h=f6c1ee3164a5cbcf525d754ca07fad61c474a9b4;hb=aafcf704892f03308a84407e898d9e8b19496a1c;hp=9648a3ddc5cef5e27fb8f5bac2c5809fe7e77fd0;hpb=ed26f73bb07b412b350ed1e42ec92bbb4f08a0e0;p=tools%2Fmedcoupling.git diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 9648a3ddc..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; } } @@ -1193,48 +1234,64 @@ std::list::iterator QuadraticPolygon::CheckInList(Node *n, s return iEnd; } -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) +/*! +* Compute the remaining parts of the intersection of mesh1 by mesh2. +* The general algorithm is to : +* - either return full cells from pol1 that were simply not touched by mesh2 +* - or to: +* - concatenate pieces from pol1 into consecutive pieces (a bit like zipConsecutiveSegments()) +* - complete those pieces by edges found in edgesInPol2OnBoundary, which are edges from pol2 located on the boundary of the previously built +* intersecting cells +*/ +void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::set& notUsedInPol1, const std::set& edgesInPol2OnBoundary, + 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(); - for(std::set::const_iterator it9=notUsedInPol1.begin();it9!=notUsedInPol1.end();it9++) - { (*it9)->initLocs(); (*it9)->declareOn(); } - for(std::set::const_iterator itA=edgesInPol2OnBoundary.begin();itA!=edgesInPol2OnBoundary.end();itA++) - { (*itA)->initLocs(); (*itA)->declareIn(); } + for(std::set::const_iterator it1=notUsedInPol1.begin();it1!=notUsedInPol1.end();it1++) + { (*it1)->initLocs(); (*it1)->declareOn(); } + for(std::set::const_iterator it2=edgesInPol2OnBoundary.begin();it2!=edgesInPol2OnBoundary.end();it2++) + { (*it2)->initLocs(); (*it2)->declareIn(); } //// std::set notUsedInPol1L(notUsedInPol1); - IteratorOnComposedEdge it(const_cast(&pol1)); + IteratorOnComposedEdge itPol1(const_cast(&pol1)); int sz=pol1.size(); std::list pol1Zip; + // 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: if(pol1.size()==(int)notUsedInPol1.size() && edgesInPol2OnBoundary.empty()) { pol1.appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1); return ; } + // Zip consecutive edges found in notUsedInPol1L and which are not overlapping with boundary edge from edgesInPol2OnBoundary: while(!notUsedInPol1L.empty()) { - for(int i=0;igetStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1);i++) - it.nextLoop(); - if(it.current()->getStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1) + // If all nodes are IN or edge is ON (i.e. as initialised at the begining of the method) then error + for(int i=0;igetStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1);i++) + itPol1.nextLoop(); + if(itPol1.current()->getStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1) 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 !"); QuadraticPolygon *tmp1=new QuadraticPolygon; do { - Edge *ee=it.current()->getPtr(); + Edge *ee=itPol1.current()->getPtr(); if(ee->getLoc()==FULL_ON_1) { ee->incrRef(); notUsedInPol1L.erase(ee); - tmp1->pushBack(new ElementaryEdge(ee,it.current()->getDirection())); + tmp1->pushBack(new ElementaryEdge(ee,itPol1.current()->getDirection())); } - it.nextLoop(); + itPol1.nextLoop(); } - while(it.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty()); + while(itPol1.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty()); pol1Zip.push_back(tmp1); } + //// std::list retPolsUnderContruction; std::list edgesInPol2OnBoundaryL(edgesInPol2OnBoundary.begin(),edgesInPol2OnBoundary.end()); - std::map > pol1ZipConsumed; + std::map > pol1ZipConsumed; // for memory management only. std::size_t maxNbOfTurn=edgesInPol2OnBoundaryL.size(),nbOfTurn=0,iiMNT=0; for(std::list::const_iterator itMNT=pol1Zip.begin();itMNT!=pol1Zip.end();itMNT++,iiMNT++) nbOfTurn+=(*itMNT)->size(); @@ -1242,59 +1299,92 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: // [ABN] at least 3 turns for very small cases (typically one (quad) edge against one (quad or lin) edge forming a new cell)! maxNbOfTurn = maxNbOfTurn<3 ? 3 : maxNbOfTurn; nbOfTurn=0; - while(nbOfTurn::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();) + // retPolsUnderConstruction initially empty -> see if(!pol1Zip.empty()) below ... + for(std::list::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();) { - if((*it1)->getStartNode()==(*it1)->getEndNode()) - { - it1++; - continue; - } - Node *curN=(*it1)->getEndNode(); - bool smthHappened=false; + 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(); (*it1)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); } - else if(curN==(*it2)->getEndNode()) - { (*it2)->incrRef(); (*it1)->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; } } - if(smthHappened) + + // 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 it3=pol1Zip.begin();it3!=pol1Zip.end();) + for(std::list::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();) { - if(curN==(*it3)->getStartNode()) + if(curN==(*it2)->getStartNode()) { - for(std::list::const_iterator it4=(*it3)->_sub_edges.begin();it4!=(*it3)->_sub_edges.end();it4++) - { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*it1)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); } + (*it2)->incrRef(); + (*itConstr)->pushBack(new ElementaryEdge(*it2,true)); + curN=(*it2)->getEndNode(); smthHappened=true; - pol1ZipConsumed[*it1].push_back(*it3); - curN=(*it3)->getEndNode(); - it3=pol1Zip.erase(it3); + it2=edgesInPol2OnBoundaryL.erase(it2); } else - it3++; + 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::const_iterator it5=(*it1)->_sub_edges.begin();it5!=(*it1)->_sub_edges.end();it5++) + for(std::list::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();) + { + if(curN==(*itZip)->getStartNode()) // we found a matching piece to append in pol1Zip. Append all of it to the current polygon being reconstr + { + for(std::list::const_iterator it4=(*itZip)->_sub_edges.begin();it4!=(*itZip)->_sub_edges.end();it4++) + { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*itConstr)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); } + pol1ZipConsumed[*itConstr].push_back(*itZip); + curN=(*itZip)->getEndNode(); + itZip=pol1Zip.erase(itZip); // one zipped piece has been consumed + break; // we can stop here, pieces in pol1Zip are not connected, by definition. + } + else + itZip++; + } + } + else // Nothing happened. + { + for(std::list::const_iterator it5=(*itConstr)->_sub_edges.begin();it5!=(*itConstr)->_sub_edges.end();it5++) { Edge *ee=(*it5)->getPtr(); if(edgesInPol2OnBoundary.find(ee)!=edgesInPol2OnBoundary.end()) edgesInPol2OnBoundaryL.push_back(ee); } - for(std::list::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++) + for(std::list::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++) pol1Zip.push_front(*it6); - pol1ZipConsumed.erase(*it1); - delete *it1; - it1=retPolsUnderContruction.erase(it1); + pol1ZipConsumed.erase(*itConstr); + delete *itConstr; + itConstr=retPolsUnderContruction.erase(itConstr); } } - if(!pol1Zip.empty()) + if(!pol1Zip.empty()) // the filling process of retPolsUnderConstruction starts here { QuadraticPolygon *tmp=new QuadraticPolygon; QuadraticPolygon *first=*(pol1Zip.begin()); @@ -1304,6 +1394,8 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: retPolsUnderContruction.push_back(tmp); pol1Zip.erase(pol1Zip.begin()); } + else + break; nbOfTurn++; } if(nbOfTurn==maxNbOfTurn) @@ -1312,14 +1404,21 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: oss << " Number of turns is = " << nbOfTurn << " !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } - for(std::list::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();it1++) + // Convert to integer connectivity: + for(std::list::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();itConstr++) { - if((*it1)->getStartNode()==(*it1)->getEndNode()) + if((*itConstr)->getStartNode()==(*itConstr)->getEndNode()) // take only fully closed reconstructed polygon { - (*it1)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1); - for(std::list::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++) + (*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; - delete *it1; + delete *itConstr; + } + else + { + std::ostringstream oss; oss << "Internal error during reconstruction of residual of cell! Non fully closed polygon built!"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); } } }