X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling%2FMEDCouplingUMesh_intersection.cxx;h=f938bb103dd89d0f5b34d36151a7d98d6c0e6221;hb=0cc7b8ec634eb17df8b0b70c26e7473378e16a5b;hp=6f29181487dc2b3382be96e39c9d0d784ba96e45;hpb=305e7a33835941b70e0c8037351c92476ed4a08f;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx index 6f2918148..f938bb103 100644 --- a/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh_intersection.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D +// Copyright (C) 2007-2023 CEA, EDF // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -50,13 +50,13 @@ using namespace MEDCoupling; /// @cond INTERNAL -int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter) +mcIdType InternalAddPoint(const INTERP_KERNEL::Edge *e, mcIdType id, const double *coo, mcIdType startId, mcIdType endId, DataArrayDouble& addCoo, mcIdType& nodesCnter) { if(id!=-1) return id; else { - int ret(nodesCnter++); + mcIdType ret(nodesCnter++); double newPt[2]; e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt); addCoo.insertAtTheEnd(newPt,newPt+2); @@ -64,13 +64,13 @@ int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, in } } -int InternalAddPointOriented(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter) +mcIdType InternalAddPointOriented(const INTERP_KERNEL::Edge *e, mcIdType id, const double *coo, mcIdType startId, mcIdType endId, DataArrayDouble& addCoo, mcIdType& nodesCnter) { if(id!=-1) return id; else { - int ret(nodesCnter++); + mcIdType ret(nodesCnter++); double newPt[2]; e->getMiddleOfPointsOriented(coo+2*startId,coo+2*endId,newPt); addCoo.insertAtTheEnd(newPt,newPt+2); @@ -79,17 +79,17 @@ int InternalAddPointOriented(const INTERP_KERNEL::Edge *e, int id, const double } -void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) +void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const mcIdType *connBg, mcIdType offset, DataArrayIdType *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) { - int tmp[3]; + mcIdType tmp[3]; int trueStart(start>=0?start:nbOfEdges+start); - tmp[0]=linOrArc?(int)INTERP_KERNEL::NORM_QPOLYG:(int)INTERP_KERNEL::NORM_POLYGON; tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp]; + tmp[0]=ToIdType(linOrArc?INTERP_KERNEL::NORM_QPOLYG:INTERP_KERNEL::NORM_POLYGON); tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp]; newConnOfCell->insertAtTheEnd(tmp,tmp+3); if(linOrArc) { if(stp-start>1) { - int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); + mcIdType tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); InternalAddPointOriented(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2); middles.push_back(tmp3+offset); } @@ -98,15 +98,15 @@ void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int st } } -void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) +void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const mcIdType *connBg, mcIdType offset, DataArrayIdType *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) { - int tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]); + mcIdType tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]); newConnOfCell->pushBackSilent(tmpEnd); if(linOrArc) { if(stp-start>1) { - int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); + mcIdType tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2); middles.push_back(tmp3+offset); } @@ -115,15 +115,15 @@ void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int s } } -void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) +void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const mcIdType *connBg, mcIdType offset, DataArrayIdType *newConnOfCell, DataArrayDouble *appendedCoords, std::vector& middles) { // only the quadratic point to deal with: if(linOrArc) { if(stp-start>1) // if we are covering more than one segment we need to create a new mid point { - int tmpSrt(connBg[start]),tmpEnd(connBg[stp % nbOfEdges]); // % to handle last seg. - int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); + mcIdType tmpSrt(connBg[start]),tmpEnd(connBg[stp % nbOfEdges]); // % to handle last seg. + mcIdType tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2); InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2); middles.push_back(tmp3+offset); } @@ -132,20 +132,20 @@ void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, } } -void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map,int>& m, int forbVal0, int forbVal1, std::vector& isect) +void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map,mcIdType>& m, mcIdType forbVal0, mcIdType forbVal1, std::vector& isect) { MCAuto nTmp(n); nTmp->incrRef(); - std::map,int>::const_iterator it(m.find(nTmp)); + std::map,mcIdType>::const_iterator it(m.find(nTmp)); if(it==m.end()) throw INTERP_KERNEL::Exception("Internal error in remapping !"); - int v((*it).second); + mcIdType v((*it).second); if(v==forbVal0 || v==forbVal1) return ; if(std::find(isect.begin(),isect.end(),v)==isect.end()) isect.push_back(v); } -bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map,int>& m, int forbVal0, int forbVal1, std::vector& isect) +bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map,mcIdType>& m, mcIdType forbVal0, mcIdType forbVal1, std::vector& isect) { int sz(c.size()); if(sz<=1) @@ -166,7 +166,7 @@ bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map< namespace MEDCoupling { - INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge2(INTERP_KERNEL::NormalizedCellType typ, const int *bg, const double *coords2D, std::map< MCAuto,int>& m) + INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge2(INTERP_KERNEL::NormalizedCellType typ, const mcIdType *bg, const double *coords2D, std::map< MCAuto,mcIdType>& m) { INTERP_KERNEL::Edge *ret(0); MCAuto n0(new INTERP_KERNEL::Node(coords2D[2*bg[0]],coords2D[2*bg[0]+1])),n1(new INTERP_KERNEL::Node(coords2D[2*bg[1]],coords2D[2*bg[1]+1])); @@ -198,9 +198,13 @@ namespace MEDCoupling return ret; } - INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map >& mapp2, const int *bg) + INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map& mapp2, const mcIdType *bg) { INTERP_KERNEL::Edge *ret=0; + + mapp2[bg[0]].second = INTERP_KERNEL::USAGE_LINEAR; + mapp2[bg[1]].second = INTERP_KERNEL::USAGE_LINEAR; + switch(typ) { case INTERP_KERNEL::NORM_SEG2: @@ -220,7 +224,8 @@ namespace MEDCoupling ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first); else ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]].first,mapp2[bg[2]].first,mapp2[bg[1]].first); - mapp2[bg[2]].second=false; + if (mapp2[bg[2]].second != INTERP_KERNEL::USAGE_LINEAR) // switch the node usage to quadratic only if it is not used as an extreme point for another edge + mapp2[bg[2]].second = INTERP_KERNEL::USAGE_QUADRATIC_ONLY; break; } default: @@ -235,47 +240,98 @@ namespace MEDCoupling * The input mesh 'mDesc' must be so that mDim==1 and spaceDim==2. * 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'. */ - INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector& candidates, - std::map& mapp) + INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector& candidates, + std::map& mapp) { mapp.clear(); - std::map > mapp2;//bool is for a flag specifying if node is boundary (true) or only a middle for SEG3. + std::map mapp2; // the last var is a flag specifying if node is an extreme node of the seg (LINEAR) or only a middle for SEG3 (QUADRATIC_ONLY). const double *coo=mDesc->getCoords()->getConstPointer(); - const int *c=mDesc->getNodalConnectivity()->getConstPointer(); - const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer(); - std::set s; - for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + const mcIdType *c=mDesc->getNodalConnectivity()->getConstPointer(); + const mcIdType *cI=mDesc->getNodalConnectivityIndex()->getConstPointer(); + std::set s; + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) s.insert(c+cI[*it]+1,c+cI[(*it)+1]); - for(std::set::const_iterator it2=s.begin();it2!=s.end();it2++) + for(std::set::const_iterator it2=s.begin();it2!=s.end();it2++) { INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]); - mapp2[*it2]=std::pair(n,true); + mapp2[*it2]=INTERP_KERNEL::NodeWithUsage(n,INTERP_KERNEL::USAGE_UNKNOWN); } INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon; - for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) { INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]]; ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1)); } - for(std::map >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++) + for(std::map::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++) { - if((*it2).second.second) + if((*it2).second.second == INTERP_KERNEL::USAGE_LINEAR) mapp[(*it2).second.first]=(*it2).first; ((*it2).second.first)->decrRef(); } return ret; } - INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo) + INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMeshWithTree(const MEDCouplingUMesh *mDesc, const std::vector& candidates, + std::map& mapp, + const BBTreePts<2,mcIdType> & nodeTree, + const std::map& mapRev) + { + mapp.clear(); + std::map mapp2; // the last var is a flag specifying if node is an extreme node of the seg (LINEAR) or only a middle for SEG3 (QUADRATIC_ONLY). + const double *coo=mDesc->getCoords()->getConstPointer(); + const mcIdType *c=mDesc->getNodalConnectivity()->getConstPointer(); + const mcIdType *cI=mDesc->getNodalConnectivityIndex()->getConstPointer(); + std::set s; + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + s.insert(c+cI[*it]+1,c+cI[(*it)+1]); + for(std::set::const_iterator it2=s.begin();it2!=s.end();it2++) + { + INTERP_KERNEL::Node *n; + // Look for a potential node to merge + std::vector candNode; + nodeTree.getElementsAroundPoint(coo+2*(*it2), candNode); + if (candNode.size() > 2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MEDCouplingUMeshBuildQPFromMeshWithTree(): some nodes are not properly merged (within eps) in input mesh!"); + bool node_created=false; + if (candNode.size()) + { + auto itt=mapRev.find(candNode[0]); + if (itt != mapRev.end()) // we might hit a node which is in the coords array but not used in the connectivity in which case it won't be in the revMap + { + node_created=true; + n = (*itt).second; + n->incrRef(); + } + } + if(!node_created) + n = new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]); + mapp2[*it2]=INTERP_KERNEL::NodeWithUsage(n,INTERP_KERNEL::USAGE_UNKNOWN); + } + INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon; + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + { + INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]]; + ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1)); // this call will set quad points to false in the map + } + for(std::map::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++) + { + if((*it2).second.second == INTERP_KERNEL::USAGE_LINEAR) + mapp[(*it2).second.first]=(*it2).first; + ((*it2).second.first)->decrRef(); + } + return ret; + } + + INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(mcIdType nodeId, const double *coo1, mcIdType offset1, const double *coo2, mcIdType offset2, const std::vector& addCoo) { if(nodeId>=offset2) { - int locId=nodeId-offset2; + mcIdType locId=nodeId-offset2; return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]); } if(nodeId>=offset1) { - int locId=nodeId-offset1; + mcIdType locId=nodeId-offset1; return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]); } return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]); @@ -284,16 +340,16 @@ namespace MEDCoupling /** * 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, - /*output*/std::map& mapp, std::map& mappRev) + void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, mcIdType offset1, const double *coo2, mcIdType offset2, const std::vector& addCoo, + const mcIdType *desc1Bg, const mcIdType *desc1End, const std::vector >& intesctEdges1, + /*output*/std::map& mapp, std::map& mappRev) { - for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++) + for(const mcIdType *desc1=desc1Bg;desc1!=desc1End;desc1++) { - int eltId1=abs(*desc1)-1; - for(std::vector::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++) + mcIdType eltId1=std::abs(*desc1)-1; + for(std::vector::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++) { - std::map::const_iterator it=mappRev.find(*it1); + std::map::const_iterator it=mappRev.find(*it1); if(it==mappRev.end()) { INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo); @@ -310,26 +366,30 @@ namespace MEDCoupling /*! * Returns true if a colinearization has been found in the given cell. If false is returned the content pushed in \a newConnOfCell is equal to [ \a connBg , \a connEnd ) . * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair). + * \param forbiddenPoints the list of points that should not be removed in the process */ -bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords) +bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const mcIdType *connBg, const mcIdType *connEnd, mcIdType offset, + const std::map& forbiddenPoints, + DataArrayIdType *newConnOfCell, DataArrayDouble *appendedCoords) { std::size_t sz(std::distance(connBg,connEnd)); if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell. throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !"); sz--; - INTERP_KERNEL::AutoPtr tmpConn(new int[sz]); + INTERP_KERNEL::AutoPtr tmpConn(new mcIdType[sz]); + INTERP_KERNEL::AutoPtr tmpConn2(new mcIdType[sz]); const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0])); - unsigned nbs(cm.getNumberOfSons2(connBg+1,sz)); + unsigned nbs(cm.getNumberOfSons2(connBg+1,ToIdType(sz))); unsigned nbOfHit(0); // number of fusions operated int posBaseElt(0),posEndElt(0),nbOfTurn(0); - const unsigned int maxNbOfHit = cm.isQuadratic() ? nbs-2 : nbs-3; // a quad cell is authorized to end up with only two edges, a linear one has to keep 3 at least + const std::size_t maxNbOfHit = cm.isQuadratic() ? nbs-2 : nbs-3; // a quad cell is authorized to end up with only two edges, a linear one has to keep 3 at least INTERP_KERNEL::NormalizedCellType typeOfSon; - std::vector middles; + std::vector middles; bool ret(false); for(;(nbOfTurn+nbOfHit),int> m; + cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,ToIdType(sz),tmpConn,typeOfSon); + std::map,mcIdType> m; INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m)); posEndElt = posBaseElt+1; @@ -339,8 +399,13 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg { for(unsigned i=1;iareColinears(); if(isColinear) @@ -353,14 +418,21 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg eCand->decrRef(); if(!isColinear) break; + // Update last connectivity + std::copy((mcIdType *)tmpConn2, tmpConn2+sz, (mcIdType *)tmpConn); } } // Now move forward: const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt); // the first element to be inspected going forward for(unsigned j=fwdStart+1;jareColinears()); if(isColinear) @@ -373,18 +445,20 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg eCand->decrRef(); if(!isColinear) break; + // Update last connectivity + std::copy((mcIdType *)tmpConn2, tmpConn2+sz, (mcIdType *)tmpConn); } //push [posBaseElt,posEndElt) in newConnOfCell using e // The if clauses below are (voluntary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its beginning! if(nbOfTurn==0) // at the beginning of the connectivity (insert type) - EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); + EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); else if((nbOfHit+nbOfTurn) != (nbs-1)) // in the middle - EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); + EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); if ((nbOfHit+nbOfTurn) == (nbs-1)) // at the end (only quad points to deal with) - EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); + EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles); posBaseElt=posEndElt; e->decrRef(); } @@ -395,14 +469,14 @@ bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg -bool IsColinearOfACellOf(const std::vector< std::vector >& intersectEdge1, const std::vector& candidates, int start, int stop, int& retVal) +bool IsColinearOfACellOf(const std::vector< std::vector >& intersectEdge1, const std::vector& candidates, mcIdType start, mcIdType stop, mcIdType& retVal) { if(candidates.empty()) return false; - for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) { - const std::vector& pool(intersectEdge1[*it]); - int tmp[2]; tmp[0]=start; tmp[1]=stop; + const std::vector& pool(intersectEdge1[*it]); + mcIdType tmp[2]; tmp[0]=start; tmp[1]=stop; if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end()) { retVal=*it+1; @@ -434,37 +508,37 @@ bool IsColinearOfACellOf(const std::vector< std::vector >& intersectEdge1, */ void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, const std::vector& addCoo, - const std::vector< std::vector >& subDiv, std::vector< std::vector >& intersectEdge) + const std::vector< std::vector >& subDiv, std::vector< std::vector >& intersectEdge) { - int offset1=m1->getNumberOfNodes(); - int ncell=m2->getNumberOfCells(); - const int *c=m2->getNodalConnectivity()->begin(); - const int *cI=m2->getNodalConnectivityIndex()->begin(); + mcIdType offset1=m1->getNumberOfNodes(); + mcIdType ncell2=m2->getNumberOfCells(); + const mcIdType *c=m2->getNodalConnectivity()->begin(); + const mcIdType *cI=m2->getNodalConnectivityIndex()->begin(); const double *coo=m2->getCoords()->begin(); const double *cooBis=m1->getCoords()->begin(); - int offset2=offset1+m2->getNumberOfNodes(); - intersectEdge.resize(ncell); - for(int i=0;igetNumberOfNodes(); + intersectEdge.resize(ncell2); + for(mcIdType i=0;i& divs=subDiv[i]; - int nnode=cI[1]-cI[0]-1; - std::map > mapp2; - std::map mapp22; - for(int j=0;j& divs=subDiv[i]; + mcIdType nnode=cI[1]-cI[0]-1; + std::map mapp2; + std::map mapp22; + for(mcIdType j=0;j(nn,true); + mcIdType nnid=c[(*cI)+j+1]; + mapp2[nnid]=INTERP_KERNEL::NodeWithUsage(nn,INTERP_KERNEL::USAGE_UNKNOWN); mapp22[nn]=nnid+offset1; } INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1); - for(std::map >::const_iterator it=mapp2.begin();it!=mapp2.end();it++) + for(std::map::const_iterator it=mapp2.begin();it!=mapp2.end();it++) ((*it).second.first)->decrRef(); std::vector addNodes(divs.size()); - std::map mapp3; + std::map mapp3; for(std::size_t j=0;j >& intersectEdge2, const DataArrayDouble *coords1, const std::vector& addCoo, const std::map& mergedNodes, const std::vector< std::vector >& colinear2, const std::vector< std::vector >& intersectEdge1, - MCAuto& idsInRetColinear, MCAuto& idsInMesh1DForIdsInRetColinear) +/* + * Build the final 1D mesh resulting from the newly created points after intersection with the segments of the descending 2D mesh. + * @param[out] idsInRetColinear IDs of edges in the result (ret) that are colinears to one of the segment of the descending 2D mesh. Indexing scheme + * is the one of the ret 1D mesh. + * @param[out] idsInMesh1DForIdsInRetColinear same IDs as above in the descending 2D mesh + */ +MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector >& intersectEdge2, + const DataArrayDouble *coords1, const std::vector& addCoo, const std::map& mergedNodes, + const std::vector< std::vector >& colinear2, const std::vector< std::vector >& intersectEdge1, + MCAuto& idsInRetColinear, MCAuto& idsInMesh1DForIdsInRetColinear) { - idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1); - idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1); - int nCells(mesh1D->getNumberOfCells()); - if(nCells!=(int)intersectEdge2.size()) + idsInRetColinear=DataArrayIdType::New(); idsInRetColinear->alloc(0,1); + idsInMesh1DForIdsInRetColinear=DataArrayIdType::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1); + mcIdType nCells=mesh1D->getNumberOfCells(); + if(nCells!=ToIdType(intersectEdge2.size())) throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !"); const DataArrayDouble *coo2(mesh1D->getCoords()); - const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin()); + const mcIdType *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin()); const double *coo2Ptr(coo2->begin()); - int offset1(coords1->getNumberOfTuples()); - int offset2(offset1+coo2->getNumberOfTuples()); - int offset3(offset2+addCoo.size()/2); + mcIdType offset1(coords1->getNumberOfTuples()); + mcIdType offset2(offset1+coo2->getNumberOfTuples()); + mcIdType offset3(offset2+ToIdType(addCoo.size())/2); std::vector addCooQuad; - MCAuto cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0); - int tmp[4],cicnt(0),kk(0); - for(int i=0;i cOut(DataArrayIdType::New()),ciOut(DataArrayIdType::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0); + mcIdType tmp[4],cicnt(0),kk(0); + for(mcIdType i=0;i,int> m; + std::map,mcIdType> m; INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m)); - const std::vector& subEdges(intersectEdge2[i]); - int nbSubEdge(subEdges.size()/2); - for(int j=0;j& subEdges(intersectEdge2[i]); + mcIdType nbSubEdge=ToIdType(subEdges.size()/2); + for(mcIdType j=0;j n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)); + MCAuto n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)), + n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)); MCAuto e2(e->buildEdgeLyingOnMe(n1,n2)); INTERP_KERNEL::Edge *e2Ptr(e2); - std::map::const_iterator itm; + std::map::const_iterator itm; if(dynamic_cast(e2Ptr)) { tmp[0]=INTERP_KERNEL::NORM_SEG3; @@ -518,7 +601,7 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j]; itm=mergedNodes.find(subEdges[2*j+1]); tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1]; - tmp[3]=offset3+(int)addCooQuad.size()/2; + tmp[3]=offset3+ToIdType(addCooQuad.size()/2); double tmp2[2]; e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2); cicnt+=4; @@ -536,7 +619,7 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: cOut->insertAtTheEnd(tmp,tmp+3); ciOut->pushBackSilent(cicnt); } - int tmp00; + mcIdType tmp00; if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00)) { idsInRetColinear->pushBackSilent(kk); @@ -545,11 +628,12 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: } e->decrRef(); } + MCAuto ret(MEDCouplingUMesh::New(mesh1D->getName(),1)); ret->setConnectivity(cOut,ciOut,true); MCAuto arr3(DataArrayDouble::New()); - arr3->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2); - MCAuto arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,C_DEALLOC,(int)addCooQuad.size()/2,2); + arr3->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,addCoo.size()/2,2); + MCAuto arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,DeallocType::C_DEALLOC,addCooQuad.size()/2,2); std::vector coordss(4); coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4; MCAuto arr(DataArrayDouble::Aggregate(coordss)); @@ -557,12 +641,12 @@ MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std:: return ret.retn(); } -MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector >& intersectEdge1) +MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector >& intersectEdge1) { - std::vector allEdges; - for(const int *it2(descBg);it2!=descEnd;it2++) + std::vector allEdges; + for(const mcIdType *it2(descBg);it2!=descEnd;it2++) { - const std::vector& edge1(intersectEdge1[std::abs(*it2)-1]); + const std::vector& edge1(intersectEdge1[std::abs(*it2)-1]); if(*it2>0) allEdges.insert(allEdges.end(),edge1.begin(),edge1.end()); else @@ -575,31 +659,31 @@ MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const MCAuto ret(MEDCouplingUMesh::New("",2)); ret->setCoords(coords); ret->allocateCells(1); - std::vector connOut(nbOfEdgesOf2DCellSplit); + std::vector connOut(nbOfEdgesOf2DCellSplit); for(std::size_t kk=0;kkinsertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]); + ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,ToIdType(connOut.size()),&connOut[0]); return ret.retn(); } -MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector >& intersectEdge1) +MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, mcIdType cellIdInMesh2D, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector >& intersectEdge1) { - const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin()); + const mcIdType *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin()); const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]])); - std::size_t ii(0); + int ii(0); unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1)); if(sz!=std::distance(descBg,descEnd)) throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !"); - INTERP_KERNEL::AutoPtr tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]); - std::vector allEdges,centers; + INTERP_KERNEL::AutoPtr tmpPtr(new mcIdType[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]); + std::vector allEdges,centers; const double *coordsPtr(coords->begin()); MCAuto addCoo(DataArrayDouble::New()); addCoo->alloc(0,1); - int offset(coords->getNumberOfTuples()); - for(const int *it2(descBg);it2!=descEnd;it2++,ii++) + mcIdType offset(coords->getNumberOfTuples()); + for(const mcIdType *it2(descBg);it2!=descEnd;it2++,ii++) { INTERP_KERNEL::NormalizedCellType typeOfSon; cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon); - const std::vector& edge1(intersectEdge1[std::abs(*it2)-1]); + const std::vector& edge1(intersectEdge1[std::abs(*it2)-1]); if(*it2>0) allEdges.insert(allEdges.end(),edge1.begin(),edge1.end()); else @@ -608,11 +692,11 @@ MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, con centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center. else {//the current edge has been subsplit -> create corresponding centers. - std::size_t nbOfCentersToAppend(edge1.size()/2); - std::map< MCAuto,int> m; + mcIdType nbOfCentersToAppend=ToIdType(edge1.size()/2); + std::map< MCAuto,mcIdType> m; MCAuto ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m)); - std::vector::const_iterator it3(allEdges.end()-edge1.size()); - for(std::size_t k=0;k::const_iterator it3(allEdges.end()-edge1.size()); + for(mcIdType k=0;ksetCoords(addCoo); } ret->allocateCells(1); - std::vector connOut(nbOfEdgesOf2DCellSplit); + std::vector connOut(nbOfEdgesOf2DCellSplit); for(std::size_t kk=0;kkinsertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]); + ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,ToIdType(connOut.size()),&connOut[0]); return ret.retn(); } @@ -651,7 +735,7 @@ MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, con * * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()] */ -MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector >& intersectEdge1) +MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, mcIdType cellIdInMesh2D, const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector >& intersectEdge1) { const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D))); if(!cm.isQuadratic()) @@ -660,7 +744,7 @@ MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCou return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1); } -void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector& conn, const std::vector< MCAuto >& edges) +void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector& conn, const std::vector< MCAuto >& edges) { bool isQuad(false); for(std::vector< MCAuto >::const_iterator it=edges.begin();it!=edges.end();it++) @@ -670,25 +754,25 @@ void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector& conn, con isQuad=true; } if(!isQuad) - mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]); + mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,ToIdType(conn.size()),&conn[0]); else { const double *coo(mesh2D->getCoords()->begin()); std::size_t sz(conn.size()); std::vector addCoo; - std::vector conn2(conn); - int offset(mesh2D->getNumberOfNodes()); + std::vector conn2(conn); + mcIdType offset(mesh2D->getNumberOfNodes()); for(std::size_t i=0;igetMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i addCoo.insert(addCoo.end(),tmp,tmp+2); - conn2.push_back(offset+(int)i); + conn2.push_back(offset+ToIdType(i)); } mesh2D->getCoords()->rearrange(1); mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size()); mesh2D->getCoords()->rearrange(2); - mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]); + mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,ToIdType(conn2.size()),&conn2[0]); } } @@ -698,23 +782,23 @@ void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector& conn, con * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using * a set of edges defined in \a splitMesh1D. */ -void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector& edge1Bis, const std::vector< MCAuto >& edge1BisPtr, - std::vector< std::vector >& out0, std::vector< std::vector< MCAuto > >& out1) +void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector& edge1Bis, const std::vector< MCAuto >& edge1BisPtr, + std::vector< std::vector >& out0, std::vector< std::vector< MCAuto > >& out1) { std::size_t nb(edge1Bis.size()/2); std::size_t nbOfEdgesOf2DCellSplit(nb/2); - int iEnd(splitMesh1D->getNumberOfCells()); + mcIdType iEnd=splitMesh1D->getNumberOfCells(); if(iEnd==0) throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !"); std::size_t ii,jj; - const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin()); + const mcIdType *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin()); for(ii=0;ii single output cell out0.resize(1); out1.resize(1); - std::vector& connOut(out0[0]); + std::vector& connOut(out0[0]); connOut.resize(nbOfEdgesOf2DCellSplit); std::vector< MCAuto >& edgesPtr(out1[0]); edgesPtr.resize(nbOfEdgesOf2DCellSplit); @@ -728,25 +812,25 @@ void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vec { // [i,iEnd[ contains the out0.resize(2); out1.resize(2); - std::vector& connOutLeft(out0[0]); - std::vector& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1] + std::vector& connOutLeft(out0[0]); + std::vector& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1] std::vector< MCAuto >& eleft(out1[0]); std::vector< MCAuto >& eright(out1[1]); for(std::size_t k=ii;k > ees(iEnd); - for(int ik=0;ik,int> m; + std::map< MCAuto,mcIdType> m; MCAuto ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m)); ees[ik]=ee; } - for(int ik=iEnd-1;ik>=0;ik--) + for(mcIdType ik=iEnd-1;ik>=0;ik--) connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]); for(std::size_t k=jj+1;k& edges, const std::vector< MCAuto >& edgesPtr); + CellInfo(const std::vector& edges, const std::vector< MCAuto >& edgesPtr); public: - std::vector _edges; + std::vector _edges; std::vector< MCAuto > _edges_ptr; }; -CellInfo::CellInfo(const std::vector& edges, const std::vector< MCAuto >& edgesPtr) +CellInfo::CellInfo(const std::vector& edges, const std::vector< MCAuto >& edgesPtr) { std::size_t nbe(edges.size()); - std::vector edges2(2*nbe); std::vector< MCAuto > edgesPtr2(2*nbe); + std::vector edges2(2*nbe); std::vector< MCAuto > edgesPtr2(2*nbe); for(std::size_t i=0;i& edges, const std::vector< MCAuto& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { } - EdgeInfo(int istart, int iend, int pos, const MCAuto& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { } - bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; } - void somethingHappendAt(int pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight); - void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const; + EdgeInfo(mcIdType istart, mcIdType iend, const MCAuto& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { } + EdgeInfo(mcIdType istart, mcIdType iend, mcIdType pos, const MCAuto& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { } + bool isInMyRange(mcIdType pos) const { return pos>=_istart && pos<_iend; } + void somethingHappendAt(mcIdType pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight); + void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, mcIdType offset, mcIdType neighbors[2]) const; private: - int _istart; - int _iend; + mcIdType _istart; + mcIdType _iend; MCAuto _mesh; MCAuto _edge; - int _left; - int _right; + mcIdType _left; // index (local numbering) of the left 2D cell bordering the edge '_edge' + mcIdType _right; // same as above, right side. }; -void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight) +/* + * Update indices of left and right 2D cell bordering the current edge. + */ +void EdgeInfo::somethingHappendAt(mcIdType pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight) { const MEDCouplingUMesh *mesh(_mesh); if(mesh) @@ -802,6 +889,8 @@ void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAutopos) { _left++; _right++; return ; } + if (_right > pos && _left != pos) + { _right++; return ; } if(_right==pos) { bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end()); @@ -834,7 +923,7 @@ void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto barys(mesh->computeCellCenterOfMass()); - int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps)); + mcIdType cellId(mesh2D->getCellContainingPoint(barys->begin(),eps)); if(cellId==-1) throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !"); neighbors[0]=offset+cellId; neighbors[1]=offset+cellId; @@ -862,32 +951,32 @@ void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int of class VectorOfCellInfo { public: - VectorOfCellInfo(const std::vector& edges, const std::vector< MCAuto >& edgesPtr); + VectorOfCellInfo(const std::vector& edges, const std::vector< MCAuto >& edgesPtr); std::size_t size() const { return _pool.size(); } - int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const; - void setMeshAt(std::size_t pos, const MCAuto& mesh, int istart, int iend, const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MCAuto > >& edgePtrs); - const std::vector& getConnOf(int pos) const { return get(pos)._edges; } - const std::vector< MCAuto >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; } + mcIdType getPositionOf(double eps, const MEDCouplingUMesh *mesh) const; + void setMeshAt(mcIdType pos, const MCAuto& mesh, mcIdType istart, mcIdType iend, const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MCAuto > >& edgePtrs); + const std::vector& getConnOf(mcIdType pos) const { return get(pos)._edges; } + const std::vector< MCAuto >& getEdgePtrOf(mcIdType pos) const { return get(pos)._edges_ptr; } MCAuto getZeMesh() const { return _ze_mesh; } - void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const; + void feedEdgeInfoAt(double eps, mcIdType pos, mcIdType offset, mcIdType neighbors[2]) const; private: - int getZePosOfEdgeGivenItsGlobalId(int pos) const; - void updateEdgeInfo(int pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight); - const CellInfo& get(int pos) const; - CellInfo& get(int pos); + mcIdType getZePosOfEdgeGivenItsGlobalId(mcIdType pos) const; + void updateEdgeInfo(mcIdType pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight); + const CellInfo& get(mcIdType pos) const; + CellInfo& get(mcIdType pos); private: - std::vector _pool; - MCAuto _ze_mesh; - std::vector _edge_info; + std::vector _pool; // for a newly created 2D cell, the list of edges ToIdType( and edges ptr constiuing it + MCAuto _ze_mesh; // the aggregated mesh + std::vector _edge_info; // for each new edge added when cuting the 2D cell, the information on left and right bordering 2D cell }; -VectorOfCellInfo::VectorOfCellInfo(const std::vector& edges, const std::vector< MCAuto >& edgesPtr):_pool(1) +VectorOfCellInfo::VectorOfCellInfo(const std::vector& edges, const std::vector< MCAuto >& edgesPtr):_pool(1) { _pool[0]._edges=edges; _pool[0]._edges_ptr=edgesPtr; } -int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const +mcIdType VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const { if(_pool.empty()) throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !"); @@ -900,7 +989,9 @@ int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) co return zeMesh->getCellContainingPoint(barys->begin(),eps); } -void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto& mesh, int istart, int iend, const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MCAuto > >& edgePtrs) +void VectorOfCellInfo::setMeshAt(mcIdType pos, const MCAuto& mesh, mcIdType istart, mcIdType iend, + const MCAuto& mesh1DInCase, const std::vector< std::vector >& edges, + const std::vector< std::vector< MCAuto > >& edgePtrs) { get(pos);//to check pos bool isFast(pos==0 && _pool.size()==1); @@ -912,11 +1003,11 @@ void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back())); // std::vector pool(_pool.size()-1+sz); - for(std::size_t i=0;i _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2); } -void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const +void VectorOfCellInfo::feedEdgeInfoAt(double eps, mcIdType pos, mcIdType offset, mcIdType neighbors[2]) const { _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors); } -int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const +mcIdType VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(mcIdType pos) const { if(pos<0) throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !"); - int ret(0); + mcIdType ret(0); for(std::vector::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++) { if((*it).isInMyRange(pos)) @@ -965,9 +1056,9 @@ int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !"); } -void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight) +void VectorOfCellInfo::updateEdgeInfo(mcIdType pos, const std::vector< MCAuto >& newLeft, const std::vector< MCAuto >& newRight) { - get(pos);//to check; + get(pos);//to perform the sanity check; if(_edge_info.empty()) return ; std::size_t sz(_edge_info.size()-1); @@ -975,16 +1066,16 @@ void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto=(int)_pool.size()) + if(pos<0 || pos>=ToIdType(_pool.size())) throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !"); return _pool[pos]; } -CellInfo& VectorOfCellInfo::get(int pos) +CellInfo& VectorOfCellInfo::get(mcIdType pos) { - if(pos<0 || pos>=(int)_pool.size()) + if(pos<0 || pos>=ToIdType(_pool.size())) throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !"); return _pool[pos]; } @@ -1000,43 +1091,43 @@ CellInfo& VectorOfCellInfo::get(int pos) * * \param [in] allEdges a list of pairs (beginNode, endNode). Represents all edges (already cut) in the single 2D cell being handled here. Linked with \a allEdgesPtr to get the equation of edge. */ -MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector& allEdges, const std::vector< MCAuto >& allEdgesPtr, int offset, - MCAuto& idsLeftRight) +MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector& allEdges, const std::vector< MCAuto >& allEdgesPtr, mcIdType offset, + MCAuto& idsLeftRight) { - int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells()); + mcIdType nbCellsInSplitMesh1D=splitMesh1D->getNumberOfCells(); if(nbCellsInSplitMesh1D==0) throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !"); - const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin()); + const mcIdType *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin()); std::size_t nb(allEdges.size()),jj; if(nb%2!=0) throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !"); - std::vector edge1Bis(nb*2); + std::vector edge1Bis(nb*2); std::vector< MCAuto > edge1BisPtr(nb*2); std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()); std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb); std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()); std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb); // - idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2); - int *idsLeftRightPtr(idsLeftRight->getPointer()); + idsLeftRight=DataArrayIdType::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2); + mcIdType *idsLeftRightPtr(idsLeftRight->getPointer()); VectorOfCellInfo pool(edge1Bis,edge1BisPtr); // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start // of the connectivity. - MCAuto renumb(DataArrayInt::New()); + MCAuto renumb(DataArrayIdType::New()); renumb->alloc(nbCellsInSplitMesh1D,1); - const int * renumbP(renumb->begin()); + const mcIdType * renumbP(renumb->begin()); - int i, first=cSplitPtr[1]; + mcIdType i, first=cSplitPtr[1]; // Follow 1D line backward as long as it is connected: for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--) first=cSplitPtr[ciSplitPtr[i]+1]; if (i < nbCellsInSplitMesh1D-1) { // Build circular permutation to shift consecutive edges together - renumb->iota(i+1); + renumb->iota(nbCellsInSplitMesh1D-1-i); renumb->applyModulus(nbCellsInSplitMesh1D); splitMesh1D->renumberCells(renumbP, false); cSplitPtr = splitMesh1D->getNodalConnectivity()->begin(); @@ -1045,10 +1136,15 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh else renumb->iota(); // - - for(int iStart=0;iStart partOfSplitMesh1D(static_cast(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true))); - int pos(pool.getPositionOf(eps,partOfSplitMesh1D)); + mcIdType pos(pool.getPositionOf(eps,partOfSplitMesh1D)); // MCAutoretTmp(MEDCouplingUMesh::New("",2)); retTmp->setCoords(splitMesh1D->getCoords()); retTmp->allocateCells(); - std::vector< std::vector > out0; + std::vector< std::vector > out0; std::vector< std::vector< MCAuto > > out1; BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1); @@ -1077,7 +1173,7 @@ MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh // iStart=iEnd; } - for(int mm=0;mm >& intersectEdge1, int offset, - MCAuto& idsLeftRight) +MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, mcIdType cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D, + const mcIdType *descBg, const mcIdType *descEnd, const std::vector< std::vector >& intersectEdge1, mcIdType offset, + MCAuto& idsLeftRight) { - const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin()); + const mcIdType *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin()); // - std::vector allEdges; + std::vector allEdges; std::vector< MCAuto > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D - for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode + for(const mcIdType *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode { - int edgeId(std::abs(*it)-1); - std::map< MCAuto,int> m; + mcIdType edgeId(std::abs(*it)-1); + std::map< MCAuto,mcIdType> m; MCAuto ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m)); - const std::vector& edge1(intersectEdge1[edgeId]); + const std::vector& edge1(intersectEdge1[edgeId]); if(*it>0) allEdges.insert(allEdges.end(),edge1.begin(),edge1.end()); else @@ -1112,7 +1208,7 @@ MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCo return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight); } -bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps) +bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const mcIdType *conn1, const INTERP_KERNEL::CellModel& typ2, const mcIdType *conn2, double eps) { if(!typ1.isQuadratic() && !typ2.isQuadratic()) {//easy case comparison not @@ -1155,26 +1251,26 @@ bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, con * * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction. */ -int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps) +mcIdType FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const mcIdType *candidatesIn2DBg, const mcIdType *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, mcIdType cellIdInMesh1DSplitRelative, double eps) { if(candidatesIn2DEnd==candidatesIn2DBg) throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !"); const double *coo(mesh2DSplit->getCoords()->begin()); if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1) return *candidatesIn2DBg; - int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1); + mcIdType edgeId(std::abs(cellIdInMesh1DSplitRelative)-1); MCAuto cur1D(static_cast(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true))); if(cellIdInMesh1DSplitRelative<0) cur1D->changeOrientationOfCells(); - const int *c1D(cur1D->getNodalConnectivity()->begin()); + const mcIdType *c1D(cur1D->getNodalConnectivity()->begin()); const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0])); - for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++) + for(const mcIdType *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++) { MCAuto cur2D(static_cast(mesh2DSplit->buildPartOfMySelf(it,it+1,true))); - const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin()); + const mcIdType *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin()); const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]])); unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1)); - INTERP_KERNEL::AutoPtr tmpPtr(new int[ci[1]-ci[0]]); + INTERP_KERNEL::AutoPtr tmpPtr(new mcIdType[ci[1]-ci[0]]); for(unsigned it2=0;it2 >& intersectEdge1, std::vector< std::vector >& colinear2, std::vector< std::vector >& subDiv2, std::vector& addCoo, std::map& mergedNodes) + std::vector< std::vector >& intersectEdge1, std::vector< std::vector >& colinear2, std::vector< std::vector >& subDiv2, std::vector& addCoo, std::map& mergedNodes) { static const int SPACEDIM=2; INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); - INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); - const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin()); + const mcIdType *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin()); // Build BB tree of all edges in the tool mesh (second mesh) MCAuto bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps)); const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); - int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells()); + mcIdType nDescCell1=m1Desc->getNumberOfCells(),nDescCell2=m2Desc->getNumberOfCells(); intersectEdge1.resize(nDescCell1); colinear2.resize(nDescCell2); subDiv2.resize(nDescCell2); - BBTree myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps); + BBTree myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps); + BBTreePts treeNodes2(m2Desc->getCoords()->begin(),0,0,m2Desc->getCoords()->getNumberOfTuples(),eps); - std::vector candidates1(1); - int offset1(m1Desc->getNumberOfNodes()); - int offset2(offset1+m2Desc->getNumberOfNodes()); - for(int i=0;i candidates1(1); + mcIdType offset1(m1Desc->getNumberOfNodes()); + mcIdType offset2(offset1+m2Desc->getNumberOfNodes()); + for(mcIdType i=0;i candidates2; // edges of mesh2 candidate for intersection + std::vector 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; + std::map map1,map2; + std::map revMap2; // 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(m2Desc,candidates2,map2); + // Build revMap2 + for (auto& kv : map2) + revMap2[kv.second] = kv.first; candidates1[0]=i; - INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1); + // In the construction of pol1 we might reuse nodes from pol2, that we have identified as to be merged. + INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMeshWithTree(m1Desc,candidates1,map1,treeNodes2, revMap2); // 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; @@ -1255,23 +1356,23 @@ void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const M * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs(). */ void MEDCouplingUMesh::IntersectDescending2DMeshes(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< std::vector >& intersectEdge1, std::vector< std::vector >& colinear2, std::vector< std::vector >& subDiv2, + MEDCouplingUMesh *& m1Desc, DataArrayIdType *&desc1, DataArrayIdType *&descIndx1, DataArrayIdType *&revDesc1, DataArrayIdType *&revDescIndx1, std::vector& addCoo, - MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2) + MEDCouplingUMesh *& m2Desc, DataArrayIdType *&desc2, DataArrayIdType *&descIndx2, DataArrayIdType *&revDesc2, DataArrayIdType *&revDescIndx2) { // Build desc connectivity - desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New(); - desc2=DataArrayInt::New(); - descIndx2=DataArrayInt::New(); - revDesc2=DataArrayInt::New(); - revDescIndx2=DataArrayInt::New(); - MCAuto dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1); - MCAuto dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2); + desc1=DataArrayIdType::New(); descIndx1=DataArrayIdType::New(); revDesc1=DataArrayIdType::New(); revDescIndx1=DataArrayIdType::New(); + desc2=DataArrayIdType::New(); + descIndx2=DataArrayIdType::New(); + revDesc2=DataArrayIdType::New(); + revDescIndx2=DataArrayIdType::New(); + MCAuto dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1); + MCAuto dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2); m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1); m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2); MCAuto dd9(m1Desc),dd10(m2Desc); - std::map notUsedMap; + std::map notUsedMap; Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap); m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef(); m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef(); @@ -1281,42 +1382,44 @@ void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, c * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the * (newly created) nodes corresponding to the edge intersections. * Output params: - * @param[out] cr, crI connectivity of the resulting mesh - * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2 + * @param[out] cr connectivity of the resulting mesh + * @param[out] crI connectivity of the resulting mesh + * @param[out] cNb1 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2 + * @param[out] cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2 * TODO: describe input parameters */ -void MEDCouplingUMesh::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, +void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const mcIdType *desc1, const mcIdType *descIndx1, + const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, + const MEDCouplingUMesh *m2, const mcIdType *desc2, const mcIdType *descIndx2, 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& addCoordsQuadratic, std::vector& cr, std::vector& crI, std::vector& cNb1, std::vector& cNb2) { static const int SPACEDIM=2; const double *coo1(m1->getCoords()->begin()); - const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin()); - int offset1(m1->getNumberOfNodes()); + const mcIdType *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin()); + mcIdType offset1(m1->getNumberOfNodes()); const double *coo2(m2->getCoords()->begin()); - const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin()); - int offset2(offset1+m2->getNumberOfNodes()); - int offset3(offset2+((int)addCoords.size())/2); + const mcIdType *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin()); + mcIdType offset2(offset1+m2->getNumberOfNodes()); + mcIdType offset3(offset2+ToIdType(addCoords.size())/2); MCAuto bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps)); const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin()); // Here a BBTree on 2D-cells, not on segments: - BBTree myTree(bbox2,0,0,m2->getNumberOfCells(),eps); - int ncell1(m1->getNumberOfCells()); + BBTree myTree(bbox2,0,0,m2->getNumberOfCells(),eps); + mcIdType ncell1=m1->getNumberOfCells(); crI.push_back(0); - for(int i=0;i candidates2; + std::vector candidates2; myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2); - std::map mapp; - std::map mappRev; + std::map mapp; + std::map mappRev; 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 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); // @@ -1326,10 +1429,12 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo for(it1.first();!it1.finished();it1.next()) edges1.insert(it1.current()->getPtr()); // - std::map > edgesIn2ForShare; // common edges + std::map > edgesIn2ForShare; // common edges std::vector pol2s(candidates2.size()); - int ii=0; - for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) + mcIdType ii=0; + // Build, for each intersecting cell candidate from mesh2, the corresponding QP. + // Again all the additional intersecting nodes are there. + for(std::vector::const_iterator 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); @@ -1339,8 +1444,14 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo 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); } + // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges: + for (auto &p: pol2s) + p.cleanDegeneratedConsecutiveEdges(); + edgesIn2ForShare.clear(); // removing temptation to use it further since it might now contain invalid edges. + /// ii=0; - for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) + // Now rebuild intersected cells from all this: + for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++) { INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]); pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2); @@ -1361,25 +1472,25 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo throw INTERP_KERNEL::Exception(oss.str()); } } - for(std::map::const_iterator it=mappRev.begin();it!=mappRev.end();it++) + for(std::map::const_iterator it=mappRev.begin();it!=mappRev.end();it++) (*it).second->decrRef(); } } -void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector& conn, const double *coords, double eps) +void InsertNodeInConnIfNecessary(mcIdType nodeIdToInsert, std::vector& conn, const double *coords, double eps) { - std::vector::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert)); + std::vector::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert)); if(it!=conn.end()) return ; std::size_t sz(conn.size()); std::size_t found(std::numeric_limits::max()); for(std::size_t i=0;i(),1./normm)); - std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies(),1./normm)); + std::transform(v1,v1+3,v1,std::bind(std::multiplies(),std::placeholders::_1,1./normm)); + std::transform(v2,v2+3,v2,std::bind(std::multiplies(),std::placeholders::_1,1./normm)); double v3[3]; v3[0]=v1[1]*v2[2]-v1[2]*v2[1]; v3[1]=v1[2]*v2[0]-v1[0]*v2[2]; v3[2]=v1[0]*v2[1]-v1[1]*v2[0]; double normm2(sqrt(v3[0]*v3[0]+v3[1]*v3[1]+v3[2]*v3[2])),dotTest(v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]); @@ -1395,13 +1506,13 @@ void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector& conn, con conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert); } -void SplitIntoToPart(const std::vector& conn, int pt0, int pt1, std::vector& part0, std::vector& part1) +void SplitIntoToPart(const std::vector& conn, mcIdType pt0, mcIdType pt1, std::vector& part0, std::vector& part1) { std::size_t sz(conn.size()); - std::vector *curPart(&part0); + std::vector *curPart(&part0); for(std::size_t i=0;i& conn, int pt0, int pt1, std::vector /*! * this method method splits cur cells 3D Surf in sub cells 3DSurf using the previous subsplit. This method is the last one used to clip. */ -void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair >& cut3DSurf, - const int *desc, const int *descIndx, const double *coords, double eps, - std::vector >& res) const +void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair >& cut3DSurf, + const mcIdType *desc, const mcIdType *descIndx, const double *coords, double eps, + std::vector >& res) const { checkFullyDefined(); if(getMeshDimension()!=3 || getSpaceDimension()!=3) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!"); - const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin()); - int nbOfCells(getNumberOfCells()); + const mcIdType *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin()); + mcIdType nbOfCells=getNumberOfCells(); if(nbOfCells!=1) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !"); - for(int i=0;i& p=cut3DSurf[desc[offset+j]]; + const std::pair& p=cut3DSurf[desc[offset+j]]; const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]])); - int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1; - INTERP_KERNEL::AutoPtr tmp(new int[sz]); + mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1; + INTERP_KERNEL::AutoPtr tmp(new mcIdType[sz]); INTERP_KERNEL::NormalizedCellType cmsId; unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId)); - std::vector elt((int *)tmp,(int *)tmp+nbOfNodesSon); + std::vector elt((mcIdType *)tmp,(mcIdType *)tmp+nbOfNodesSon); if(p.first!=-1 && p.second!=-1) { if(p.first!=-2) { InsertNodeInConnIfNecessary(p.first,elt,coords,eps); InsertNodeInConnIfNecessary(p.second,elt,coords,eps); - std::vector elt1,elt2; + std::vector elt1,elt2; SplitIntoToPart(elt,p.first,p.second,elt1,elt2); res.push_back(elt1); res.push_back(elt2); @@ -1465,22 +1576,23 @@ void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pairgetNumberOfTuples()); - MCAuto c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach); - const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin()); - int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer()); - int prevPosOfCi(ciPtr[0]); - for(int i=0;igetNumberOfTuples()); + MCAuto c(DataArrayIdType::New()); c->alloc((std::size_t)lgthToReach); + const mcIdType *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin()); + mcIdType *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer()); + mcIdType prevPosOfCi(ciPtr[0]); + for(mcIdType i=0;igetNumberOfTuples()),nodesCnt(getNumberOfNodes()); - MCAuto c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach); + mcIdType ncells=getNumberOfCells(); + mcIdType lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()); + mcIdType nodesCnt(getNumberOfNodes()); + MCAuto c(DataArrayIdType::New()); c->alloc((std::size_t)lgthToReach); MCAuto addCoo(DataArrayDouble::New()); addCoo->alloc(0,1); - const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin()); - const int *midPtr(mid->begin()),*midIPtr(midI->begin()); + const mcIdType *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin()); + const mcIdType *midPtr(mid->begin()),*midIPtr(midI->begin()); const double *oldCoordsPtr(getCoords()->begin()); - int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer()); - int prevPosOfCi(ciPtr[0]); - for(int i=0;igetPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer()); + mcIdType prevPosOfCi(ciPtr[0]); + for(mcIdType i=0;i e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns)); - for(int k=0;kcheckFullyDefined(); m2->checkFullyDefined(); INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); - INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!"); // Step 1: compute all edge intersections (new nodes) - std::vector< std::vector > intersectEdge1, colinear2, subDiv2; + std::vector< std::vector > intersectEdge1, colinear2, subDiv2; MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes - DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0; + DataArrayIdType *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0; std::vector addCoo,addCoordsQuadratic; // coordinates of newly created nodes IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2, m1Desc,desc1,descIndx1,revDesc1,revDescIndx1, addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2); revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef(); - MCAuto dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2); + MCAuto dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2); MCAuto dd5(m1Desc),dd6(m2Desc); // Step 2: re-order newly created nodes according to the ordering found in m2 - std::vector< std::vector > intersectEdge2; + std::vector< std::vector > intersectEdge2; BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2); subDiv2.clear(); dd5=0; dd6=0; // Step 3: - std::vector cr,crI; //no DataArrayInt because interface with Geometric2D - std::vector cNb1,cNb2; //no DataArrayInt because interface with Geometric2D + std::vector cr,crI; //no DataArrayIdType because interface with Geometric2D + std::vector cNb1,cNb2; //no DataArrayIdType because interface with Geometric2D BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo, /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2); // Step 4: Prepare final result: MCAuto addCooDa(DataArrayDouble::New()); - addCooDa->alloc((int)(addCoo.size())/2,2); + addCooDa->alloc(addCoo.size()/2,2); std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer()); MCAuto addCoordsQuadraticDa(DataArrayDouble::New()); - addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2); + addCoordsQuadraticDa->alloc(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; MCAuto coo(DataArrayDouble::Aggregate(coordss)); MCAuto ret(MEDCouplingUMesh::New("Intersect2D",2)); - MCAuto conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer()); - MCAuto connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer()); - MCAuto c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer()); - MCAuto c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer()); + MCAuto conn(DataArrayIdType::New()); conn->alloc(cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer()); + MCAuto connI(DataArrayIdType::New()); connI->alloc(crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer()); + MCAuto c1(DataArrayIdType::New()); c1->alloc(cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer()); + MCAuto c2(DataArrayIdType::New()); c2->alloc(cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer()); ret->setConnectivity(conn,connI,true); ret->setCoords(coo); cellNb1=c1.retn(); cellNb2=c2.retn(); @@ -1652,6 +1768,10 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 * and finally, in case of quadratic polygon the centers of edges new nodes. * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input. * + * \b WARNING: the 2D mesh should be correctly oriented for this method to work properly. Methods changeSpaceDimension() and + * orientCorrectly2DCells() can be used for this. + * \b WARNING: the two meshes should be "clean" (no un-merged nodes, no non-conformal cells) + * * \param [in] mesh2D - the 2D mesh (spacedim=meshdim=2) to be intersected using \a mesh1D tool. The mesh must be so that each point in the space covered by \a mesh2D * must be covered exactly by one entity, \b no \b more. If it is not the case, some tools are available to heal the mesh (conformize2D, mergeNodes) * \param [in] mesh1D - the 1D mesh (spacedim=2 meshdim=1) the is the tool that will be used to intersect \a mesh2D. \a mesh1D must be ordered consecutively. If it is not the case @@ -1667,7 +1787,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 * * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes */ -void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D) +void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayIdType *&cellIdInMesh2D, DataArrayIdType *&cellIdInMesh1D) { if(!mesh2D || !mesh1D) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !"); @@ -1677,23 +1797,22 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !"); // Step 1: compute all edge intersections (new nodes) - std::vector< std::vector > intersectEdge1, colinear2, subDiv2; + std::vector< std::vector > intersectEdge1, colinear2, subDiv2; std::vector addCoo,addCoordsQuadratic; // coordinates of newly created nodes INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); - INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); // // Build desc connectivity - DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New()); - MCAuto dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1); + DataArrayIdType *desc1(DataArrayIdType::New()),*descIndx1(DataArrayIdType::New()),*revDesc1(DataArrayIdType::New()),*revDescIndx1(DataArrayIdType::New()); + MCAuto dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1); MCAuto m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1)); - std::map mergedNodes; + std::map mergedNodes; Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes); // use mergeNodes to fix intersectEdge1 - for(std::vector< std::vector >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++) + for(std::vector< std::vector >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++) { std::size_t n((*it0).size()/2); - int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]); - std::map::const_iterator it1; + mcIdType eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]); + std::map::const_iterator it1; it1=mergedNodes.find(eltStart); if(it1!=mergedNodes.end()) (*it0)[0]=(*it1).second; @@ -1703,46 +1822,92 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, } // MCAuto addCooDa(DataArrayDouble::New()); - addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2); + addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,addCoo.size()/2,2); // Step 2: re-order newly created nodes according to the ordering found in m2 - std::vector< std::vector > intersectEdge2; + std::vector< std::vector > intersectEdge2; BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2); subDiv2.clear(); // Step 3: compute splitMesh1D - MCAuto idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear; - MCAuto ret2(DataArrayInt::New()); ret2->alloc(0,1); + MCAuto idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear; + MCAuto ret2(DataArrayIdType::New()); ret2->alloc(0,1); MCAuto ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1, - idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear)); - MCAuto ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits::max()); ret3->rearrange(2); - MCAuto idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells())); + idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear)); + + // ### Colinearity fix : + // if a node in ret1 has been merged with an already existing node in mesh2D, + // we might end up with edges in ret1 that are colinear with some edges in mesh2D + // even if none of the edges of the two original meshes were colinear. + // this procedure detects such edges and adds them in idsInRet1Colinear/idsInDescMesh2DForIdsInRetColinear + // a- find edges in ret1 that are in contact with 2 cells + MCAuto centerOfMassRet1(ret1->computeCellCenterOfMass()); + MCAuto cells,cellsIndex; + mesh2D->getCellsContainingPoints(centerOfMassRet1->begin(),centerOfMassRet1->getNumberOfTuples(),eps,cells,cellsIndex); + MCAuto cellsIndex2(DataArrayIdType::New()); cellsIndex2->alloc(0,1); + if (cellsIndex->getNumberOfTuples() > 1) + cellsIndex2 = cellsIndex->deltaShiftIndex(); + MCAuto idsInRet1With2Contacts(cellsIndex2->findIdsEqual(2)); + + MCAuto realIdsInDesc2D(desc1->deepCopy()); + realIdsInDesc2D->abs(); realIdsInDesc2D->applyLin(1,-1); + const mcIdType *cRet1(ret1->getNodalConnectivity()->begin()),*ciRet1(ret1->getNodalConnectivityIndex()->begin()); + for(const mcIdType *it=idsInRet1With2Contacts->begin();it!=idsInRet1With2Contacts->end();it++) + { + // b- find the edge that the 2 cells in m1Desc have in common: + // this is the edge which is colinear with the one in ret1 + const mcIdType* cellId1 = cells->begin() + cellsIndex->begin()[*it]; + const mcIdType* cellId2 = cells->begin() + cellsIndex->begin()[*it]+1; + + std::set s1(realIdsInDesc2D->begin()+dd2->begin()[*cellId1], realIdsInDesc2D->begin()+dd2->begin()[*cellId1+1]); + std::set s2(realIdsInDesc2D->begin()+dd2->begin()[*cellId2], realIdsInDesc2D->begin()+dd2->begin()[*cellId2+1]); + + std::vector commonEdgeId; + std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(), std::back_inserter(commonEdgeId)); + + // c- find correct orientation for commonEdgeId + const mcIdType* firstNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+1; + const mcIdType* secondNodeInColinearEdgeRet1 = cRet1 + ciRet1[*it]+2; + mcIdType commonEdgeIdCorrectlyOriented; + if(IsColinearOfACellOf(intersectEdge1, commonEdgeId, *firstNodeInColinearEdgeRet1,*secondNodeInColinearEdgeRet1, commonEdgeIdCorrectlyOriented)) + { + idsInRet1Colinear->pushBackSilent(*it); + idsInDescMesh2DForIdsInRetColinear->pushBackSilent(commonEdgeIdCorrectlyOriented); + } + } + // ### End colinearity fix + + MCAuto ret3(DataArrayIdType::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits::max()); ret3->rearrange(2); + MCAuto idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells())); // deal with cells in mesh2D that are not cut but only some of their edges are - MCAuto idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy()); + MCAuto idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy()); idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1); idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique(); - MCAuto out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of the edges of those cells + + MCAuto out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of the edges of those cells if(!idsInDesc2DToBeRefined->empty()) { - DataArrayInt *out0(0),*outi0(0); - MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0); - MCAuto outi0s(outi0); + DataArrayIdType *out0(0),*outi0(0); + DataArrayIdType::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0); + MCAuto outi0s(outi0); out0s=out0; out0s=out0s->buildUnique(); out0s->sort(true); } // MCAuto ret1NonCol(static_cast(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end()))); - MCAuto baryRet1(ret1NonCol->computeCellCenterOfMass()); - MCAuto elts,eltsIndex; - mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex); - MCAuto eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1); + MCAuto baryRet1(centerOfMassRet1->selectByTupleId(idsInRet1NotColinear->begin(), idsInRet1NotColinear->end())); + DataArrayIdType *out(0),*outi(0); + DataArrayIdType::ExtractFromIndexedArrays(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end(),cells,cellsIndex,out,outi); + MCAuto elts(out),eltsIndex(outi); + + MCAuto eltsIndex2(DataArrayIdType::New()); eltsIndex2->alloc(0,1); if (eltsIndex->getNumberOfTuples() > 1) eltsIndex2 = eltsIndex->deltaShiftIndex(); - MCAuto eltsIndex3(eltsIndex2->findIdsEqual(1)); + MCAuto eltsIndex3(eltsIndex2->findIdsEqual(1)); if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells()) throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !"); - MCAuto cellsToBeModified(elts->buildUnique()); - MCAuto untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells())); - if((DataArrayInt *)out0s) + MCAuto cellsToBeModified(elts->buildUnique()); + MCAuto untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells())); + if((DataArrayIdType *)out0s) untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one std::vector< MCAuto > outMesh2DSplit; // OK all is ready to insert in ret2 mesh @@ -1752,26 +1917,26 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, outMesh2DSplit.back()->setCoords(ret1->getCoords()); ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end()); } - if((DataArrayInt *)out0s) + if((DataArrayIdType *)out0s) {// here dealing with cells in out0s but not in cellsToBeModified - MCAuto fewModifiedCells(out0s->buildSubstraction(cellsToBeModified)); - const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin()); - for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++) + MCAuto fewModifiedCells(out0s->buildSubstraction(cellsToBeModified)); + const mcIdType *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin()); + for(const mcIdType *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++) { outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1)); ret1->setCoords(outMesh2DSplit.back()->getCoords()); } - int offset(ret2->getNumberOfTuples()); + mcIdType offset(ret2->getNumberOfTuples()); ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end()); - MCAuto partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1); - partOfRet3->fillWithValue(std::numeric_limits::max()); partOfRet3->rearrange(2); - int kk(0),*ret3ptr(partOfRet3->getPointer()); - for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++) + MCAuto partOfRet3(DataArrayIdType::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1); + partOfRet3->fillWithValue(std::numeric_limits::max()); partOfRet3->rearrange(2); + mcIdType kk(0),*ret3ptr(partOfRet3->getPointer()); + for(const mcIdType *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++) { - int faceId(std::abs(*it)-1); - for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++) + mcIdType faceId(std::abs(*it)-1); + for(const mcIdType *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++) { - int tmp(fewModifiedCells->findIdFirstEqual(*it2)); + mcIdType tmp(fewModifiedCells->findIdFirstEqual(*it2)); if(tmp!=-1) { if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1]) @@ -1799,17 +1964,17 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, } } cellsToBeModified=cellsToBeModified->buildUniqueNotSorted(); - for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++) + for(const mcIdType *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++) { - MCAuto idsNonColPerCell(elts->findIdsEqual(*it)); + MCAuto idsNonColPerCell(elts->findIdsEqual(*it)); idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end()); - MCAuto idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end())); + MCAuto idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end())); MCAuto partOfMesh1CuttingCur2DCell(static_cast(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end()))); - MCAuto partOfRet3; + MCAuto partOfRet3; MCAuto splitOfOneCell(BuildMesh2DCutFrom(eps,*it,m1Desc,partOfMesh1CuttingCur2DCell,dd1->begin()+dd2->getIJ(*it,0),dd1->begin()+dd2->getIJ((*it)+1,0),intersectEdge1,ret2->getNumberOfTuples(),partOfRet3)); ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true); outMesh2DSplit.push_back(splitOfOneCell); - for(std::size_t i=0;igetNumberOfCells();i++) + for(mcIdType i=0;igetNumberOfCells();i++) ret2->pushBackSilent(*it); } // @@ -1820,16 +1985,16 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, // ret1->getCoords()->setInfoOnComponents(compNames); MCAuto ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp)); - // To finish - filter ret3 - std::numeric_limits::max() -> -1 - negate values must be resolved. + // To finish - filter ret3 - std::numeric_limits::max() -> -1 - negate values must be resolved. ret3->rearrange(1); - MCAuto edgesToDealWith(ret3->findIdsStrictlyNegative()); - for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++) + MCAuto edgesToDealWith(ret3->findIdsStrictlyNegative()); + for(const mcIdType *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++) { - int old2DCellId(-ret3->getIJ(*it,0)-1); - MCAuto candidates(ret2->findIdsEqual(old2DCellId)); + mcIdType old2DCellId(-ret3->getIJ(*it,0)-1); + MCAuto candidates(ret2->findIdsEqual(old2DCellId)); ret3->setIJ(*it,0,FindRightCandidateAmong(ret2D,candidates->begin(),candidates->end(),ret1,*it%2==0?-((*it)/2+1):(*it)/2+1,eps));// div by 2 because 2 components natively in ret3 } - ret3->changeValue(std::numeric_limits::max(),-1); + ret3->changeValue(std::numeric_limits::max(),-1); ret3->rearrange(2); // splitMesh1D=ret1.retn(); @@ -1841,7 +2006,7 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, /*! * \b WARNING this method is \b potentially \b non \b const (if returned array is empty). * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) ! - * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this this method + * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this method * will suppress such edges to use sub edges in \a this. So this method does not add nodes in \a this if merged edges are both linear (INTERP_KERNEL::NORM_SEG2). * In the other cases new nodes can be created. If any are created, they will be appended at the end of the coordinates object before the invocation of this method. * @@ -1853,7 +2018,7 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method. * * \param [in] eps the relative error to detect merged edges. - * \return DataArrayInt * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array + * \return DataArrayIdType * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array * that the user is expected to deal with. * * \throw If \a this is not coherent. @@ -1861,31 +2026,30 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, * \throw If \a this has not meshDim equal to 2. * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells */ -DataArrayInt *MEDCouplingUMesh::conformize2D(double eps) +DataArrayIdType *MEDCouplingUMesh::conformize2D(double eps) { static const int SPACEDIM=2; checkConsistencyLight(); if(getSpaceDimension()!=2 || getMeshDimension()!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !"); - MCAuto desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New()); + MCAuto desc1(DataArrayIdType::New()),descIndx1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescIndx1(DataArrayIdType::New()); MCAuto mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1)); - const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin()); + const mcIdType *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin()); MCAuto bboxArr(mDesc->getBoundingBoxForBBTree(eps)); const double *bbox(bboxArr->begin()),*coords(getCoords()->begin()); - int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells()); - std::vector< std::vector > intersectEdge(nDescCell),overlapEdge(nDescCell); + mcIdType nCell=getNumberOfCells(),nDescCell=mDesc->getNumberOfCells(); + std::vector< std::vector > intersectEdge(nDescCell),overlapEdge(nDescCell); std::vector addCoo; - BBTree myTree(bbox,0,0,nDescCell,-eps); + BBTree myTree(bbox,0,0,nDescCell,-eps); INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); - INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); - for(int i=0;i candidates; + std::vector candidates; myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates); - for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) - if(*it>i) + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + if(*it>i) // we're dealing with pair of edges, no need to treat the same pair twice { - std::map,int> m; + std::map,mcIdType> m; INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)), *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m)); INTERP_KERNEL::MergePoints merge; @@ -1899,24 +2063,24 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps) } } // splitting done. sort intersect point in intersectEdge. - std::vector< std::vector > middle(nDescCell); - int nbOf2DCellsToBeSplit(0); + std::vector< std::vector > middle(nDescCell); + mcIdType nbOf2DCellsToBeSplit(0); bool middleNeedsToBeUsed(false); std::vector cells2DToTreat(nDescCell,false); - for(int i=0;i& isect(intersectEdge[i]); - int sz((int)isect.size()); + std::vector& isect(intersectEdge[i]); + std::size_t sz(isect.size()); if(sz>1) { - std::map,int> m; + std::map,mcIdType> m; INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)); e->sortSubNodesAbs(coords,isect); e->decrRef(); } if(sz!=0) { - int idx0(rdi[i]),idx1(rdi[i+1]); + mcIdType idx0(rdi[i]),idx1(rdi[i+1]); if(idx1-idx0!=1) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !"); if(!cells2DToTreat[rd[idx0]]) @@ -1925,22 +2089,22 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps) nbOf2DCellsToBeSplit++; } // try to reuse at most eventual 'middle' of SEG3 - std::vector& mid(middle[i]); + std::vector& mid(middle[i]); mid.resize(sz+1,-1); if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3) { middleNeedsToBeUsed=true; - const std::vector& candidates(overlapEdge[i]); - std::vector trueCandidates; - for(std::vector::const_iterator itc=candidates.begin();itc!=candidates.end();itc++) + const std::vector& candidates(overlapEdge[i]); + std::vector trueCandidates; + for(std::vector::const_iterator itc=candidates.begin();itc!=candidates.end();itc++) if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3) trueCandidates.push_back(*itc); - int stNode(c[ci[i]+1]),endNode(isect[0]); - for(int j=0;j::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++) + for(std::vector::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++) { - int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]); + mcIdType tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]); if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode)) { mid[j]=*itc; break; } } @@ -1950,29 +2114,29 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps) } } } - MCAuto ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1); + MCAuto ret(DataArrayIdType::New()),notRet(DataArrayIdType::New()); ret->alloc(nbOf2DCellsToBeSplit,1); if(nbOf2DCellsToBeSplit==0) return ret.retn(); // - int *retPtr(ret->getPointer()); - for(int i=0;igetPointer()); + for(mcIdType i=0;i mSafe,nSafe,oSafe,pSafe,qSafe,rSafe; - DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0); - MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n; - DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p; + MCAuto mSafe,nSafe,oSafe,pSafe,qSafe,rSafe; + DataArrayIdType *m(0),*n(0),*o(0),*p(0),*q(0),*r(0); + DataArrayIdType::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n; + DataArrayIdType::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p; if(middleNeedsToBeUsed) - { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; } + { DataArrayIdType::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; } MCAuto modif(static_cast(buildPartOfMySelf(ret->begin(),ret->end(),true))); - int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe)); + mcIdType nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe)); setCoords(modif->getCoords());//if nbOfNodesCreated==0 modif and this have the same coordinates pointer so this line has no effect. But for quadratic cases this line is important. setPartOfMySelf(ret->begin(),ret->end(),*modif); { - bool areNodesMerged; int newNbOfNodes; + bool areNodesMerged; mcIdType newNbOfNodes; if(nbOfNodesCreated!=0) - MCAuto tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes)); + MCAuto tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes)); } return ret.retn(); } @@ -1990,7 +2154,7 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps) * * If \a this is constituted by only linear 2D cells, this method is close to the computation of the convex hull of each cells in \a this. * - * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized. + * \return DataArrayIdType * - The list of cellIds in \a this that have at least one edge colinearized. * * \throw If \a this is not coherent. * \throw If \a this has not spaceDim equal to 2. @@ -1998,23 +2162,58 @@ DataArrayInt *MEDCouplingUMesh::conformize2D(double eps) * * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D. */ -DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps) +DataArrayIdType *MEDCouplingUMesh::colinearize2D(double eps) +{ + return internalColinearize2D(eps, false); +} + +/*! + * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells. + * In a given 2D cell, if two edges are colinear and the junction point is used by a third edge, the two edges will not be + * merged, contrary to colinearize2D(). + * + * \sa MEDCouplingUMesh::colinearize2D + */ +DataArrayIdType *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps) +{ + return internalColinearize2D(eps, true); +} + + +/*! + * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it + */ +DataArrayIdType *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform) { - MCAuto ret(DataArrayInt::New()); ret->alloc(0,1); + MCAuto ret(DataArrayIdType::New()); ret->alloc(0,1); checkConsistencyLight(); if(getSpaceDimension()!=2 || getMeshDimension()!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !"); INTERP_KERNEL::QuadraticPlanarPrecision prec(eps); - INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps); - int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()); - const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin()); - MCAuto newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0); + mcIdType nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()); + const mcIdType *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin()); + MCAuto newc(DataArrayIdType::New()),newci(DataArrayIdType::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0); + std::map forbiddenPoints; // list of points that can not be removed (or it will break conformity) + if(stayConform) // + { + // A point that is used by more than 2 edges can not be removed without breaking conformity: + MCAuto desc1(DataArrayIdType::New()),descI1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescI1(DataArrayIdType::New()); + MCAuto mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1)); + MCAuto desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New()); + MCAuto mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2)); + MCAuto dsi(revDescI2->deltaShiftIndex()); + MCAuto ids(dsi->findIdsGreaterThan(2)); + const mcIdType * cPtr(mDesc0D->getNodalConnectivity()->begin()); + for(auto it = ids->begin(); it != ids->end(); it++) + forbiddenPoints[cPtr[2*(*it)+1]] = true; // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...] + } + MCAuto appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug. const double *coords(_coords->begin()); - int *newciptr(newci->getPointer()); - for(int i=0;igetPointer()); + for(mcIdType i=0;ipushBackSilent(i); newciptr[1]=newc->getNumberOfTuples(); } @@ -2039,19 +2238,19 @@ DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps) * startNode, endNode in global numbering *\return true if the segment is indeed split */ -bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode, - const int * c, const int * cI, const int *idsBg, const int *endBg, - std::vector & pointIds, std::vector & hitSegs) +bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, mcIdType startNode, mcIdType endNode, + const mcIdType * c, const mcIdType * cI, const mcIdType *idsBg, const mcIdType *endBg, + std::vector & pointIds, std::vector & hitSegs) { using namespace std; const int SPACEDIM=3; - typedef pair PairDI; + typedef pair PairDI; set< PairDI > x; - for (const int * it = idsBg; it != endBg; ++it) + for (const mcIdType * it = idsBg; it != endBg; ++it) { assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2); - int start = c[cI[*it]+1], end = c[cI[*it]+2]; + mcIdType start = c[cI[*it]+1], end = c[cI[*it]+2]; x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords x.insert(make_pair(coo[end*SPACEDIM], end)); } @@ -2061,10 +2260,10 @@ bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, i pointIds.reserve(xx.size()); // Keep what is inside [startNode, endNode]: - int go = 0; + mcIdType go = 0; for (vector::const_iterator it=xx.begin(); it != xx.end(); ++it) { - const int idx = (*it).second; + const mcIdType idx = (*it).second; if (!go) { if (idx == startNode) go = 1; @@ -2077,7 +2276,7 @@ bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, i break; } -// vector pointIds2(pointIds.size()+2); +// vector pointIds2(pointIds.size()+2); // copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1); // pointIds2[0] = startNode; // pointIds2[pointIds2.size()-1] = endNode; @@ -2086,12 +2285,12 @@ bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, i reverse(pointIds.begin(), pointIds.end()); // Now identify smaller segments that are not sub-divided - those won't need any further treatment: - for (const int * it = idsBg; it != endBg; ++it) + for (const mcIdType * it = idsBg; it != endBg; ++it) { - int start = c[cI[*it]+1], end = c[cI[*it]+2]; - vector::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start); + mcIdType start = c[cI[*it]+1], end = c[cI[*it]+2]; + vector::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start); if (itStart == pointIds.end()) continue; - vector::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end); + vector::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end); if (itEnd == pointIds.end()) continue; if (abs(distance(itEnd, itStart)) != 1) continue; hitSegs.push_back(*it); // segment is undivided. @@ -2100,23 +2299,23 @@ bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, i return (pointIds.size() > 2); // something else apart start and end node } -void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode, - const std::vector& insidePoints, std::vector& modifiedFace) +void MEDCouplingUMesh::ReplaceEdgeInFace(const mcIdType * sIdxConn, const mcIdType * sIdxConnE, mcIdType startNode, mcIdType endNode, + const std::vector& insidePoints, std::vector& modifiedFace) { using namespace std; - int dst = distance(sIdxConn, sIdxConnE); + size_t dst = distance(sIdxConn, sIdxConnE); modifiedFace.reserve(dst + insidePoints.size()-2); modifiedFace.resize(dst); copy(sIdxConn, sIdxConnE, modifiedFace.data()); - vector::iterator shortEnd = modifiedFace.begin()+dst; - vector::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode); + vector::iterator shortEnd = modifiedFace.begin()+dst; + vector::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode); if (startPos == shortEnd) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!"); - vector::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode); + vector::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode); if (endPos == shortEnd) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!"); - int d = distance(startPos, endPos); + size_t d = distance(startPos, endPos); if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ... modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end()); // insidePoints also contains start and end node. Those don't need to be inserted. else @@ -2137,16 +2336,20 @@ void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxC * This method expects that all nodes in \a this are not closer than \a eps. * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method. * - * \param [in] eps the relative error to detect merged edges. - * \return DataArrayInt * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array + * \b WARNING this method only works for 'partition-like' non-conformities. When joining adjacent faces, the set of all small faces must fit exactly into the + * neighboring bigger face. No real face intersection is actually computed. + * + * \param [in] eps the absolute (geometric) error to detect merged edges. + * \return DataArrayIdType * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array * that the user is expected to deal with. * * \throw If \a this is not coherent. * \throw If \a this has not spaceDim equal to 3. * \throw If \a this has not meshDim equal to 3. + * \throw If the 'partition-like' condition described above is not respected. * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly() */ -DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) +DataArrayIdType *MEDCouplingUMesh::conformize3D(double eps) { using namespace std; @@ -2159,23 +2362,24 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) MCAuto connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex())); const double * coo(_coords->begin()); - MCAuto ret(DataArrayInt::New()); + MCAuto ret(DataArrayIdType::New()); ret->alloc(0,1); + const double carMeshSz = getCaracteristicDimension(); { /************************* * STEP 1 -- faces *************************/ - MCAuto descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New()); + MCAuto descDNU(DataArrayIdType::New()),descIDNU(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New()); MCAuto mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI)); - const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); - const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); + const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); + const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); MCAuto connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity())); // Build BBTree MCAuto bboxArr(mDesc->getBoundingBoxForBBTree(eps)); const double *bbox(bboxArr->begin()); getCoords()->begin(); - int nDescCell(mDesc->getNumberOfCells()); - BBTree myTree(bbox,0,0,nDescCell,-eps); + mcIdType nDescCell=mDesc->getNumberOfCells(); + BBTree myTree(bbox,0,0,nDescCell,-eps); // Surfaces - handle biggest first MCAuto surfF = mDesc->getMeasureField(true); DataArrayDouble * surfs = surfF->getArray(); @@ -2185,35 +2389,35 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) const double * normalsP = normals->getConstPointer(); // Sort faces by decreasing surface: - vector< pair > S; - for(std::size_t i=0;i < surfs->getNumberOfTuples();i++) + vector< pair > S; + for(mcIdType i=0;i < surfs->getNumberOfTuples();i++) { - pair p = make_pair(surfs->begin()[i], i); + pair p = make_pair(surfs->begin()[i], i); S.push_back(p); } sort(S.rbegin(),S.rend()); // reverse sort vector hit(nDescCell); fill(hit.begin(), hit.end(), false); - vector hitPoly; // the final result: which 3D cells have been modified. + vector hitPoly; // the final result: which 3D cells have been modified. - for( vector >::const_iterator it = S.begin(); it != S.end(); it++) + for( vector >::const_iterator it = S.begin(); it != S.end(); it++) { - int faceIdx = (*it).second; + mcIdType faceIdx = (*it).second; if (hit[faceIdx]) continue; - vector candidates, cands2; + vector candidates, cands2; myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates); // Keep only candidates whose normal matches the normal of current face - for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) + for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) { - bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps); + bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps/carMeshSz); // using rough relative epsilon if (*it2 != faceIdx && col) cands2.push_back(*it2); } if (!cands2.size()) continue; - // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later + // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes one day INTERP_KERNEL::TranslationRotationMatrix rotation; INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]), coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]), @@ -2223,16 +2427,16 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) MCAuto mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called double * cooPartRef(mPartRef->_coords->getPointer()); double * cooPartCand(mPartCand->_coords->getPointer()); - for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) rotation.transform_vector(cooPartRef+SPACEDIM*ii); - for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) + for (mcIdType ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) rotation.transform_vector(cooPartCand+SPACEDIM*ii); // Localize faces in 2D thanks to barycenters MCAuto baryPart = mPartCand->computeCellCenterOfMass(); - vector compo; compo.push_back(2); + vector compo; compo.push_back(2); MCAuto baryPartZ = baryPart->keepSelectedComponents(compo); - MCAuto idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps); + MCAuto idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps); if (!idsGoodPlane->getNumberOfTuples()) continue; @@ -2241,15 +2445,15 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) compo[0] = 0; compo.push_back(1); MCAuto baryPartXY = baryPart->keepSelectedComponents(compo); mPartRef->changeSpaceDimension(2,0.0); - MCAuto cc(DataArrayInt::New()), ccI(DataArrayInt::New()); + MCAuto cc(DataArrayIdType::New()), ccI(DataArrayIdType::New()); mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI); if (!cc->getNumberOfTuples()) continue; - MCAuto dsi = ccI->deltaShiftIndex(); + MCAuto dsi = ccI->deltaShiftIndex(); { - MCAuto tmp = dsi->findIdsInRange(0, 2); + MCAuto tmp = dsi->findIdsInRange(0, 2); if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples()) { ostringstream oss; @@ -2258,16 +2462,16 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) } } - MCAuto ids = dsi->findIdsEqual(1); + MCAuto ids = dsi->findIdsEqual(1); // Boundary face: if (!ids->getNumberOfTuples()) continue; double checkSurf=0.0; - const int * idsGoodPlaneP(idsGoodPlane->begin()); - for (const int * ii = ids->begin(); ii != ids->end(); ii++) + const mcIdType * idsGoodPlaneP(idsGoodPlane->begin()); + for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++) { - int faceIdx2 = cands2[idsGoodPlaneP[*ii]]; + mcIdType faceIdx2 = cands2[idsGoodPlaneP[*ii]]; hit[faceIdx2] = true; checkSurf += surfs->begin()[faceIdx2]; } @@ -2279,18 +2483,18 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) } // For all polyhedrons using this face, replace connectivity: - vector polyIndices, packsIds, facePack; - for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++) + vector polyIndices, packsIds, facePack; + for (mcIdType ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++) polyIndices.push_back(revDescP[ii]); ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size()); // Current face connectivity - const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1; - const int * sIdxConnE = cDesc + cIDesc[faceIdx+1]; + const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1; + const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1]; connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds); // Deletion of old faces - int jj=0; - for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj) + mcIdType jj=0; + for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj) { if (packsIds[jj] == -1) // The below should never happen - if a face is used several times, with a different layout of the nodes @@ -2301,22 +2505,22 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) connSla->deletePack(*it2, packsIds[jj]); } // Insertion of new faces: - for (const int * ii = ids->begin(); ii != ids->end(); ii++) + for (const mcIdType * ii = ids->begin(); ii != ids->end(); ii++) { // Build pack from the face to insert: - int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)]; - int facePack2Sz; - const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type! + mcIdType faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)]; + mcIdType facePack2Sz; + const mcIdType * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type! // Insert it in all hit polyhedrons: - for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2) + for (vector::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2) connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz); // without the type } } } // end step1 // Set back modified connectivity - MCAuto cAuto; cAuto.takeRef(_nodal_connec); - MCAuto cIAuto; cIAuto.takeRef(_nodal_connec_index); + MCAuto cAuto; cAuto.takeRef(_nodal_connec); + MCAuto cIAuto; cIAuto.takeRef(_nodal_connec_index); connSla->convertToPolyhedronConn(cAuto, cIAuto); { @@ -2326,36 +2530,36 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) // Now we have a face-conform mesh. // Recompute descending - MCAuto desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New()); + MCAuto desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),revDesc(DataArrayIdType::New()),revDescI(DataArrayIdType::New()); // Rebuild desc connectivity with orientation this time!! MCAuto mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI)); - const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); - const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer()); - const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); - MCAuto ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy()); - MCAuto cDeepC(mDesc->getNodalConnectivity()->deepCopy()); + const mcIdType *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer()); + const mcIdType *descIP(descI->getConstPointer()), *descP(desc->getConstPointer()); + const mcIdType *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin()); + MCAuto ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy()); + MCAuto cDeepC(mDesc->getNodalConnectivity()->deepCopy()); MCAuto connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC)); - MCAuto desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New()); + MCAuto desc2(DataArrayIdType::New()),descI2(DataArrayIdType::New()),revDesc2(DataArrayIdType::New()),revDescI2(DataArrayIdType::New()); MCAuto mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2); // std::cout << "writing!\n"; // mDesc->writeVTK("/tmp/toto_desc_confInter.vtu"); // mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu"); - const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer()); - const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin()); + const mcIdType *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer()); + const mcIdType *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin()); MCAuto bboxArr(mDesc2->getBoundingBoxForBBTree(eps)); const double *bbox2(bboxArr->begin()); - int nDesc2Cell=mDesc2->getNumberOfCells(); - BBTree myTree2(bbox2,0,0,nDesc2Cell,-eps); + mcIdType nDesc2Cell=mDesc2->getNumberOfCells(); + BBTree myTree2(bbox2,0,0,nDesc2Cell,-eps); // Edges - handle longest first MCAuto lenF = mDesc2->getMeasureField(true); DataArrayDouble * lens = lenF->getArray(); // Sort edges by decreasing length: - vector > S; - for(std::size_t i=0;i < lens->getNumberOfTuples();i++) + vector > S; + for(mcIdType i=0;i < lens->getNumberOfTuples();i++) { - pair p = make_pair(lens->getIJ(i, 0), i); + pair p = make_pair(lens->getIJ(i, 0), i); S.push_back(p); } sort(S.rbegin(),S.rend()); // reverse sort @@ -2363,26 +2567,30 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) vector hit(nDesc2Cell); fill(hit.begin(), hit.end(), false); - for( vector >::const_iterator it = S.begin(); it != S.end(); it++) + for( vector >::const_iterator it = S.begin(); it != S.end(); it++) { - int eIdx = (*it).second; + mcIdType eIdx = (*it).second; if (hit[eIdx]) continue; - vector candidates, cands2; + vector candidates, cands2; myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates); // Keep only candidates colinear with current edge double vCurr[3]; - unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2]; - for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar? + mcIdType start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2]; + for (mcIdType i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar? vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3]; - for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) + double carSz = INTERP_KERNEL::caracteristicDimVector(vCurr), eps2 = eps*carSz*carSz*carSz; + for(vector::const_iterator it2=candidates.begin();it2!=candidates.end();it2++) { double vOther[3]; - unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2]; - for (int i3=0; i3 < 3; i3++) + mcIdType start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2]; + for (mcIdType i3=0; i3 < 3; i3++) vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3]; - bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps); + // isColinear() expect unit vecotr, and relative (non geometrical) precision. + // relative prec means going from eps -> eps/curSz + // normalizing vCurr and vOther means multiplying by v^4, knowing that inside isColinear3D() the square (^2) of a dot product (^2) is computed + bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps2); // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine()) if (col) @@ -2392,69 +2600,68 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) continue; // Now rotate edges to bring them on Ox - int startNode = cDesc2[cIDesc2[eIdx]+1]; - int endNode = cDesc2[cIDesc2[eIdx]+2]; + mcIdType startNode = cDesc2[cIDesc2[eIdx]+1]; + mcIdType endNode = cDesc2[cIDesc2[eIdx]+2]; INTERP_KERNEL::TranslationRotationMatrix rotation; INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation); MCAuto mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called - MCAuto mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called - MCAuto nodeMap(mPartCand->zipCoordsTraducer()); - int nbElemsNotM1; + MCAuto mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is NOT called + MCAuto nodeMap(mPartCand->zipCoordsTraducer()); + mcIdType nbElemsNotM1; { - MCAuto tmp(nodeMap->findIdsNotEqual(-1)); + MCAuto tmp(nodeMap->findIdsNotEqual(-1)); nbElemsNotM1 = tmp->getNbOfElems(); } - MCAuto nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1); + MCAuto nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1); double * cooPartRef(mPartRef->_coords->getPointer()); double * cooPartCand(mPartCand->_coords->getPointer()); - for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) + for (mcIdType ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++) rotation.transform_vector(cooPartRef+SPACEDIM*ii); - for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) + for (mcIdType ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++) rotation.transform_vector(cooPartCand+SPACEDIM*ii); - - // Eliminate all edges for which y or z is not null + // Eliminate all edges for which y or z is not null - those are colinear edges which are simply parallels, but not on the same line MCAuto baryPart = mPartCand->computeCellCenterOfMass(); - vector compo; compo.push_back(1); + vector compo; compo.push_back(1); MCAuto baryPartY = baryPart->keepSelectedComponents(compo); compo[0] = 2; MCAuto baryPartZ = baryPart->keepSelectedComponents(compo); - MCAuto idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps); - MCAuto idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps); - MCAuto idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2); + MCAuto idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps); + MCAuto idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps); + MCAuto idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2); if (!idsGoodLine->getNumberOfTuples()) continue; // Now the ordering along the Ox axis: - std::vector insidePoints, hitSegs; + std::vector insidePoints, hitSegs; bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode], mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(), idsGoodLine->begin(), idsGoodLine->end(), /*out*/insidePoints, hitSegs); // Optim: smaller segments completely included in eIdx and not split won't need any further treatment: - for (vector::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its) + for (vector::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its) hit[cands2[*its]] = true; if (!isSplit) // current segment remains in one piece continue; // Get original node IDs in global coords array - for (std::vector::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit) + for (std::vector::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit) *iit = nodeMapInv->begin()[*iit]; - vector polyIndices, packsIds, facePack; + vector polyIndices, packsIds, facePack; // For each face implying this edge - for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++) + for (mcIdType ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++) { - int faceIdx = revDescP2[ii]; + mcIdType faceIdx = revDescP2[ii]; // each cell where this face is involved connectivity will be modified: ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]); // Current face connectivity - const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1; - const int * sIdxConnE = cDesc + cIDesc[faceIdx+1]; + const mcIdType * sIdxConn = cDesc + cIDesc[faceIdx] + 1; + const mcIdType * sIdxConnE = cDesc + cIDesc[faceIdx+1]; - vector modifiedFace; + vector modifiedFace; ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace); modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON); connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size()); @@ -2463,22 +2670,23 @@ DataArrayInt *MEDCouplingUMesh::conformize3D(double eps) // Rebuild 3D connectivity from descending: MCAuto newConn(MEDCouplingSkyLineArray::New()); - MCAuto superIdx(DataArrayInt::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0); - MCAuto idx(DataArrayInt::New()); idx->alloc(1); idx->fillWithValue(0); - MCAuto vals(DataArrayInt::New()); vals->alloc(0); + MCAuto superIdx(DataArrayIdType::New()); superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0); + MCAuto idx(DataArrayIdType::New()); idx->alloc(1); idx->fillWithValue(0); + MCAuto vals(DataArrayIdType::New()); vals->alloc(0); newConn->set3(superIdx, idx, vals); - for(std::size_t ii = 0; ii < getNumberOfCells(); ii++) - for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++) + mcIdType nbCells=getNumberOfCells(); + for(mcIdType ii = 0; ii < nbCells; ii++) + for (mcIdType jj=descIP[ii]; jj < descIP[ii+1]; jj++) { - int sz, faceIdx = abs(descP[jj])-1; + mcIdType sz, faceIdx = abs(descP[jj])-1; bool orient = descP[jj]>0; - const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz); + const mcIdType * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz); if (orient) newConn->pushBackPack(ii, p+1, p+sz); // +1 to skip type else { - vector rev(sz-1); - for (int kk=0; kk rev(sz-1); + for (mcIdType kk=0; kkpushBackPack(ii, rev.data(), rev.data()+sz-1); } }