From: abn Date: Mon, 12 May 2014 14:36:37 +0000 (+0200) Subject: Intersection of 2D mesh with a 1D line: Intersect2DMeshWith1DLine X-Git-Tag: intersect_1D_OK~6 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=45bf6a4ab7cbb85ca99ccffc8aecf0b045552d75;p=tools%2Fmedcoupling.git Intersection of 2D mesh with a 1D line: Intersect2DMeshWith1DLine First working version. Still some doc and code re-factoring to be done. --- diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx index e3a738530..2d3a1d918 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx @@ -137,7 +137,7 @@ void ComposedEdge::initLocations() const } /** - * Reset the status of all edges (OUT, IN, ON) because they were potentially assignated + * Reset the status of all edges (OUT, IN, ON) because they were potentially assigned * by the previous candidate processing. */ void ComposedEdge::initLocationsWithOther(const ComposedEdge& other) const diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx index e014d5942..142a7cc8d 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx @@ -37,6 +37,11 @@ namespace INTERP_KERNEL class ElementaryEdge; class IteratorOnComposedEdge; + /** + * A set of quadratic or linear edges, described mainly by their connectivity + * The set is assumed to be connected, but not necessarily closed (i.e. not necessarily forming a closed polygon). + * Some methods however requires a closed form. + */ class ComposedEdge { friend class IteratorOnComposedEdge; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx index 45088dbef..4db25de03 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx @@ -223,10 +223,9 @@ void ElementaryEdge::fillGlobalInfoAbs2(const std::map& edgesThis, std::vector& addCoo, std::map mapAddCoo) const; INTERPKERNEL_EXPORT void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const; - INTERPKERNEL_EXPORT static ElementaryEdge *BuildEdgeFromCrudeDataArray(bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end); + INTERPKERNEL_EXPORT static ElementaryEdge *BuildEdgeFromStartEndDir(bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end); private: bool _direction; Edge *_ptr; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 5df47c7fc..255f9a080 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -408,6 +408,57 @@ void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const s } } +void QuadraticPolygon::appendEdgeFromCrudeDataArray2(const std::map& mapp, bool isQuad, + const int *nodalBg, const double *coords, + const int segId, const std::vector >& intersectEdges) +{ + if(!isQuad) + { + bool direct=true; + int edgeId=segId; + const std::vector& subEdge=intersectEdges[edgeId]; + std::size_t nbOfSubEdges=subEdge.size()/2; + for(std::size_t j=0;j& subEdge=intersectEdges[edgeId]; + std::size_t nbOfSubEdges=subEdge.size()/2; + if(colinearity) + { + for(std::size_t j=0;jdecrRef(); + } + st0->decrRef(); endd0->decrRef(); middle0->decrRef(); + } +} + + void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, int edgeId, const std::vector& subEdge, const std::map& mapp) { std::size_t nbOfSubEdges=subEdge.size()/2; @@ -415,7 +466,7 @@ void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size {//it is not a quadratic subedge Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second; Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second; - ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end); + ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end); pushBack(e); } else @@ -429,8 +480,10 @@ 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 + * This method builds a QP from the descending conn of a cell stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order, to give the * orientation of edge. + * The information stored in intersectEdges1 and colinear1 is also used to know which intermediate nodes (intersection points with the other polygon) should be considered + * on the edges. */ 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, @@ -530,7 +583,7 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, const int segId, 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) +{ + bool direct=true; + int edgeId=segId;//current edge id of pol2 + std::map >::const_iterator it1=alreadyExistingIn2.find(segId); + if(it1!=alreadyExistingIn2.end()) + { + const std::vector& edgesAlreadyBuilt=(*it1).second; + for(std::vector::const_iterator it3=edgesAlreadyBuilt.begin();it3!=edgesAlreadyBuilt.end();it3++) + { + Edge *ee=(*it3)->getPtr(); ee->incrRef(); + pushBack(new ElementaryEdge(ee,(*it3)->getDirection())); + } + return; + } + bool directos=colinear1[edgeId].empty(); + std::vector > > idIns1; + int 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]; + 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 + //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 + } + 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(); + appendEdgeFromCrudeDataArray2(mapp,isQuad,nodalBg,coords,segId,intersectEdges); + std::size_t newSz=_sub_edges.size(); + std::size_t zeSz=newSz-oldSz; + alreadyExistingIn2[segId].resize(zeSz); + std::list::const_reverse_iterator it5=_sub_edges.rbegin(); + for(std::size_t p=0;p& subEdge=intersectEdges[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++) + { + int 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; + offset2=0; + for(std::size_t k=0;k build new Edge instance + //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges); + Node *start=(*mapp.find(idBg)).second; + Node *end=(*mapp.find(idEnd)).second; + ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end); + pushBack(e); + alreadyExistingIn2[segId].push_back(e); + } + else + {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1 + ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)]; + Edge *ee=e->getPtr(); + ee->incrRef(); + ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11)); + pushBack(e2); + alreadyExistingIn2[segId].push_back(e2); + } + } + } +} + +/** + * @param edgeId current edge id of pol2 + */ +void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray(const int edgeId, 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) const +{ + std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1); + const std::vector& c=colinear1[edgeId]; + if(c.empty()) + return; + const std::vector& subEdge=intersectEdges[edgeId]; + std::size_t nbOfSubEdges=subEdge.size()/2; + // + int 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; + 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 + } +} + + /*! * Method expected to be called on pol2. Every params not suffixed by numbered are supposed to refer to pol2 (this). - * Method to find edges that are ON. + * Method to find edges that are ON, and mark them correspondingly on pol1. */ void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(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) const { + std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1); std::size_t nbOfSeg=std::distance(descBg,descEnd); for(std::size_t i=0;i& subEdge=intersectEdges[edgeId]; std::size_t nbOfSubEdges=subEdge.size()/2; // - std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1); int offset1=0; for(std::size_t j=0;j& 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, int offset, std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI) const { int nbOfNodesInPg=0; bool presenceOfQuadratic=presenceOfQuadraticEdge(); @@ -622,7 +825,7 @@ void QuadraticPolygon::appendCrudeData(const std::map for(std::list::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++) { INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf(); - node->unApplySimilarity(xBary,yBary,fact); + //node->unApplySimilarity(xBary,yBary,fact); addCoordsQuadratic.push_back((*node)[0]); addCoordsQuadratic.push_back((*node)[1]); conn.push_back(off+j); @@ -633,21 +836,26 @@ void QuadraticPolygon::appendCrudeData(const std::map } /*! - * This method make the hypothesis that 'this' and 'other' are splited at the minimum into edges that are fully IN, OUT or ON. + * This method make the hypothesis that 'this' and 'other' are split at the minimum into edges that are fully IN, OUT or ON. * This method returns newly created polygons in 'conn' and 'connI' and the corresponding ids ('idThis','idOther') are stored respectively into 'nbThis' and 'nbOther'. - * @param [in,out] edgesThis, parameter that keep informed the caller abount the edges in this not shared by the result of intersection of \a this with \a other - * @param [in,out] edgesBoundaryOther, parameter that strores 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) +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) { double xBaryBB, yBaryBB; double fact=normalizeExt(&other, xBaryBB, yBaryBB); - //Locate 'this' relative to 'other' + //Locate 'this' relative to 'other' (edges of 'this', aka 'pol1' are marked as IN or OUT) other.performLocatingOperationSlow(*this); // without any assumption - std::vector res=buildIntersectionPolygons(other,*this); + std::vector res = BuildIntersectionPolygons(other,*this); + unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact); + for(std::vector::iterator it=res.begin();it!=res.end();it++) { - (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI); + (*it)->appendCrudeData(mapp,offset,addCoordsQuadratic,conn,connI); INTERP_KERNEL::IteratorOnComposedEdge it1(*it); for(it1.first();!it1.finished();it1.next()) { @@ -666,9 +874,64 @@ void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set & zip_list, + std::set & edgesThis, std::set & edgesBoundaryOther, + const std::map& mapp, + const int idThis, const std::vector > & mapZipTo2, + int offset, std::vector& addCoordsQuadratic, + std::vector& conn,std::vector& connI, + std::vector& nbThis, std::vector& nbOther, std::vector& nbOtherI) +{ + std::vector res; + if(!zip_list.empty()) + ClosePolygonsSimple(pol1, zip_list, res, mapZipTo2, nbOther,nbOtherI); + + for(std::vector::iterator it=res.begin();it!=res.end();it++) + { + (*it)->appendCrudeData(mapp,offset,addCoordsQuadratic,conn,connI); + INTERP_KERNEL::IteratorOnComposedEdge it1(*it); + for(it1.first();!it1.finished();it1.next()) + { + Edge *e=it1.current()->getPtr(); + if(edgesThis.find(e)!=edgesThis.end()) + edgesThis.erase(e); + else + { + if(edgesBoundaryOther.find(e)!=edgesBoundaryOther.end()) + edgesBoundaryOther.erase(e); + else + edgesBoundaryOther.insert(e); + } + } + // Mapping to pol1 (mapping to pol2s has been handled in ClosePolygonsSimple). + nbThis.push_back(idThis); + delete (*it); + } +} + +/** @sa ComposedEdge::initLocationsWithOther() + */ +void QuadraticPolygon::initLocationsWithSeveralOthers(const std::vector & others) const +{ + std::set s1,s2; + for(std::list::const_iterator it1=_sub_edges.begin();it1!=_sub_edges.end();it1++) + s1.insert((*it1)->getPtr()); + std::vector::const_iterator it_qp; + for(it_qp=others.begin(); it_qp!=others.end(); it_qp++) + for(std::list::const_iterator it2=(*it_qp)._sub_edges.begin();it2!=(*it_qp)._sub_edges.end();it2++) + s2.insert((*it2)->getPtr()); + initLocations(); + for(it_qp=others.begin(); it_qp!=others.end(); it_qp++) + (*it_qp).initLocations(); + std::vector s3; + std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(),std::back_insert_iterator< std::vector >(s3)); + for(std::vector::const_iterator it3=s3.begin();it3!=s3.end();it3++) + (*it3)->declareOn(); +} + + /*! * 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. @@ -876,7 +1139,7 @@ std::vector QuadraticPolygon::intersectMySelfWith(const Quad SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits); //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done. performLocatingOperation(cpyOfOther); - return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther); + return BuildIntersectionPolygons(cpyOfThis,cpyOfOther); } /*! @@ -942,36 +1205,35 @@ void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol2) const } } -void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) const +void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol1) const { - IteratorOnComposedEdge it(&pol2); + IteratorOnComposedEdge it(&pol1); for(it.first();!it.finished();it.next()) { ElementaryEdge *cur=it.current(); - cur->locateFullyMySelfAbsolute(*this); + cur->locateFullyMySelfAbsolute(*this); // *this=pol2=other } } /*! * Given 2 polygons 'pol1' and 'pol2' (localized) the resulting polygons are returned. * - * this : pol2 simplified. * @param pol1 pol1 split. * @param pol2 pol2 split. */ -std::vector QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const +std::vector QuadraticPolygon::BuildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) { std::vector ret; std::list pol2Zip=pol2.zipConsecutiveInSegments(); if(!pol2Zip.empty()) - closePolygons(pol2Zip,pol1,ret); + ClosePolygons(pol2Zip,pol1,pol2,ret); else {//borders of pol2 do not cross pol1,and pol2 borders are outside of pol1. That is to say, either pol2 and pol1 //do not overlap or pol1 is fully inside pol2. So in the first case no intersection, in the other case //the intersection is pol1. ElementaryEdge *e1FromPol1=pol1[0]; TypeOfEdgeLocInPolygon loc=FULL_ON_1; - loc=e1FromPol1->locateFullyMySelf(*this,loc); + loc=e1FromPol1->locateFullyMySelf(pol2,loc); if(loc==FULL_IN_1) ret.push_back(new QuadraticPolygon(pol1)); } @@ -1014,14 +1276,61 @@ std::list QuadraticPolygon::zipConsecutiveInSegments() const return ret; } +std::list QuadraticPolygon::ZipConsecutiveSegments2(const std::vector & candidates2, + const std::vector & pol2s, + std::vector > & mapZipTo2) +{ + int prevI = -2; + TypeOfEdgeLocInPolygon prevLoc = FULL_UNKNOWN; + std::list ret; + std::vector::const_iterator it2; int ii=0; + QuadraticPolygon *tmpQp = 0; + for(it2 = candidates2.begin(), ii = 0; it2 != candidates2.end(); it2++, ii++) + { + IteratorOnComposedEdge it_ii(const_cast(static_cast(&pol2s[ii]))); // hmmm ... :-) + bool firstEdge = true; + for (; !it_ii.finished(); it_ii.next()) + { + if(it_ii.current()->getLoc() == FULL_IN_1) + { + // only keep consecutive items from mesh2 (consecutive seg IDs and consecutive IN pieces) + if(prevLoc != FULL_IN_1 || (firstEdge && *it2 != prevI+1)) + { + // the start node of a new QP has to be ON (otherwise we hit a piece of line starting in the middle of a cell) + if (it_ii.current()->getStartNode()->getLoc() == ON_1) + { + tmpQp = new QuadraticPolygon; + ret.push_back(tmpQp); + mapZipTo2.push_back(std::vector(0)); + } + else + continue; // firstEdge remains to its current value + } + tmpQp->pushBack(it_ii.current()->clone()); + // Mapping - one quadratic polygon in ret is linked to several 1D cells: + if (std::find(mapZipTo2.back().begin(),mapZipTo2.back().end(), *it2) == mapZipTo2.back().end()) + mapZipTo2.back().push_back(*it2); + } + prevLoc = it_ii.current()->getLoc(); + firstEdge = false; + } + prevI = *it2; + } + // Now clean line pieces not ending ON: + for(std::list::iterator it3 = ret.begin(); it3 != ret.end(); it3++) + if ((*it3)->back()->getEndNode()->getLoc() != ON_1) + ret.erase(it3); + + return ret; +} + /*! - * 'this' should be considered as pol2Simplified. - * @param pol2zip is a list of set of edges (openned polygon) coming from split polygon 2. + * @param pol2zip is a list of set of edges (=an opened polygon) coming from split polygon 2. * @param pol1 is split pol1. * @param results the resulting \b CLOSED polygons. */ -void QuadraticPolygon::closePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, - std::vector& results) const +void QuadraticPolygon::ClosePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2, + std::vector& results) { bool directionKnownInPol1=false; bool directionInPol1; @@ -1036,7 +1345,7 @@ void QuadraticPolygon::closePolygons(std::list& pol2Zip, con } if(!directionKnownInPol1) { - if(!(*iter)->amIAChanceToBeCompletedBy(pol1,*this,directionInPol1)) + if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol1)) { delete *iter; iter=pol2Zip.erase(iter); continue; } else directionKnownInPol1=true; @@ -1052,10 +1361,89 @@ void QuadraticPolygon::closePolygons(std::list& pol2Zip, con } } +void QuadraticPolygon::ClosePolygonsSimple(const QuadraticPolygon& pol1, const std::list & zip_list, + std::vector& results,const std::vector > & mapZipTo2, + std::vector& nbOther, std::vector& nbOtherI) +{ + // Register all edges of pol1 as 'alive' to start with. They will be consumed as they are + // added to the newly formed cells. Edges of pol2s (in zip_list) might be used more than once. + std::set edges1_alive; + IteratorOnComposedEdge it(const_cast(&pol1)); + for(; !it.finished(); it.next()) + edges1_alive.insert(it.current()); + + int inf_loop_detect, i; + std::list::const_iterator itt; + for (inf_loop_detect = pol1.recursiveSize() + 1, itt = zip_list.begin(); itt != zip_list.end(); itt++) + inf_loop_detect += 2*(*itt)->recursiveSize(); + + for (i = 0, it.first(); !edges1_alive.empty() && i <= inf_loop_detect;) + { + QuadraticPolygon * qp = new QuadraticPolygon(); + results.push_back(qp); + // Mapping initialization + int nbElemFrom2 = 0; + // Start anywhere valid (=alive) on pol1 + ElementaryEdge * startE; + for (startE = it.current(); edges1_alive.find(it.current()) == edges1_alive.end(); it.nextLoop(), startE = it.current()); + + // Close the new cell: + do + { + ElementaryEdge *tmp_e = it.current(); + qp->pushBack(tmp_e->clone()); + edges1_alive.erase(tmp_e); + Node * nodeToTest = tmp_e->getEndNode(); + it.nextLoop(); i++; + // Try to match the end node of the latest added edge from pol1 with one of the node from + // one of the ComposedEdge in zip_list + std::list::const_iterator ret; + bool direction; + int idx; // for mapping + ret = CheckInList2(nodeToTest, zip_list, direction, idx); + if (ret == zip_list.end()) + continue; + size_t oldSz = nbOther.size(); + size_t addSz = std::distance(mapZipTo2[idx].begin(), mapZipTo2[idx].end()); + nbOther.resize(oldSz+addSz); + std::copy(mapZipTo2[idx].begin(), mapZipTo2[idx].end(), &nbOther[oldSz]); + nbElemFrom2 += addSz; + if (!direction) + (*ret)->reverse(); // TODO: discuss with Anthony - API of IteratorOnComposedEdge should allow an easy reverse looping? + // Switch to zip_list + IteratorOnComposedEdge it_pol2(*ret); + ElementaryEdge *tmp_e2; + for(/* it_pol2.first() */; !it_pol2.finished(); it_pol2.next(), i++) + { + tmp_e2 = it_pol2.current(); + qp->pushBack(tmp_e2->clone()); + } + // search where we landed on pol1 to resume looping on it. + if (!edges1_alive.empty()) + { + int j; + for(j = 0; it.current()->getStartNode() != tmp_e2->getEndNode() && j < inf_loop_detect; j++) // better way of doing this? + it.nextLoop(); + if (j >= inf_loop_detect) // could be a bit more precise and stop before ... + throw Exception("QuadraticPolygon::ClosePolygonsSimple() - Internal error - 1D intersection was looping infinitely!"); + // sanity check - if we are not done with the cell completion, we should land + // on a live edge. + if(edges1_alive.find(it.current()) == edges1_alive.end() && it.current() != startE) + throw Exception("QuadraticPolygon::ClosePolygonsSimple() - Internal error ! Should never happen."); + } + } + while(it.current() != startE && !edges1_alive.empty() && i < inf_loop_detect); + // Complete mapping + nbOtherI.push_back(nbElemFrom2+nbOtherI.back()); + } + if (i >= inf_loop_detect) + throw Exception("QuadraticPolygon::ClosePolygonsSimple() - Internal error - 1D intersection was looping infinitely!"); +} + /*! * 'this' is expected to be set of edges (not closed) of pol2 split. */ -bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction) +bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction) { IteratorOnComposedEdge it(const_cast(&pol1Splitted)); bool found=false; @@ -1069,7 +1457,7 @@ bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Spl it.next(); } if(!found) - throw Exception("Internal error : polygons uncompatible each others. Should never happend"); + throw Exception("Internal error: polygons incompatible with each others. Should never happen!"); //Ok we found correspondance between this and pol1. Searching for right direction to close polygon. ElementaryEdge *e=_sub_edges.back(); if(e->getLoc()==FULL_ON_1) @@ -1099,7 +1487,7 @@ bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Spl } /*! - * This method fills as much as possible 'this' (part of pol2 split) with edges of 'pol1Splitted'. + * This method fills as much as possible 'this' (a sub-part of pol2 split) with edges of 'pol1Splitted'. */ std::list::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted, std::list::iterator iStart, @@ -1151,6 +1539,28 @@ std::list::iterator QuadraticPolygon::CheckInList(Node *n, s return iEnd; } +std::list::const_iterator QuadraticPolygon::CheckInList2(Node *n, const std::list & zip_list, + bool & direction, int & index) +{ + direction = true; + index = 0; + for(std::list::const_iterator iter = zip_list.begin(); iter!=zip_list.end(); iter++, index++) + { + if ((*iter)->front()->getStartNode() == n) + return iter; + else + { + if ((*iter)->back()->getEndNode() == n) + { + direction = false; + return iter; + } + } + } + return zip_list.end(); +} + + 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) { @@ -1166,7 +1576,7 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: std::list pol1Zip; 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); + pol1.appendCrudeData(mapp,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1); return ; } while(!notUsedInPol1L.empty()) @@ -1272,7 +1682,7 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std:: { if((*it1)->getStartNode()==(*it1)->getEndNode()) { - (*it1)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1); + (*it1)->appendCrudeData(mapp,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++) delete *it6; delete *it1; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx index 670ac6745..2935032b2 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx @@ -71,13 +71,31 @@ namespace INTERP_KERNEL const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector >& intersectEdges1, const std::vector< std::vector >& colinear1, std::map >& alreadyExistingIn2); + INTERPKERNEL_EXPORT void buildFromCrudeDataArray3(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, const int segId, 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); + INTERPKERNEL_EXPORT void updateLocOfEdgeFromCrudeDataArray(const int edgeId, 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) const; INTERPKERNEL_EXPORT void updateLocOfEdgeFromCrudeDataArray2(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) const; INTERPKERNEL_EXPORT void appendEdgeFromCrudeDataArray(std::size_t edgeId, const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector >& intersectEdges); + INTERPKERNEL_EXPORT void appendEdgeFromCrudeDataArray2(const std::map& mapp, bool isQuad, + const int *nodalBg, const double *coords, + const int segId, const std::vector >& intersectEdges); INTERPKERNEL_EXPORT void appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, int edgeId, const std::vector& subEdge, const std::map& mapp); - INTERPKERNEL_EXPORT void appendCrudeData(const std::map& mapp, double xBary, double yBary, double fact, int offset, std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI) const; + INTERPKERNEL_EXPORT void appendCrudeData(const std::map& mapp, int offset, std::vector& addCoordsQuadratic, std::vector& conn, std::vector& connI) const; INTERPKERNEL_EXPORT void 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& nb1, std::vector& nb2); + INTERPKERNEL_EXPORT static void BuildPartitionFromZipList(const QuadraticPolygon & pol1, const std::list & zip_list, + std::set & edgesThis, std::set & edgesBoundaryOther, + const std::map& mapp, + const int idThis, const std::vector > & mapZipTo2, + int offset, std::vector& addCoordsQuadratic, + std::vector& conn,std::vector& connI, + std::vector& nbThis, std::vector& nbOther, std::vector& nbOtherI); + INTERPKERNEL_EXPORT void initLocationsWithSeveralOthers(const std::vector & others) const; // INTERPKERNEL_EXPORT double intersectWith(const QuadraticPolygon& other) const; INTERPKERNEL_EXPORT double intersectWith(const QuadraticPolygon& other, double* barycenter) const; @@ -89,23 +107,31 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT void performLocatingOperation(QuadraticPolygon& pol2) const; INTERPKERNEL_EXPORT void performLocatingOperationSlow(QuadraticPolygon& pol2) const; INTERPKERNEL_EXPORT static void SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits); - INTERPKERNEL_EXPORT std::vector buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const; - INTERPKERNEL_EXPORT bool amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction); + INTERPKERNEL_EXPORT static std::vector BuildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2); + INTERPKERNEL_EXPORT bool haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction); INTERPKERNEL_EXPORT static void 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); + INTERPKERNEL_EXPORT static std::list ZipConsecutiveSegments2(const std::vector & candidates2, const std::vector & pol2s, + std::vector > & mapZipTo2); protected: std::list zipConsecutiveInSegments() const; void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; - void closePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, std::vector& results) const; + static void ClosePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2, std::vector& results); + static void ClosePolygonsSimple(const QuadraticPolygon& pol1, const std::list & zip_list, + std::vector& results, const std::vector > & mapZipTo2, + std::vector& nbOther, std::vector& nbOtherI); template static void UpdateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2, const EDGES *e1, const EDGES *e2); + std::list::iterator fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted, std::list::iterator iStart, std::list::iterator iEnd, bool direction); static std::list::iterator CheckInList(Node *n, std::list::iterator iStart, std::list::iterator iEnd); + static std::list::const_iterator CheckInList2(Node *n, const std::list & zip_list, + bool & direction, int & index); }; } diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index e0113eba9..9842c1aa1 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -4250,7 +4250,7 @@ namespace ParaMEDMEM } /** - * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI). + * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI). */ void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo, const int *desc1Bg, const int *desc1End, const std::vector >& intesctEdges1, @@ -4271,6 +4271,27 @@ namespace ParaMEDMEM } } } + + /** + * Same as above for a 1D cell made of a single segment + * TODO: factorize with the above + */ + void MEDCouplingUMeshBuildQPFromMesh4(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo, + const int segId, const std::vector >& intesctEdges1, + /*output*/std::map& mapp, std::map& mappRev) + { + for(std::vector::const_iterator it1=intesctEdges1[segId].begin();it1!=intesctEdges1[segId].end();it1++) + { + std::map::const_iterator it=mappRev.find(*it1); + if(it==mappRev.end()) + { + INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo); + mapp[node]=*it1; + mappRev[*it1]=node; + } + } + } + } /// @endcond @@ -8685,7 +8706,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New("Intersect2D",2); MEDCouplingAutoRefCountObjectPtr conn=DataArrayInt::New(); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer()); MEDCouplingAutoRefCountObjectPtr connI=DataArrayInt::New(); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer()); - MEDCouplingAutoRefCountObjectPtr c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer()); + MEDCouplingAutoRefCountObjectPtr c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer()); MEDCouplingAutoRefCountObjectPtr c2=DataArrayInt::New(); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer()); ret->setConnectivity(conn,connI,true); ret->setCoords(coo); @@ -8693,6 +8714,73 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 return ret.retn(); } +MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, + double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2) +{ + m1->checkFullyDefined(); + m2->checkFullyDefined(); + if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=1 || m2->getSpaceDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works on unstructured meshes m1 with (meshdim, spacedim) = (2,2) and m2 with (meshdim, spacedim) = (1,2)!"); + + // Check that m2 is a line (and not a more complex 1D mesh) - each point is used at most by 2 segments: + DataArrayInt * d = DataArrayInt::New(), * dI = DataArrayInt::New(); + DataArrayInt *rD = DataArrayInt::New(), * rDI = DataArrayInt::New(); + MEDCouplingUMesh * m_points; + m_points = m2->buildDescendingConnectivity(d, dI, rD, rDI); + DataArrayInt * dsi = rDI->deltaShiftIndex(); + MEDCouplingAutoRefCountObjectPtr dsii = dsi->getIdsNotInRange(0,3); + m_points->decrRef(); d->decrRef(); dI->decrRef(); rD->decrRef(); rDI->decrRef(); dsi->decrRef(); + if (dsii->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine only work with mesh2 being a (piecewise) connected line!"); + + // Step 1: compute all edge intersections (new nodes) + std::vector< std::vector > intersectEdge1, colinear2, subDiv2; + MEDCouplingUMesh *m1Desc=0; + DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0; + std::vector addCoo,addCoordsQuadratic; // coordinates of newly created nodes + INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; + INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps; + IntersectDescending2DMeshWith1DMesh(m1,m2,eps,intersectEdge1,colinear2, subDiv2, + m1Desc,desc1,descIndx1,revDesc1,revDescIndx1, + addCoo); + revDesc1->decrRef(); revDescIndx1->decrRef(); + MEDCouplingAutoRefCountObjectPtr dd1(desc1),dd2(descIndx1); + MEDCouplingAutoRefCountObjectPtr dd5(m1Desc); + + // Step 2: re-order newly created nodes according to the ordering found in m2 + std::vector< std::vector > intersectEdge2; + BuildIntersectEdges(m1Desc,m2,addCoo,subDiv2,intersectEdge2); + subDiv2.clear(); dd5=0; + + // Step 3: this is where we have a significant difference with Intersect2DMeshes + std::vector cr,crI; //no DataArrayInt because interface with Geometric2D + std::vector cNb1,cNb2,cNbI2; //no DataArrayInt because interface with Geometric2D + cNbI2.push_back(0); + Build2DCellsFrom1DCut(eps,m1,desc1->getConstPointer(),descIndx1->getConstPointer(),intersectEdge1, + colinear2,m2,intersectEdge2,addCoo, + /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2,cNbI2); + + // Step 4: Prepare final result: + MEDCouplingAutoRefCountObjectPtr addCooDa=DataArrayDouble::New(); + addCooDa->alloc((int)(addCoo.size())/2,2); + std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer()); + MEDCouplingAutoRefCountObjectPtr addCoordsQuadraticDa=DataArrayDouble::New(); + addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2); + std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer()); + std::vector coordss(4); + coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa; + MEDCouplingAutoRefCountObjectPtr coo=DataArrayDouble::Aggregate(coordss); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New("Intersect2D",2); + MEDCouplingAutoRefCountObjectPtr conn=DataArrayInt::New(); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer()); + MEDCouplingAutoRefCountObjectPtr connI=DataArrayInt::New(); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer()); + MEDCouplingAutoRefCountObjectPtr c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer()); + MEDCouplingAutoRefCountObjectPtr c2=DataArrayInt::New(); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer()); + MEDCouplingAutoRefCountObjectPtr cI2=DataArrayInt::New(); cI2->alloc((int)cNbI2.size(),1); std::copy(cNbI2.begin(),cNbI2.end(),cI2->getPointer()); + ret->setConnectivity(conn,connI,true); + ret->setCoords(coo); + cellNb1=c1.retn(); cellNb2=c2.retn();cellNbI2=cI2.retn(); + return ret.retn(); +} /** * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the @@ -8735,7 +8823,7 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ); // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects: MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev); - // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes. + // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes. pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1, desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1); // @@ -8753,7 +8841,8 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]]; const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2); // Complete mapping with elements coming from the current cell it2 in mesh2: - MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev); + MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2], + desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev); // pol2 is the new QP in the final merged result. pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2, pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare); @@ -8785,6 +8874,119 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo } } +void MEDCouplingUMesh::Build2DCellsFrom1DCut(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, + const MEDCouplingUMesh *m2, const std::vector >& intesctEdges2, + const std::vector& addCoords, + std::vector& addCoordsQuadratic, std::vector& cr, std::vector& crI, + std::vector& cNb1, std::vector& cNb2, std::vector& cNbI2) +{ + static const int SPACEDIM=2; + const double *coo1=m1->getCoords()->getConstPointer(); + const int *conn1=m1->getNodalConnectivity()->getConstPointer(); + const int *connI1=m1->getNodalConnectivityIndex()->getConstPointer(); + int offset1=m1->getNumberOfNodes(); + const double *coo2=m2->getCoords()->getConstPointer(); + const int *conn2=m2->getNodalConnectivity()->getConstPointer(); + const int *connI2=m2->getNodalConnectivityIndex()->getConstPointer(); + int offset2=offset1+m2->getNumberOfNodes(); + int offset3=offset2+((int)addCoords.size())/2; + MEDCouplingAutoRefCountObjectPtr bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree()); + const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); + // Here a BBTree on 1D-segments from the tool mesh, not on segments: + BBTree myTree(bbox2,0,0,m2->getNumberOfCells(),eps); + int ncell1=m1->getNumberOfCells(); + crI.push_back(0); + // for all 2D cells in base mesh, find intersecting segments + for(int i=0;i candidates2; + myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2); + // Sort the candidates: this will help reconstructing connected lines traversing cell(i) in m1. + std::sort( candidates2.begin(), candidates2.end() ); + std::map mapp; + std::map mappRev; + // conversion to the geometric DS format: + INTERP_KERNEL::QuadraticPolygon pol1; + INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]]; + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ); + // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects: + MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1, + /* output-->*/mapp,mappRev); + // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes. + pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1, + desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1); + // + std::set edges1;// store all edges of pol1 that are NOT consumed by intersect cells. + // If any after the for-loop over candidates2 -> it remains unhandled parts of pol1 which should appear in the result (sse computeResidual below) + std::set edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1. + INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1); + for(it1.first();!it1.finished();it1.next()) + edges1.insert(it1.current()->getPtr()); + // + std::map > edgesIn2ForShare; // common edges + std::vector pol2s(candidates2.size()); + std::vector::iterator it2; int ii=0; + for(it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) + { + INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]]; + const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2); + // Complete mapping with elements coming from the current segment it2 in linear mesh2: + MEDCouplingUMeshBuildQPFromMesh4(coo1,offset1,coo2,offset2,addCoords,*it2,intesctEdges2, + /* output */mapp,mappRev); + // pol2 is the new QP in the final merged result. + pol2s[ii].buildFromCrudeDataArray3(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,*it2,intesctEdges2, + pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare); + } + pol1.initLocationsWithSeveralOthers(pol2s); + for(it2 = candidates2.begin(), ii = 0; it2 != candidates2.end(); it2++, ii++) + { + pol2s[ii].updateLocOfEdgeFromCrudeDataArray(*it2 ,intesctEdges2,pol1, desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2); + double xBaryBB, yBaryBB; + double fact=pol1.normalizeExt(&pol2s[ii], xBaryBB, yBaryBB); + //Locate pol2s[ii] relative to pol1 (this is done the other way around in the 2D/2D intersection) + pol1.performLocatingOperationSlow(pol2s[ii]); + pol1.unApplyGlobalSimilarityExt(pol2s[ii],xBaryBB,yBaryBB,fact); + } + // At this point all edges in (all elements of) pol2s are fully localized. + // Zip consecutive IN segments from the list of candidates, and discard the pieces which are ending nowhere (e.g. an IN polyline + // with an end-point in the middle of the cell). + std::vector > mapZipTo2; + std::list zip_list = INTERP_KERNEL::QuadraticPolygon::ZipConsecutiveSegments2(candidates2, pol2s, mapZipTo2); + + // Now the core of the algo - main output is in cr, crI, cNb1, cNb2 and cNbI2. + INTERP_KERNEL::QuadraticPolygon::BuildPartitionFromZipList(pol1, zip_list, edges1, edgesBoundary2, mapp, + i, mapZipTo2, + offset3,addCoordsQuadratic,cr,crI, cNb1,cNb2, cNbI2); + // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched + // by m2 but that we still want to keep in the final result. + if(!edges1.empty()) + { + size_t oldSz = cNb2.size(); + try + { + INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2); + } + catch(INTERP_KERNEL::Exception& e) + { + std::ostringstream oss; oss << "Error when computing residual of cell #" << i << " in source/m1 mesh ! Maybe the neighbours of this cell in mesh are not well connected !\n" << "The deep reason is the following : " << e.what(); + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + // ComputeResidual has pushed some -1 at the back of cNb2. Replace this with void entries in cNbI2: + // TODO: discuss the format with Anthony. + size_t newSz = cNb2.size(); + cNb2.resize(oldSz); + int val = (oldSz == 0) ? 0 : cNbI2.back(); + for(size_t sss=0; sss::const_iterator it=mappRev.begin();it!=mappRev.end();it++) + (*it).second->decrRef(); + for (std::list::const_iterator itt = zip_list.begin(); itt != zip_list.end(); itt++) + delete(*itt); + } +} + + void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map& m, int forbVal0, int forbVal1, std::vector& isect) { std::map::const_iterator it(m.find(n)); @@ -9080,7 +9282,7 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c intersectEdge1.resize(nDescCell1); colinear2.resize(nDescCell2); subDiv2.resize(nDescCell2); - BBTree myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps); + BBTree myTree(bbox2,0,0,nDescCell2,-eps); std::vector candidates1(1); int offset1=m1->getNumberOfNodes(); @@ -9118,6 +9320,67 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef(); } +void MEDCouplingUMesh::IntersectDescending2DMeshWith1DMesh(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, + std::vector< std::vector >& intersectEdge1, std::vector< std::vector >& colinear2, std::vector< std::vector >& subDiv2, + MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1, + std::vector& addCoo) +{ + static const int SPACEDIM=2; + // Build desc connectivity of first mesh only + desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); + revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1); + m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1); + MEDCouplingAutoRefCountObjectPtr dd9(m1Desc); + const int *c1=m1Desc->getNodalConnectivity()->getConstPointer(); + const int *ci1=m1Desc->getNodalConnectivityIndex()->getConstPointer(); + + // Build BB tree of all edges in the tool mesh (second mesh) + MEDCouplingAutoRefCountObjectPtr bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree()); + const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); + int nDescCell1=m1Desc->getNumberOfCells(); + int nCell2=m2->getNumberOfCells(); + intersectEdge1.resize(nDescCell1); + colinear2.resize(nCell2); + subDiv2.resize(nCell2); + BBTree myTree(bbox2,0,0,nCell2,-eps); + + std::vector candidates1(1); + int offset1=m1->getNumberOfNodes(); + int offset2=offset1+m2->getNumberOfNodes(); + for(int i=0;i candidates2; // edges of mesh2 candidate for intersection + myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2); + if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1 + { + std::map map1,map2; + // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format + INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2,candidates2,map2); + candidates1[0]=i; + INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1); + // This following part is to avoid that some removed nodes (for example due to a merge between pol1 and pol2) are replaced by a newly created one + // This trick guarantees that Node * are discriminant (i.e. form a unique identifier) + std::set nodes; + pol1->getAllNodes(nodes); pol2->getAllNodes(nodes); + std::size_t szz(nodes.size()); + std::vector< MEDCouplingAutoRefCountObjectPtr > nodesSafe(szz); + std::set::const_iterator itt(nodes.begin()); + for(std::size_t iii=0;iiiincrRef(); nodesSafe[iii]=*itt; } + // end of protection + // Performs egde cutting: + pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo); + delete pol2; + delete pol1; + } + else + intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i+1]); + } + m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef(); +} + + /*! * This method performs the 2nd step of Partition of 2D mesh. * This method has 4 inputs : diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index e45707c60..8a3c98a27 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -240,6 +240,8 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT static void ComputeVecAndPtOfFace(double eps, const double *coords, const int *begin, const int *end, double *v, double *p); MEDCOUPLING_EXPORT static void TryToCorrectPolyhedronOrientation(int *begin, int *end, const double *coords); MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2); + MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, + double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2); MEDCOUPLING_EXPORT static bool BuildConvexEnvelopOf2DCellJarvis(const double *coords, const int *nodalConnBg, const int *nodalConnEnd, DataArrayInt *nodalConnecOut); MEDCOUPLING_EXPORT static bool RemoveIdsFromIndexedArrays(const int *idsToRemoveBg, const int *idsToRemoveEnd, DataArrayInt *arr, DataArrayInt *arrIndx, int offsetForRemoval=0); MEDCOUPLING_EXPORT static void ExtractFromIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn, @@ -311,11 +313,20 @@ namespace ParaMEDMEM MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1, std::vector& addCoo, MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2); + static void IntersectDescending2DMeshWith1DMesh(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, + std::vector< std::vector >& intersectEdge1, std::vector< std::vector >& colinear2, std::vector< std::vector >& subDiv2, + MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1, + std::vector& addCoo); static void BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, const std::vector& addCoo, const std::vector< std::vector >& subDiv, std::vector< std::vector >& intersectEdge); static void BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector >& intesctEdges2, const std::vector& addCoords, std::vector& addCoordsQuadratic, std::vector& cr, std::vector& crI, std::vector& cNb1, std::vector& cNb2); + static void Build2DCellsFrom1DCut(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, + const MEDCouplingUMesh *m2, const std::vector >& intesctEdges2, + const std::vector& addCoords, + std::vector& addCoordsQuadratic, std::vector& cr, std::vector& crI, + std::vector& cNb1, std::vector& cNb2, std::vector& cNbI2); static void AssemblyForSplitFrom3DCurve(const std::vector& cut3DCurve, std::vector& nodesOnPlane, const int *nodal3DSurf, const int *nodalIndx3DSurf, const int *nodal3DCurve, const int *nodalIndx3DCurve, const int *desc, const int *descIndx, std::vector< std::pair >& cut3DSurf); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx index d26230463..26b1ef7e7 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx @@ -2057,3 +2057,79 @@ void MEDCouplingBasicsTest5::testSimplexize3() // m->decrRef(); } + +/** + * Two square cells intersected by a vertical line crossing only one cell. + */ +void MEDCouplingBasicsTest5::testIntersect2DMeshWith1DLine1() +{ + MEDCouplingCMesh *m1c = MEDCouplingCMesh::New(); + DataArrayDouble *coordX = DataArrayDouble::New(); + const double arrX[3] = {-1., 1., 2}; + coordX->alloc(3,1); + std::copy(arrX,arrX+3,coordX->getPointer()); + m1c->setCoordsAt(0,coordX); + DataArrayDouble *coordY = DataArrayDouble::New(); + const double arrY[4] = {0., 2.}; + coordY->alloc(2,1); + std::copy(arrY,arrY+2,coordY->getPointer()); + m1c->setCoordsAt(1,coordY); + MEDCouplingUMesh *m1 = m1c->buildUnstructured(); + + // A simple line: + MEDCouplingUMesh *m2 = MEDCouplingUMesh::New("bla", 1); + const int conn2A[7] = {INTERP_KERNEL::NORM_SEG2,0,1,INTERP_KERNEL::NORM_SEG3,1,2,3}; + const int connI2A[3] = {0,3,7}; + const double coo2A[8] = {0.,-1.0, 0.,1., 0.,3., 0.5,2.2}; + DataArrayDouble *coord2 = DataArrayDouble::New(); + DataArrayInt *conn2 = DataArrayInt::New(); + DataArrayInt *connI2 = DataArrayInt::New(); + coord2->alloc(4,2); conn2->alloc(7,1); connI2->alloc(3, 1); + std::copy(coo2A ,coo2A+8, coord2->getPointer()); + std::copy(conn2A ,conn2A+7, conn2->getPointer()); + std::copy(connI2A,connI2A+3, connI2->getPointer()); + m2->setCoords(coord2); + m2->setConnectivity(conn2, connI2); + + // End of construction of input meshes m1bis and m2 -> start of specific part of the test + DataArrayInt * map1 = 0, * map2 = 0, * mapI2 = 0; + MEDCouplingUMesh *m3 = MEDCouplingUMesh::Intersect2DMeshWith1DLine(m1, m2, 1e-10, map1, map2, mapI2); + bool areNodesMerged; int newNbNodes; + DataArrayInt * d_tmp = m3->mergeNodes(1.0e-8, areNodesMerged, newNbNodes); + d_tmp->decrRef(); + CPPUNIT_ASSERT_EQUAL(3,m3->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(20,m3->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(2,m3->getSpaceDimension()); + + const int exp1[3] = {0,0,1}; + const int exp2[4] = {0,1,0,1}; + const int expI2[4] = {0,2,4,4}; + CPPUNIT_ASSERT_EQUAL(3, map1->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(4, map2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(4, mapI2->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(exp1,exp1+3, map1->getConstPointer())); + CPPUNIT_ASSERT(std::equal(exp2,exp2+4, map2->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expI2,expI2+4, mapI2->getConstPointer())); + + const int expConn[27] = {32, 1, 10, 7, 11, 4, 12, 13, 14, 15, 16, 32, 10, 0, 3, 11, 7, 17, 18, 19, 14, 13, 5, 2, 1, 4, 5}; + const int expConnI[4] = {0, 11, 22, 27}; + CPPUNIT_ASSERT_EQUAL(27, m3->getNodalConnectivity()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(4, m3->getNodalConnectivityIndex()->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expConn,expConn+27, m3->getNodalConnectivity()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expConnI,expConnI+4, m3->getNodalConnectivityIndex()->getConstPointer())); + + map1->decrRef(); + map2->decrRef(); + mapI2->decrRef(); + m3->decrRef(); + // + m2->decrRef(); + m1->decrRef(); + coordX->decrRef(); + coordY->decrRef(); + m1c->decrRef(); + coord2->decrRef(); + conn2->decrRef(); + connI2->decrRef(); +} + diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx index 14b222900..87dcf1051 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx @@ -76,6 +76,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( testDAIBuildSubstractionOptimized1 ); CPPUNIT_TEST( testDAIIsStrictlyMonotonic1 ); CPPUNIT_TEST( testSimplexize3 ); + CPPUNIT_TEST( testIntersect2DMeshWith1DLine1 ); CPPUNIT_TEST_SUITE_END(); public: void testUMeshTessellate2D1(); @@ -118,6 +119,7 @@ namespace ParaMEDMEM void testDAIBuildSubstractionOptimized1(); void testDAIIsStrictlyMonotonic1(); void testSimplexize3(); + void testIntersect2DMeshWith1DLine1(); }; } diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 1f58aec18..6cab661aa 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -14724,9 +14724,99 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(m.getNodalConnectivityIndex().isEqual(DataArrayInt([0,5]))) pass + def testIntersect2DMeshWith1DLine1(self): + """ Basic test for Intersect2DMeshWith1DLine: a vertical line intersecting a square. """ + m1c = MEDCouplingCMesh() + coordX = DataArrayDouble([-1., 1., 2]) + m1c.setCoordsAt(0,coordX) + coordY = DataArrayDouble([0., 2.]) + m1c.setCoordsAt(1,coordY); + m1 = m1c.buildUnstructured() + + # A simple line: + m2 = MEDCouplingUMesh("bla", 1) + coord2 = DataArrayDouble([0.,-1.0, 0.,1., 0.,3., 0.5,2.2], 4, 2) + conn2 = DataArrayInt([NORM_SEG2,0,1,NORM_SEG3,1,2,3]) + connI2 = DataArrayInt([0,3,7]) + m2.setCoords(coord2) + m2.setConnectivity(conn2, connI2) + + # End of construction of input meshes m1bis and m2 -> start of specific part of the test + m3, map1, map2, mapI2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10) + m3.mergeNodes(1.0e-8) + + self.assertEqual(3,m3.getNumberOfCells()) + self.assertEqual(20,m3.getNumberOfNodes()) + self.assertEqual(2,m3.getSpaceDimension()) + # Mapping + exp1, exp2, expI2 = [0,0,1], [0,1,0,1], [0,2,4,4] + self.assertEqual(3, map1.getNumberOfTuples()) + self.assertEqual(4, map2.getNumberOfTuples()) + self.assertEqual(4, mapI2.getNumberOfTuples()) + self.assertEqual(exp1, map1.getValues()) + self.assertEqual(exp2, map2.getValues()) + self.assertEqual(expI2, mapI2.getValues()) + + expConn = [32, 1, 10, 7, 11, 4, 12, 13, 14, 15, 16, 32, 10, 0, 3, 11, 7, 17, 18, 19, 14, 13, 5, 2, 1, 4, 5] # 27tpl + expConnI = [0, 11, 22, 27] + self.assertEqual(27, m3.getNodalConnectivity().getNumberOfTuples()) + self.assertEqual(4, m3.getNodalConnectivityIndex().getNumberOfTuples()) + self.assertEqual(expConn, m3.getNodalConnectivity().getValues()) + self.assertEqual(expConnI, m3.getNodalConnectivityIndex().getValues()) + + def testSwig2Intersect2DMeshWith1DLine2(self): + """ Star pattern (a triangle intersecting another one upside down) """ + coords1 = DataArrayDouble([-2.,1., 2.,1., 0.,-2.], 3,2) + coords2 = DataArrayDouble([0.,2., 2.,-1., -2.,-1., 0.,3.], 4,2) + m1 = MEDCouplingUMesh("triangle", 2) + m2 = MEDCouplingUMesh("tri_line", 1) + m1.setCoords(coords1) + m2.setCoords(coords2) + m1.setConnectivity(DataArrayInt([NORM_TRI3, 0,1,2]), DataArrayInt([0,4])) + m2.setConnectivity(DataArrayInt([NORM_SEG2,0,1,NORM_SEG2,1,2,NORM_SEG2,2,3]), DataArrayInt([0,3,6,9])) + + # End of construction of input meshes m1bis and m2 -> start of specific part of the test + m3, map1, map2, mapI2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10) + # Mapping + exp1, exp2, expI2 = [0,0,0,0], [2, 0,1,2, 0, 1], [0,1,4,5,6] + self.assertEqual(exp1, map1.getValues()) + self.assertEqual(exp2, map2.getValues()) + self.assertEqual(expI2, mapI2.getValues()) + expConn = [5, 0, 8, 12, 5, 8, 7, 9, 10, 11, 12, 5, 7, 1, 9, 5, 10, 2, 11] + expConnI = [0, 4, 11, 15, 19] + self.assertEqual(expConn, m3.getNodalConnectivity().getValues()) + self.assertEqual(expConnI, m3.getNodalConnectivityIndex().getValues()) + + def testSwig2Intersect2DMeshWith1DLine3(self): + """ Line pieces ending (or fully located) in the middle of a cell """ + m1c = MEDCouplingCMesh() + m1c.setCoordsAt(0,DataArrayDouble([-1., 1.])) + m1c.setCoordsAt(1,DataArrayDouble([-1., 1.])); + m1 = m1c.buildUnstructured() + coords2 = DataArrayDouble([0.,0., 0.,1.5, -1.5,0., 0.6,1.1, 1.1,0.6, 0.5,0.0, 0.0,-0.5], 7,2) + m2 = MEDCouplingUMesh("piecewise_line", 1) + m2.setCoords(coords2) + c = DataArrayInt([NORM_SEG2,0,1, NORM_SEG2,1,2, NORM_SEG2,3,4, NORM_SEG2,5,6]) + cI = DataArrayInt([0,3,6,9,12]) + m2.setConnectivity(c, cI) + + # End of construction of input meshes m1bis and m2 -> start of specific part of the test + m3, map1, map2, mapI2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10) + self.assertEqual(3,m3.getNumberOfCells()) + # Mapping + exp1, exp2, expI2 = [0,0,0], [1,2, 1, 2], [0,2,3,4] + self.assertEqual(exp1, map1.getValues()) + self.assertEqual(exp2, map2.getValues()) + self.assertEqual(expI2, mapI2.getValues()) + expConn = [5, 1, 0, 11, 13, 12, 14, 15, 5, 11, 2, 13, 5, 14, 3, 15] + expConnI = [0, 8, 12, 16] + self.assertEqual(expConn, m3.getNodalConnectivity().getValues()) + self.assertEqual(expConnI, m3.getNodalConnectivityIndex().getValues()) + def setUp(self): pass pass if __name__ == '__main__': unittest.main() + diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 1f6e195b2..467474d25 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -2392,6 +2392,18 @@ namespace ParaMEDMEM return ret; } + static PyObject *Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps) throw(INTERP_KERNEL::Exception) + { + DataArrayInt *cellNb1=0,*cellNb2=0, *cellNbI2=0; + MEDCouplingUMesh *mret=MEDCouplingUMesh::Intersect2DMeshWith1DLine(m1,m2,eps,cellNb1,cellNb2, cellNbI2); + PyObject *ret=PyTuple_New(4); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(mret),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb1),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,3,SWIG_NewPointerObj(SWIG_as_voidptr(cellNbI2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + return ret; + } + PyObject *buildSlice3D(PyObject *origin, PyObject *vec, double eps) const throw(INTERP_KERNEL::Exception) { int spaceDim=self->getSpaceDimension();