From bc460c6b18a13410f9a14ba8cf633e503fd6d147 Mon Sep 17 00:00:00 2001 From: geay Date: Fri, 4 Jul 2014 13:25:20 +0200 Subject: [PATCH] All adrien tests are OK with this version. --- src/MEDCoupling/MEDCouplingUMesh.cxx | 521 +++++++++++++++--- src/MEDCoupling/MEDCouplingUMesh.hxx | 1 + src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 146 +++++ src/MEDCoupling_Swig/MEDCouplingCommon.i | 2 + 4 files changed, 588 insertions(+), 82 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 1f5006af9..dc0dcf5b1 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -8792,6 +8792,8 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 return ret.retn(); } +/// @cond INTERNAL + bool IsColinearOfACellOf(const std::vector< std::vector >& intersectEdge1, const std::vector& candidates, int start, int stop, int& retVal) { if(candidates.empty()) @@ -8815,7 +8817,8 @@ bool IsColinearOfACellOf(const std::vector< std::vector >& intersectEdge1, return false; } -//tony to put in private of MEDCouplingUMesh +/// @endcond + 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, MEDCouplingAutoRefCountObjectPtr& idsInRetColinear, MEDCouplingAutoRefCountObjectPtr& idsInMesh1DForIdsInRetColinear) { @@ -8948,41 +8951,312 @@ void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector& conn, con } } -MEDCouplingUMesh *BuildMesh2DCutFrom(int cellIdInMesh2D, const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh2DDesc, const MEDCouplingUMesh *splitMesh1D, - const int *descBg, const int *descEnd, const std::vector< std::vector >& intersectEdge1, int offset, - MEDCouplingAutoRefCountObjectPtr& idsLeftRight) +/*! + * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve. + */ +void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector& edge1Bis, const std::vector< MEDCouplingAutoRefCountObjectPtr >& edge1BisPtr, + std::vector< std::vector >& out0, std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr > >& out1) { - idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(splitMesh1D->getNumberOfCells(),2); - int *idsLeftRightPtr(idsLeftRight->getPointer()); - MEDCouplingAutoRefCountObjectPtr ret(MEDCouplingUMesh::New("",2)); - ret->setCoords(splitMesh1D->getCoords()); - ret->allocateCells(); - int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells()); - if(nbCellsInSplitMesh1D==0) - throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error ! input 1D mesh must have at least one cell !"); - const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin()), - *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin()); + std::size_t nb(edge1Bis.size()/2); + std::size_t nbOfEdgesOf2DCellSplit(nb/2); + int 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()); + for(ii=0;ii allEdges; - std::vector< MEDCouplingAutoRefCountObjectPtr > allEdgesPtr; - for(const int *it(descBg);it!=descEnd;it++) + if(jj==nb) + {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell + out0.resize(1); out1.resize(1); + std::vector& connOut(out0[0]); + connOut.resize(nbOfEdgesOf2DCellSplit); + std::vector< MEDCouplingAutoRefCountObjectPtr >& edgesPtr(out1[0]); + edgesPtr.resize(nbOfEdgesOf2DCellSplit); + for(std::size_t kk=0;kk,int> m; - MEDCouplingAutoRefCountObjectPtr ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m)); - const std::vector& edge1(intersectEdge1[edgeId]); - if(*it>0) - allEdges.insert(allEdges.end(),edge1.begin(),edge1.end()); + // [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< MEDCouplingAutoRefCountObjectPtr >& eleft(out1[0]); + std::vector< MEDCouplingAutoRefCountObjectPtr >& eright(out1[1]); + for(std::size_t k=ii;k > ees(iEnd); + for(int ik=iEnd-1;ik>=0;ik--) + { + std::map< MEDCouplingAutoRefCountObjectPtr,int> m; + MEDCouplingAutoRefCountObjectPtr ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m)); + ees[iEnd-1-ik]=ee; + } + for(int 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< MEDCouplingAutoRefCountObjectPtr >& edgesPtr); +public: + std::vector _edges; + std::vector< MEDCouplingAutoRefCountObjectPtr > _edges_ptr; +}; + +CellInfo::CellInfo(const std::vector& edges, const std::vector< MEDCouplingAutoRefCountObjectPtr >& edgesPtr) +{ + std::size_t nbe(edges.size()); + std::vector edges2(2*nbe); std::vector< MEDCouplingAutoRefCountObjectPtr > edgesPtr2(2*nbe); + for(std::size_t i=0;i& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { } + EdgeInfo(int istart, int iend, int pos, const MEDCouplingAutoRefCountObjectPtr& 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< MEDCouplingAutoRefCountObjectPtr >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr >& newRight); + void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const; +private: + int _istart; + int _iend; + MEDCouplingAutoRefCountObjectPtr _mesh; + MEDCouplingAutoRefCountObjectPtr _edge; + int _left; + int _right; +}; + +void EdgeInfo::somethingHappendAt(int pos, const std::vector< MEDCouplingAutoRefCountObjectPtr >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr >& newRight) +{ + const MEDCouplingUMesh *mesh(_mesh); + if(mesh) + return ; + if(_rightpos) + { _left++; _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()); + if((isLeft && isRight) || (!isLeft && !isRight)) + throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !"); + if(isLeft) + return ; + if(isRight) + { + _right++; + return ; + } + } + if(_left==pos) + { + bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end()); + if((isLeft && isRight) || (!isLeft && !isRight)) + throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !"); + if(isLeft) + { + _right++; + return ; + } + if(isRight) + { + _left++; + _right++; + return ; + } + } +} + +void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const +{ + const MEDCouplingUMesh *mesh(_mesh); + if(!mesh) + { + neighbors[0]=offset+_left; neighbors[1]=offset+_right; + } + else + {// not fully splitting cell case + if(mesh2D->getNumberOfCells()==1) + {//little optimization. 1 cell no need to find in which cell mesh is ! + neighbors[0]=offset; neighbors[1]=offset; + return; + } else - allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend()); - std::size_t sz(edge1.size()); - for(std::size_t cnt=0;cnt barys(mesh->getBarycenterAndOwner()); + int 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; + } } - std::size_t nb(allEdges.size()),ii,jj; +} + +class VectorOfCellInfo +{ +public: + VectorOfCellInfo(const std::vector& edges, const std::vector< MEDCouplingAutoRefCountObjectPtr >& edgesPtr); + std::size_t size() const { return _pool.size(); } + int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const; + void setMeshAt(int pos, const MEDCouplingAutoRefCountObjectPtr& mesh, int istart, int iend, const MEDCouplingAutoRefCountObjectPtr& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr > >& edgePtrs); + const std::vector& getConnOf(int pos) const { return get(pos)._edges; } + const std::vector< MEDCouplingAutoRefCountObjectPtr >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; } + MEDCouplingAutoRefCountObjectPtr getZeMesh() const { return _ze_mesh; } + void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const; +private: + int getZePosOfEdgeGivenItsGlobalId(int pos) const; + void updateEdgeInfo(int pos, const std::vector< MEDCouplingAutoRefCountObjectPtr >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr >& newRight); + const CellInfo& get(int pos) const; + CellInfo& get(int pos); +private: + std::vector _pool; + MEDCouplingAutoRefCountObjectPtr _ze_mesh; + std::vector _edge_info; +}; + +VectorOfCellInfo::VectorOfCellInfo(const std::vector& edges, const std::vector< MEDCouplingAutoRefCountObjectPtr >& edgesPtr):_pool(1) +{ + _pool[0]._edges=edges; + _pool[0]._edges_ptr=edgesPtr; +} + +int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const +{ + if(_pool.empty()) + throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !"); + if(_pool.size()==1) + return 0; + const MEDCouplingUMesh *zeMesh(_ze_mesh); + if(!zeMesh) + throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !"); + MEDCouplingAutoRefCountObjectPtr barys(mesh->getBarycenterAndOwner()); + return zeMesh->getCellContainingPoint(barys->begin(),eps); +} + +void VectorOfCellInfo::setMeshAt(int pos, const MEDCouplingAutoRefCountObjectPtr& mesh, int istart, int iend, const MEDCouplingAutoRefCountObjectPtr& mesh1DInCase, const std::vector< std::vector >& edges, const std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr > >& edgePtrs) +{ + get(pos);//to check pos + bool isFast(pos==0 && _pool.size()==1); + std::size_t sz(edges.size()); + // dealing with edges + if(sz==1) + _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase)); + else + _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back())); + // + std::vector pool(_pool.size()-1+sz); + for(int i=0;i > ms; + if(pos>0) + { + MEDCouplingAutoRefCountObjectPtr elt(static_cast(_ze_mesh->buildPartOfMySelf2(0,pos,true))); + ms.push_back(elt); + } + ms.push_back(mesh); + if(pos<_ze_mesh->getNumberOfCells()-1) + { + MEDCouplingAutoRefCountObjectPtr elt(static_cast(_ze_mesh->buildPartOfMySelf2(pos+1,_ze_mesh->getNumberOfCells(),true))); + ms.push_back(elt); + } + std::vector< const MEDCouplingUMesh *> ms2(ms.size()); + for(std::size_t j=0;j=0 !"); + int ret(0); + for(std::vector::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++) + { + if((*it).isInMyRange(pos)) + return ret; + } + throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !"); +} + +void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MEDCouplingAutoRefCountObjectPtr >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr >& newRight) +{ + get(pos);//to check; + if(_edge_info.empty()) + return ; + std::size_t sz(_edge_info.size()-1); + for(std::size_t i=0;i=(int)_pool.size()) + throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !"); + return _pool[pos]; +} + +CellInfo& VectorOfCellInfo::get(int pos) +{ + if(pos<0 || pos>=(int)_pool.size()) + throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !"); + return _pool[pos]; +} + +MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *splitMesh1D, const std::vector& allEdges, const std::vector< MEDCouplingAutoRefCountObjectPtr >& allEdgesPtr, int offset, + MEDCouplingAutoRefCountObjectPtr& idsLeftRight) +{ + int 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()); + std::size_t nb(allEdges.size()),jj; if(nb%2!=0) throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !"); - std::size_t nbOfEdgesOf2DCellSplit(nb/2); std::vector edge1Bis(nb*2); std::vector< MEDCouplingAutoRefCountObjectPtr > edge1BisPtr(nb*2); std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()); @@ -8990,8 +9264,11 @@ MEDCouplingUMesh *BuildMesh2DCutFrom(int cellIdInMesh2D, const MEDCouplingUMesh 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()); + VectorOfCellInfo pool(edge1Bis,edge1BisPtr); for(int iStart=0;iStart partOfSplitMesh1D(static_cast(splitMesh1D->buildPartOfMySelf2(iStart,iEnd,1,true))); + int pos(pool.getPositionOf(eps,partOfSplitMesh1D)); // - if(jj==nb) - {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell - std::vector connOut(nbOfEdgesOf2DCellSplit); - std::vector< MEDCouplingAutoRefCountObjectPtr > edgesPtr(nbOfEdgesOf2DCellSplit); - for(std::size_t kk=0;kk connOutLeft,connOutRight;//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1] - std::vector< MEDCouplingAutoRefCountObjectPtr > eleft,eright; - for(std::size_t k=ii;k=iStart;ik--) - { - connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]); - std::map< MEDCouplingAutoRefCountObjectPtr,int> m; - MEDCouplingAutoRefCountObjectPtr ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m)); - eleft.push_back(ee); - } - for(std::size_t k=jj+1;k,int> m; - connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]); - MEDCouplingAutoRefCountObjectPtr ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m)); - eright.push_back(ee); - } - AddCellInMesh2D(ret,connOutLeft,eleft); - AddCellInMesh2D(ret,connOutRight,eright); - for(int mm=iStart;mmretTmp(MEDCouplingUMesh::New("",2)); + retTmp->setCoords(splitMesh1D->getCoords()); + retTmp->allocateCells(); + + std::vector< std::vector > out0; + std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr > > out1; + + BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1); + for(std::size_t cnt=0;cnt >& intersectEdge1, int offset, + MEDCouplingAutoRefCountObjectPtr& idsLeftRight) +{ + const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin()); + // + std::vector allEdges; + std::vector< MEDCouplingAutoRefCountObjectPtr > allEdgesPtr; + for(const int *it(descBg);it!=descEnd;it++) + { + int edgeId(std::abs(*it)-1); + std::map< MEDCouplingAutoRefCountObjectPtr,int> m; + MEDCouplingAutoRefCountObjectPtr ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m)); + const std::vector& edge1(intersectEdge1[edgeId]); + if(*it>0) + allEdges.insert(allEdges.end(),edge1.begin(),edge1.end()); + else + allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend()); + std::size_t sz(edge1.size()); + for(std::size_t cnt=0;cnt idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCpy()); idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1); idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique(); - MEDCouplingAutoRefCountObjectPtr out0s;//ids in mesh2D that are impacted by the fact that some edges of \a mesh1D are part of those cells + MEDCouplingAutoRefCountObjectPtr 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); @@ -9197,7 +9468,7 @@ void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, MEDCouplingAutoRefCountObjectPtr idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end())); MEDCouplingAutoRefCountObjectPtr partOfMesh1CuttingCur2DCell(static_cast(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end()))); MEDCouplingAutoRefCountObjectPtr partOfRet3; - MEDCouplingAutoRefCountObjectPtr splitOfOneCell(BuildMesh2DCutFrom(*it,mesh2D,m1Desc,partOfMesh1CuttingCur2DCell,dd1->begin()+dd2->getIJ(*it,0),dd1->begin()+dd2->getIJ((*it)+1,0),intersectEdge1,ret2->getNumberOfTuples(),partOfRet3)); + MEDCouplingAutoRefCountObjectPtr 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(int i=0;igetNumberOfCells();i++) @@ -9306,6 +9577,90 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo } } +/** + * Provides a renumbering of the cells of this (which has to be a piecewise connected 1D line), so that + * the segments of the line are indexed in consecutive order (i.e. cells \a i and \a i+1 are neighbors). + * This doesn't modify the mesh. + * The caller is to deal with the resulting DataArrayInt. + * \throw If the coordinate array is not set. + * \throw If the nodal connectivity of the cells is not defined. + * \throw If m1 is not a mesh of dimension 2, or m1 is not a mesh of dimension 1 + * \throw If m2 is not a (piecewise) line (i.e. if a point has more than 2 adjacent segments) + */ +DataArrayInt *MEDCouplingUMesh::orderConsecutiveCells1D() const +{ + checkFullyDefined(); + if(getMeshDimension()!=1 || getSpaceDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::orderConsecutiveCells1D works on unstructured mesh with (meshdim, spacedim) = (1,2)!"); + + // Check that this is a line (and not a more complex 1D mesh) - each point is used at most by 2 segments: + MEDCouplingAutoRefCountObjectPtr _d(DataArrayInt::New()),_dI(DataArrayInt::New()); + MEDCouplingAutoRefCountObjectPtr _rD(DataArrayInt::New()),_rDI(DataArrayInt::New()); + MEDCouplingAutoRefCountObjectPtr m_points(buildDescendingConnectivity(_d, _dI, _rD, _rDI)); + const int *d(_d->getConstPointer()), *dI(_dI->getConstPointer()); + const int *rD(_rD->getConstPointer()), *rDI(_rDI->getConstPointer()); + MEDCouplingAutoRefCountObjectPtr _dsi(_rDI->deltaShiftIndex()); + const int * dsi(_dsi->getConstPointer()); + MEDCouplingAutoRefCountObjectPtr dsii = _dsi->getIdsNotInRange(0,3); + m_points=0; + if (dsii->getNumberOfTuples()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::orderConsecutiveCells1D only work with a mesh being a (piecewise) connected line!"); + + int nc(getNumberOfCells()); + MEDCouplingAutoRefCountObjectPtr result(DataArrayInt::New()); + result->alloc(nc,1); + + // set of edges not used so far + std::set edgeSet; + for (int i=0; i linePiece; + // fills a list of consecutive segment linked to startSeg. This can go forward or backward. + for (int direction=0;direction<2;direction++) // direction=0 --> forward, direction=1 --> backward + { + // Fill the list forward (resp. backward) from the start segment: + int activeSeg = startSeg; + int prevPointId = -20; + int ptId; + while (!edgeSet.empty()) + { + if (!(direction == 1 && prevPointId==-20)) // prevent adding twice startSeg + { + if (direction==0) + linePiece.push_back(activeSeg); + else + linePiece.push_front(activeSeg); + edgeSet.erase(activeSeg); + } + + int ptId1 = d[dI[activeSeg]], ptId2 = d[dI[activeSeg]+1]; + ptId = direction ? (ptId1 == prevPointId ? ptId2 : ptId1) : (ptId2 == prevPointId ? ptId1 : ptId2); + if (dsi[ptId] == 1) // hitting the end of the line + break; + prevPointId = ptId; + int seg1 = rD[rDI[ptId]], seg2 = rD[rDI[ptId]+1]; + activeSeg = (seg1 == activeSeg) ? seg2 : seg1; + } + } + // Done, save final piece into DA: + std::copy(linePiece.begin(), linePiece.end(), result->getPointer()+newIdx); + newIdx += linePiece.size(); + + // identify next valid start segment (one which is not consumed) + if(!edgeSet.empty()) + startSeg = *(edgeSet.begin()); + } + while (!edgeSet.empty()); + return result.retn(); +} + +/// @cond INTERNAL + void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map,int>& m, int forbVal0, int forbVal1, std::vector& isect) { MEDCouplingAutoRefCountObjectPtr nTmp(n); nTmp->incrRef(); @@ -9336,6 +9691,8 @@ bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map< return presenceOfOn; } +/// @endcond + /** * This method split some of edges of 2D cells in \a this. The edges to be split are specified in \a subNodesInSeg and in \a subNodesInSegI using index storage mode. * To do the work this method can optionally needs information about middle of subedges for quadratic cases if a minimal creation of new nodes is wanted. diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index e8e4f0c1c..f43c356bf 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -266,6 +266,7 @@ namespace ParaMEDMEM DataArrayInt *& commonCellsArr, DataArrayInt *& commonCellsIArr); MEDCOUPLING_EXPORT DataArrayInt *buildUnionOf2DMesh() const; MEDCOUPLING_EXPORT DataArrayInt *buildUnionOf3DMesh() const; + MEDCOUPLING_EXPORT DataArrayInt *orderConsecutiveCells1D() const; private: MEDCouplingUMesh(); MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCopy); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 787f4f27e..d4d92e319 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -15727,6 +15727,152 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(c.isEqual(DataArrayInt([0,2,3,4,5,7,8,9,11,12,14,15,1,1,6,6,10,10,13,13]))) self.assertTrue(d.isEqual(DataArrayInt([(12,13),(14,15),(16,17),(18,19)]))) pass + + def testIntersect2DMeshWith1DLine6(self): + """ Basic test for Intersect2DMeshWith1DLine: a vertical line intersecting a square. """ + m1c = MEDCouplingCMesh() + coordX = DataArrayDouble([-1., 1., 2]) + m1c.setCoordsAt(0,coordX) + coordY = DataArrayDouble([0., 2.]) + m1c.setCoordsAt(1,coordY); + m1 = m1c.buildUnstructured() + + # A simple line: + m2 = MEDCouplingUMesh("bla", 1) + coord2 = DataArrayDouble([0.,-1.0, 0.,1., 0.,3., 0.5,2.2], 4, 2) + conn2 = DataArrayInt([NORM_SEG2,0,1,NORM_SEG3,1,2,3]) + connI2 = DataArrayInt([0,3,7]) + m2.setCoords(coord2) + m2.setConnectivity(conn2, connI2) + + # End of construction of input meshes m1bis and m2 -> start of specific part of the test + a,b,c,d = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10) + self.assertTrue(a.getNodalConnectivity().isEqual(DataArrayInt([4,2,1,4,5,32,0,3,11,7,10,14,15,16,17,18,32,4,1,10,7,11,19,20,21,22,23]))) + self.assertTrue(a.getNodalConnectivityIndex().isEqual(DataArrayInt([0,5,16,27]))) + self.assertTrue(b.getNodalConnectivity().isEqual(DataArrayInt([1,6,10,1,10,7,2,7,11,12,2,11,8,13]))) + self.assertTrue(b.getNodalConnectivityIndex().isEqual(DataArrayInt([0,3,6,10,14]))) + self.assertTrue(a.getCoords().getHiddenCppPointer()==b.getCoords().getHiddenCppPointer()) + self.assertTrue(a.getCoords()[:6].isEqual(m1.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[6:10].isEqual(m2.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[10:].isEqual(DataArrayDouble([(0.,0.),(0.5164175471673584,2.),(0.3796918047064557,1.43726403104512),(0.3796918047064557,2.56273596895488),(-1.,1.),(-0.24179122641632078,2.),(0.3796918047064558,1.4372640310451201),(0.,0.5),(-0.5,0.),(1.,1.),(0.5,0.),(0.,0.5),(0.3796918047064558,1.4372640310451201),(0.7582087735836792,2.)]),1e-12)) + self.assertTrue(c.isEqual(DataArrayInt([1,0,0]))) + self.assertTrue(d.isEqual(DataArrayInt([(-1,-1),(1,2),(1,2),(-1,-1)]))) + pass + + def testSwig2Intersect2DMeshWith1DLine7(self): + """ Star pattern (a triangle intersecting another one upside down) """ + coords1 = DataArrayDouble([-2.,1., 2.,1., 0.,-2.], 3,2) + coords2 = DataArrayDouble([0.,2., 2.,-1., -2.,-1., 0.,3.], 4,2) + m1 = MEDCouplingUMesh("triangle", 2) + m2 = MEDCouplingUMesh("tri_line", 1) + m1.setCoords(coords1) + m2.setCoords(coords2) + m1.setConnectivity(DataArrayInt([NORM_TRI3, 0,1,2]), DataArrayInt([0,4])) + m2.setConnectivity(DataArrayInt([NORM_SEG2,0,1,NORM_SEG2,1,2,NORM_SEG2,2,3]), DataArrayInt([0,3,6,9])) + # End of construction of input meshes m1bis and m2 -> start of specific part of the test + a,b,c,d = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10) + self.assertTrue(a.getNodalConnectivity().isEqual(DataArrayInt([5,1,9,7,5,2,11,10,5,0,8,12,5,7,9,10,11,12,8]))) + self.assertTrue(a.getNodalConnectivityIndex().isEqual(DataArrayInt([0,4,8,12,19]))) + self.assertTrue(b.getNodalConnectivity().isEqual(DataArrayInt([1,3,7,1,7,9,1,9,4,1,4,10,1,10,11,1,11,5,1,5,12,1,12,8,1,8,6]))) + self.assertTrue(b.getNodalConnectivityIndex().isEqual(DataArrayInt([0,3,6,9,12,15,18,21,24,27]))) + self.assertTrue(a.getCoords()[:3].isEqual(m1.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[3:7].isEqual(m2.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[7:].isEqual(DataArrayDouble([(0.6666666666666666,1.),(-1.,1.),(1.3333333333333333,1.1102230246251565e-16),(0.6666666666666665,-0.9999999999999996),(-0.6666666666666667,-1.),(-1.4285714285714284,0.14285714285714302)]),1e-12)) + self.assertTrue(a.getCoords().getHiddenCppPointer()==b.getCoords().getHiddenCppPointer()) + self.assertTrue(c.isEqual(DataArrayInt([0,0,0,0]))) + self.assertTrue(d.isEqual(DataArrayInt([(-1,-1),(0,3),(-1,-1),(-1,-1),(1,3),(-1,-1),(-1,-1),(2,3),(-1,-1)]))) + pass + + def testSwig2Intersect2DMeshWith1DLine8(self): + """ Line pieces ending (or fully located) in the middle of a cell """ + m1c = MEDCouplingCMesh() + m1c.setCoordsAt(0,DataArrayDouble([-1., 1.])) + m1c.setCoordsAt(1,DataArrayDouble([-1., 1.])); + m1 = m1c.buildUnstructured() + coords2 = DataArrayDouble([0.,0., 0.,1.5, -1.5,0., 0.5,0.0, 0.0,-0.5, 1.1,-0.6], 6,2) + m2 = MEDCouplingUMesh("piecewise_line", 1) + m2.setCoords(coords2) + c = DataArrayInt([NORM_SEG2,2,1, NORM_SEG2,1,4, NORM_SEG2,4,3, NORM_SEG2,3,5]) + cI = DataArrayInt([0,3,6,9,12]) + m2.setConnectivity(c, cI) + a,b,c,d = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10) + self.assertTrue(a.getNodalConnectivity().isEqual(DataArrayInt([5,2,11,10,5,3,13,7,8,12,5,1,0,10,11,12,8,7,13]))) + self.assertTrue(a.getNodalConnectivityIndex().isEqual(DataArrayInt([0,4,10,19]))) + self.assertTrue(b.getNodalConnectivity().isEqual(DataArrayInt([1,6,10,1,10,11,1,11,5,1,5,12,1,12,8,1,8,7,1,7,13,1,13,9]))) + self.assertTrue(b.getNodalConnectivityIndex().isEqual(DataArrayInt([0,3,6,9,12,15,18,21,24]))) + self.assertTrue(a.getCoords()[:4].isEqual(m1.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[4:10].isEqual(m2.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[10:].isEqual(DataArrayDouble([(-1.,0.5),(-0.5,1.),(0.,1.),(1.,-0.5)]),1e-12)) + self.assertTrue(a.getCoords().getHiddenCppPointer()==b.getCoords().getHiddenCppPointer()) + self.assertTrue(c.isEqual(DataArrayInt([0,0,0]))) + self.assertTrue(d.isEqual(DataArrayInt([(-1,-1),(0,2),(-1,-1),(-1,-1),(1,2),(1,2),(1,2),(-1,-1)]))) + pass + + def testSwig2Intersect2DMeshWith1DLine9(self): + """ Intersection with a line whose connectivity is not consecutive """ + m1c = MEDCouplingCMesh() + coordX = DataArrayDouble([-1., 1., 2]) + m1c.setCoordsAt(0,coordX) + coordY = DataArrayDouble([0., 2.]) + m1c.setCoordsAt(1,coordY); + m1 = m1c.buildUnstructured() + # A simple line: + m2 = MEDCouplingUMesh("bla", 1) + coord2 = DataArrayDouble([0.,1.5, 0.5,1., 0.0,0.5, 0.0,3.0, 0.0,-1.0], 5, 2) + conn2 = DataArrayInt([NORM_SEG2,3,0,NORM_SEG3,0,2,1,NORM_SEG2,2,4]) + connI2 = DataArrayInt([0,3,7,10]) + m2.setCoords(coord2) + m2.setConnectivity(conn2, connI2) + # End of construction of input meshes m1bis and m2 -> start of specific part of the test + a,b,c,d = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10) + self.assertTrue(a.getNodalConnectivity().isEqual(DataArrayInt([4,2,1,4,5,32,4,1,11,8,6,12,14,15,16,17,18,19,32,0,3,12,6,8,11,20,21,22,23,24,25]))) + self.assertTrue(a.getNodalConnectivityIndex().isEqual(DataArrayInt([0,5,18,31]))) + self.assertTrue(b.getNodalConnectivity().isEqual(DataArrayInt([1,9,12,1,12,6,2,6,8,13,1,8,11,1,11,10]))) + self.assertTrue(b.getNodalConnectivityIndex().isEqual(DataArrayInt([0,3,6,10,13,16]))) + self.assertTrue(a.getCoords()[:6].isEqual(m1.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[6:11].isEqual(m2.getCoords(),1e-12)) + self.assertTrue(a.getCoords()[11:].isEqual(DataArrayDouble([(0.,0.),(0.,2.),(0.5,1.),(1.,1.),(0.5,0.),(0.,0.25),(0.5,1.),(0.,1.75),(0.5,2.),(-1.,1.),(-0.5,2.),(0.,1.75),(0.5,1.),(0.,0.25),(-0.5,0.)]),1e-12)) + self.assertTrue(a.getCoords().getHiddenCppPointer()==b.getCoords().getHiddenCppPointer()) + self.assertTrue(c.isEqual(DataArrayInt([1,0,0]))) + self.assertTrue(d.isEqual(DataArrayInt([(-1,-1),(1,2),(1,2),(1,2),(-1,-1)]))) + pass + + def testOrderConsecutiveCells1D1(self): + """A line in several unconnected pieces:""" + m2 = MEDCouplingUMesh.New("bla", 1) + c = DataArrayInt([NORM_SEG2,0,1,NORM_SEG3,1,3,2, NORM_SEG2,3,4, + NORM_SEG3,5,7,6, NORM_SEG3,7,9,8, NORM_SEG2,9,10, + NORM_SEG2,11,12,NORM_SEG2,12,13, + NORM_SEG2,14,15]) + cI = DataArrayInt([0,3,7,10,14,18,21,24,27,30]) + coords2 = DataArrayDouble([float(i) for i in range(32)], 16,2) + m2.setCoords(coords2); + m2.setConnectivity(c, cI); + m2.checkCoherency2(1.0e-8); + + # Shuffle a bit :-) + m2.renumberCells(DataArrayInt([0,3,6,8,1,4,7,5,2]), True); + res = m2.orderConsecutiveCells1D() + expRes = [0,3,6,8,1,4,2,7,5] + self.assertEqual(m2.getNumberOfCells(),res.getNumberOfTuples()) + self.assertEqual(expRes, res.getValues()) + + # A closed line (should also work) + m3 = MEDCouplingUMesh.New("bla3", 1) + conn3A = DataArrayInt([NORM_SEG2,0,1,NORM_SEG3,1,3,2, NORM_SEG2,3,0]) + coord3 = coords2[0:5] + c.reAlloc(10) + cI.reAlloc(4) + + m3.setCoords(coord3) + m3.setConnectivity(conn3A, cI) + m3.checkCoherency2(1.0e-8) + res2 = m3.orderConsecutiveCells1D() + expRes2 = [0,1,2] + self.assertEqual(m3.getNumberOfCells(),res2.getNumberOfTuples()) + self.assertEqual(expRes2, res2.getValues()) + pass + pass if __name__ == '__main__': diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 30bcf64cc..7e5a3b181 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -295,6 +295,7 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingUMesh::ComputeRangesFromTypeDistribution; %newobject ParaMEDMEM::MEDCouplingUMesh::buildUnionOf2DMesh; %newobject ParaMEDMEM::MEDCouplingUMesh::buildUnionOf3DMesh; +%newobject ParaMEDMEM::MEDCouplingUMesh::orderConsecutiveCells1D; %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTreeFast; %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic; %newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic; @@ -1658,6 +1659,7 @@ namespace ParaMEDMEM DataArrayInt *convertNodalConnectivityToStaticGeoTypeMesh() const throw(INTERP_KERNEL::Exception); DataArrayInt *buildUnionOf2DMesh() const throw(INTERP_KERNEL::Exception); DataArrayInt *buildUnionOf3DMesh() const throw(INTERP_KERNEL::Exception); + DataArrayInt *orderConsecutiveCells1D() const throw(INTERP_KERNEL::Exception); DataArrayDouble *getBoundingBoxForBBTreeFast() const throw(INTERP_KERNEL::Exception); DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception); DataArrayDouble *getBoundingBoxForBBTree1DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception); -- 2.39.2