From 4d3a844eef692e30e399e8f2864f55731e59d26e Mon Sep 17 00:00:00 2001 From: abn Date: Wed, 9 Jan 2019 13:10:47 +0100 Subject: [PATCH] Intersection : ComputeResidual -> some refactoring for readability --- .../InterpKernelGeo2DComposedEdge.cxx | 17 +--- .../Geometric2D/InterpKernelGeo2DEdge.cxx | 16 +++ .../Geometric2D/InterpKernelGeo2DEdge.hxx | 1 + .../Geometric2D/InterpKernelGeo2DEdgeLin.cxx | 8 -- .../InterpKernelGeo2DElementaryEdge.cxx | 5 + .../InterpKernelGeo2DElementaryEdge.hxx | 1 + .../InterpKernelGeo2DQuadraticPolygon.cxx | 99 +++++++++++-------- 7 files changed, 81 insertions(+), 66 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx index 55eda23f6..3d0f0a468 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx @@ -314,25 +314,10 @@ void ComposedEdge::dumpToCout(const std::map& mapp) c { int i=0; for(std::list::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++, i++) - { - const ElementaryEdge& e(*(*iter)); - auto sI(mapp.find(e.getStartNode())), eI(mapp.find(e.getEndNode())); - int start = (sI == mapp.end() ? -1 : sI->second), end = (eI == mapp.end() ? -1 : eI->second); - std::string locs; - switch (e.getLoc()) - { - case FULL_IN_1: locs="FULL_IN_1"; break; - case FULL_ON_1: locs="FULL_ON_1"; break; - case FULL_OUT_1: locs="FULL_OUT_1"; break; - case FULL_UNKNOWN: locs="FULL_UNKNOWN"; break; - default: locs="oh my God! This is so wrong."; - } - std::cout << "Edge [" << i << "] : ("<< std::hex << &e << std::dec << ") -> (" << start << ", " << end << ")\t" << locs << std::endl; - } + (*iter)->dumpToCout(mapp, i); std::cout << std::endl; } - Node *ComposedEdge::getEndNode() const { return _sub_edges.back()->getEndNode(); diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx index 24051cbdb..e5b3c1efc 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx @@ -926,6 +926,22 @@ bool Edge::SplitOverlappedEdges(const Edge *e1, const Edge *e2, Node *nS, Node * } } +void Edge::dumpToCout(const std::map& mapp, int index) const +{ + auto sI(mapp.find(getStartNode())), eI(mapp.find(getEndNode())); + int start = (sI == mapp.end() ? -1 : sI->second), end = (eI == mapp.end() ? -1 : eI->second); + std::string locs; + switch (getLoc()) + { + case FULL_IN_1: locs="FULL_IN_1"; break; + case FULL_ON_1: locs="FULL_ON_1"; break; + case FULL_OUT_1: locs="FULL_OUT_1"; break; + case FULL_UNKNOWN: locs="FULL_UNKNOWN"; break; + default: locs="oh my God! This is so wrong."; + } + std::cout << "Edge [" << index << "] : ("<< std::hex << this << std::dec << ") -> (" << start << ", " << end << ")\t" << locs << std::endl; +} + bool Edge::isEqual(const Edge& other) const { return _start->isEqual(*other._start) && _end->isEqual(*other._end); diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx index e5501ad1a..2ed413efe 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx @@ -276,6 +276,7 @@ namespace INTERP_KERNEL static void Interpolate1DLin(const std::vector& distrib1, const std::vector& distrib2, std::map >& result); virtual void dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const = 0; + void dumpToCout(const std::map& mapp, int index) const; bool isEqual(const Edge& other) const; public: bool sortSubNodesAbs(const double *coo, std::vector& subNodes); diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx index 55eeff967..44b99017f 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx @@ -32,19 +32,11 @@ namespace INTERP_KERNEL SegSegIntersector::SegSegIntersector(const EdgeLin& e1, const EdgeLin& e2):SameTypeEdgeIntersector(e1,e2) { -// _matrix[0]=(*(e2.getStartNode()))[0]-(*(e2.getEndNode()))[0]; -// _matrix[1]=(*(e1.getEndNode()))[0]-(*(e1.getStartNode()))[0]; -// _matrix[2]=(*(e2.getStartNode()))[1]-(*(e2.getEndNode()))[1]; -// _matrix[3]=(*(e1.getEndNode()))[1]-(*(e1.getStartNode()))[1]; - _matrix[0]=(*(e1.getEndNode()))[0]-(*(e1.getStartNode()))[0]; _matrix[1]=(*(e1.getEndNode()))[1]-(*(e1.getStartNode()))[1]; _matrix[2]=(*(e2.getEndNode()))[0]-(*(e2.getStartNode()))[0]; _matrix[3]=(*(e2.getEndNode()))[1]-(*(e2.getStartNode()))[1]; - -// _col[0]=_matrix[3]*(*(e1.getStartNode()))[0]-_matrix[1]*(*(e1.getStartNode()))[1]; -// _col[1]=-_matrix[2]*(*(e2.getStartNode()))[0]+_matrix[0]*(*(e2.getStartNode()))[1]; _col[0]=_matrix[1]*(*(e1.getStartNode()))[0]-_matrix[0]*(*(e1.getStartNode()))[1]; _col[1]=_matrix[3]*(*(e2.getStartNode()))[0]-_matrix[2]*(*(e2.getStartNode()))[1]; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx index 9ea7ac851..0f62a67c0 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx @@ -183,6 +183,11 @@ void ElementaryEdge::dumpInXfigFile(std::ostream& stream, int resolution, const _ptr->dumpInXfigFile(stream,_direction,resolution,box); } +void ElementaryEdge::dumpToCout(const std::map& mapp, int index) const +{ + _ptr->dumpToCout(mapp, index); +} + bool ElementaryEdge::intresicEqual(const ElementaryEdge *other) const { return _ptr==other->_ptr; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx index 340bc63b1..18f6efafd 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx @@ -64,6 +64,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT bool intresicEqual(const ElementaryEdge *other) const; INTERPKERNEL_EXPORT bool intresicEqualDirSensitive(const ElementaryEdge *other) const; INTERPKERNEL_EXPORT void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; + INTERPKERNEL_EXPORT void dumpToCout(const std::map& mapp, int index) const; INTERPKERNEL_EXPORT bool getDirection() const { return _direction; } INTERPKERNEL_EXPORT bool intresincEqCoarse(const Edge *other) const; INTERPKERNEL_EXPORT bool isEqual(const ElementaryEdge& other) const; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 9648a3ddc..74ffdc372 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -1193,48 +1193,55 @@ 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) +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) { + // 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 ON (why ON?) 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(); @@ -1244,57 +1251,59 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: 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()) + if((*itConstr)->getStartNode()==(*itConstr)->getEndNode()) // reconstruction of a cell is finished { - it1++; + itConstr++; continue; } - Node *curN=(*it1)->getEndNode(); + Node *curN=(*itConstr)->getEndNode(); bool smthHappened=false; + // Complete a partially reconstructed polygon with boundary edges 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); } + { (*it2)->incrRef(); (*itConstr)->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); } + { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); } else it2++; } if(smthHappened) { - for(std::list::iterator it3=pol1Zip.begin();it3!=pol1Zip.end();) + for(std::list::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();) { - if(curN==(*it3)->getStartNode()) + 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=(*it3)->_sub_edges.begin();it4!=(*it3)->_sub_edges.end();it4++) - { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*it1)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); } - smthHappened=true; - pol1ZipConsumed[*it1].push_back(*it3); - curN=(*it3)->getEndNode(); - it3=pol1Zip.erase(it3); + 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 - it3++; + itZip++; } } - if(!smthHappened) + else { - for(std::list::const_iterator it5=(*it1)->_sub_edges.begin();it5!=(*it1)->_sub_edges.end();it5++) + 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()); @@ -1312,14 +1321,20 @@ 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 (?? might there be others??) { - (*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)->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()); } } } -- 2.39.2