- * 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<int>& edge1Bis, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& edge1BisPtr,
- std::vector< std::vector<int> >& out0, std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > >& out1)
-{
- 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<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
- for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
- //
- 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<int>& connOut(out0[0]);
- connOut.resize(nbOfEdgesOf2DCellSplit);
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
- edgesPtr.resize(nbOfEdgesOf2DCellSplit);
- for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
- {
- connOut[kk]=edge1Bis[2*kk];
- edgesPtr[kk]=edge1BisPtr[2*kk];
- }
- }
- else
- {
- // [i,iEnd[ contains the
- out0.resize(2); out1.resize(2);
- std::vector<int>& connOutLeft(out0[0]);
- std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& eleft(out1[0]);
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& eright(out1[1]);
- for(std::size_t k=ii;k<jj+1;k++)
- { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > ees(iEnd);
- for(int ik=0;ik<iEnd;ik++)
- {
- std::map< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int> m;
- MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> 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--)
- connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
- for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
- { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
- eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
- for(int ik=0;ik<iEnd;ik++)
- connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
- eright.insert(eright.end(),ees.begin(),ees.end());
- }
-}
-
-/// @endcond
-
-/// @cond INTERNAL
-
-struct CellInfo
-{
-public:
- CellInfo() { }
- CellInfo(const std::vector<int>& edges, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& edgesPtr);
-public:
- std::vector<int> _edges;
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > _edges_ptr;
-};
-
-CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& edgesPtr)
-{
- std::size_t nbe(edges.size());
- std::vector<int> edges2(2*nbe); std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
- for(std::size_t i=0;i<nbe;i++)
- {
- edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
- edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
- }
- _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
- std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
- std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
-}
-
-class EdgeInfo
-{
-public:
- EdgeInfo(int istart, int iend, const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
- EdgeInfo(int istart, int iend, int pos, const MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge>& 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<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& newRight);
- void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
-private:
- int _istart;
- int _iend;
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> _mesh;
- MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> _edge;
- int _left;
- int _right;
-};
-
-void EdgeInfo::somethingHappendAt(int pos, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& newRight)
-{
- const MEDCouplingUMesh *mesh(_mesh);
- if(mesh)
- return ;
- if(_right<pos)
- return ;
- if(_left>pos)
- { _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
- {
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> 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;
- }
- }
-}
-
-class VectorOfCellInfo
-{
-public:
- VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& edgesPtr);
- std::size_t size() const { return _pool.size(); }
- int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
- void setMeshAt(int pos, const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh>& mesh, int istart, int iend, const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > >& edgePtrs);
- const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
- const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> 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<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& newRight);
- const CellInfo& get(int pos) const;
- CellInfo& get(int pos);
-private:
- std::vector<CellInfo> _pool;
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> _ze_mesh;
- std::vector<EdgeInfo> _edge_info;
-};
-
-VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& 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<DataArrayDouble> barys(mesh->getBarycenterAndOwner());
- return zeMesh->getCellContainingPoint(barys->begin(),eps);
-}
-
-void VectorOfCellInfo::setMeshAt(int pos, const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh>& mesh, int istart, int iend, const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > >& 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<CellInfo> pool(_pool.size()-1+sz);
- for(int i=0;i<pos;i++)
- pool[i]=_pool[i];
- for(std::size_t j=0;j<sz;j++)
- pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
- for(int i=pos+1;i<(int)_pool.size();i++)
- pool[i+sz-1]=_pool[i];
- _pool=pool;
- //
- if(sz==2)
- updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
- //
- if(isFast)
- {
- _ze_mesh=mesh;
- return ;
- }
- //
- std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > ms;
- if(pos>0)
- {
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelf2(0,pos,true)));
- ms.push_back(elt);
- }
- ms.push_back(mesh);
- if(pos<_ze_mesh->getNumberOfCells()-1)
- {
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_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<ms2.size();j++)
- ms2[j]=ms[j];
- _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
-}
-
-void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
-{
- _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
-}
-
-int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
-{
- if(pos<0)
- throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
- int ret(0);
- for(std::vector<EdgeInfo>::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<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& 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<sz;i++)
- _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
-}
-
-const CellInfo& VectorOfCellInfo::get(int pos) const
-{
- if(pos<0 || pos>=(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];
-}
-
-/*!
- * Given :
- * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
- * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
- *
- * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
- *
- * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
- *
- * \param [in] allEdges a list of pairs (beginNode, endNode). Linked with \a allEdgesPtr to get the equation of edge.
- */
-MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt>& 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::vector<int> edge1Bis(nb*2);
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > 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());
- VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
- for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
- {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
- int iEnd(iStart);
- for(;iEnd<nbCellsInSplitMesh1D;)
- {
- for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
- if(jj!=nb)
- break;
- else
- iEnd++;
- }
- if(iEnd<nbCellsInSplitMesh1D)
- iEnd++;
- //
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelf2(iStart,iEnd,1,true)));
- int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
- //
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
- retTmp->setCoords(splitMesh1D->getCoords());
- retTmp->allocateCells();
-
- std::vector< std::vector<int> > out0;
- std::vector< std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > > out1;
-
- BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
- for(std::size_t cnt=0;cnt<out0.size();cnt++)
- AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
- pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
- //
- iStart=iEnd;
- }
- for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
- pool.feedEdgeInfoAt(eps,mm,offset,idsLeftRightPtr+2*mm);
- return pool.getZeMesh().retn();
-}
-
-MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, const MEDCouplingUMesh *splitMesh1D,
- const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt>& idsLeftRight)
-{
- const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
- //
- std::vector<int> allEdges;
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> > 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
- {
- int edgeId(std::abs(*it)-1);
- std::map< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int> m;
- MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
- const std::vector<int>& 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<sz;cnt++)
- allEdgesPtr.push_back(ee);
- }
- //
- 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)
-{
- if(!typ1.isQuadratic() && !typ2.isQuadratic())
- {//easy case comparison not
- return conn1[0]==conn2[0] && conn1[1]==conn2[1];
- }
- else if(typ1.isQuadratic() && typ2.isQuadratic())
- {
- bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
- if(!status0)
- return false;
- if(conn1[2]==conn2[2])
- return true;
- const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
- double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
- return dist<eps;
- }
- else
- {//only one is quadratic
- bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
- if(!status0)
- return false;
- const double *a(0),*bb(0),*be(0);
- if(typ1.isQuadratic())
- {
- a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
- }
- else
- {
- a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
- }
- double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
- double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
- return dist<eps;
- }
-}
-
-/*!
- * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
- * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
- *
- * \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)
-{
- 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);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
- if(cellIdInMesh1DSplitRelative<0)
- cur1D->changeOrientationOfCells();
- const int *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++)
- {
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
- const int *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<int> tmpPtr(new int[ci[1]-ci[0]]);
- for(unsigned it2=0;it2<sz;it2++)
- {
- INTERP_KERNEL::NormalizedCellType typeOfSon;
- cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
- const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
- if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
- return *it;
- }
- }
- throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
-}
-
-/// @endcond
-
-/*!
- * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
- * Thus the final result contains the aggregation of nodes of \a mesh2D, then nodes of \a mesh1D, then new nodes that are the result of the intersection
- * and finaly, 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.
- *
- * \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
- * you can invoke orderConsecutiveCells1D on \a mesh1D.
- * \param [in] eps - precision used to perform intersections and localization operations.
- * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
- * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
- * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
- * So this array has a number of tuples equal to the number of cells of \a splitMesh2D and a number of component equal to 1.
- * \param [out] cellIdInMesh1D - the array of pair that gives for each cell id \a i in \a splitMesh1D the cell in \a splitMesh2D on the left for the 1st component
- * and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
- * So this array has a number of tuples equal to the number of cells of \a splitMesh1D and a number of components equal to 2.
- *
- * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
- */
-void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
-{
- if(!mesh2D || !mesh1D)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
- mesh2D->checkFullyDefined();
- mesh1D->checkFullyDefined();
- const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
- 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<int> > intersectEdge1, colinear2, subDiv2;
- std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
- //
- // Build desc connectivity
- DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
- std::map<int,int> mergedNodes;
- Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
- // use mergeNodes to fix intersectEdge1
- for(std::vector< std::vector<int> >::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<int,int>::const_iterator it1;
- it1=mergedNodes.find(eltStart);
- if(it1!=mergedNodes.end())
- (*it0)[0]=(*it1).second;
- it1=mergedNodes.find(eltEnd);
- if(it1!=mergedNodes.end())
- (*it0)[2*n-1]=(*it1).second;
- }
- //
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCooDa(DataArrayDouble::New());
- addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
- // Step 2: re-order newly created nodes according to the ordering found in m2
- std::vector< std::vector<int> > intersectEdge2;
- BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
- subDiv2.clear();
- // Step 3: compute splitMesh1D
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
- idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
- // deal with cells in mesh2D that are not cut but only some of their edges are
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCpy());
- idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
- idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> 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);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> outi0s(outi0);
- out0s=out0;
- out0s=out0s->buildUnique();
- out0s->sort(true);
- }
- //
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> baryRet1(ret1NonCol->getBarycenterAndOwner());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> elts,eltsIndex;
- mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsIndex2(eltsIndex->deltaShiftIndex());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsIndex3(eltsIndex2->getIdsEqual(1));
- if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
- throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellsToBeModified(elts->buildUnique());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
- if((DataArrayInt *)out0s)
- untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
- std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > outMesh2DSplit;
- // OK all is ready to insert in ret2 mesh
- if(!untouchedCells->empty())
- {// the most easy part, cells in mesh2D not impacted at all
- outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
- outMesh2DSplit.back()->setCoords(ret1->getCoords());
- ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
- }
- if((DataArrayInt *)out0s)
- {// here dealing with cells in out0s but not in cellsToBeModified
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> 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++)
- {
- 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());
- ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
- partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
- int kk(0),*ret3ptr(partOfRet3->getPointer());
- for(const int *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++)
- {
- int tmp(fewModifiedCells->locateValue(*it2));
- if(tmp!=-1)
- {
- if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
- ret3ptr[2*kk]=tmp+offset;
- if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
- ret3ptr[2*kk+1]=tmp+offset;
- }
- else
- {//the current edge is shared by a 2D cell that will be split just after
- if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
- ret3ptr[2*kk]=-(*it2+1);
- if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
- ret3ptr[2*kk+1]=-(*it2+1);
- }
- }
- }
- m1Desc->setCoords(ret1->getCoords());
- ret1NonCol->setCoords(ret1->getCoords());
- ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
- if(!outMesh2DSplit.empty())
- {
- DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
- for(std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
- (*itt)->setCoords(da);
- }
- }
- cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
- for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
- {
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> idsNonColPerCell(elts->getIdsEqual(*it));
- idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> partOfRet3;
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> 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;i<splitOfOneCell->getNumberOfCells();i++)
- ret2->pushBackSilent(*it);
- }
- //
- std::size_t nbOfMeshes(outMesh2DSplit.size());
- std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
- for(std::size_t i=0;i<nbOfMeshes;i++)
- tmp[i]=outMesh2DSplit[i];
- //
- ret1->getCoords()->setInfoOnComponents(compNames);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
- // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
- ret3->rearrange(1);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> edgesToDealWith(ret3->getIdsStrictlyNegative());
- for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
- {
- int old2DCellId(-ret3->getIJ(*it,0)-1);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> candidates(ret2->getIdsEqual(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<int>::max(),-1);
- ret3->rearrange(2);
- //
- splitMesh1D=ret1.retn();
- splitMesh2D=ret2D.retn();
- cellIdInMesh2D=ret2.retn();
- cellIdInMesh1D=ret3.retn();
-}
-
-/**
- * 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 correspondance 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<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
- const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
- const std::vector<double>& addCoords,
- std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
-{
- static const int SPACEDIM=2;
- const double *coo1(m1->getCoords()->getConstPointer());
- const int *conn1(m1->getNodalConnectivity()->getConstPointer()),*connI1(m1->getNodalConnectivityIndex()->getConstPointer());
- int offset1(m1->getNumberOfNodes());
- const double *coo2(m2->getCoords()->getConstPointer());
- const int *conn2(m2->getNodalConnectivity()->getConstPointer()),*connI2(m2->getNodalConnectivityIndex()->getConstPointer());
- int offset2(offset1+m2->getNumberOfNodes());
- int offset3(offset2+((int)addCoords.size())/2);
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
- const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
- // Here a BBTree on 2D-cells, not on segments:
- BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
- int ncell1(m1->getNumberOfCells());
- crI.push_back(0);
- for(int i=0;i<ncell1;i++)
- {
- std::vector<int> candidates2;
- myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
- std::map<INTERP_KERNEL::Node *,int> mapp;
- std::map<int,INTERP_KERNEL::Node *> 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.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
- desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
- //
- std::set<INTERP_KERNEL::Edge *> edges1;// store all edges of pol1 that are NOT consumed by intersect cells. If any after iteration over candidates2 -> a part of pol1 should appear in result
- std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
- INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
- for(it1.first();!it1.finished();it1.next())
- edges1.insert(it1.current()->getPtr());
- //
- std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
- std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
- int ii=0;
- for(std::vector<int>::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);
- // Complete mapping with elements coming from the current cell it2 in mesh2:
- MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
- // pol2 is the new QP in the final merged result.
- pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
- pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
- }
- ii=0;
- for(std::vector<int>::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);
- //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
- pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
- }
- // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
- // by m2 but that we still want to keep in the final result.
- if(!edges1.empty())
- {
- try
- {
- INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
- }
- catch(INTERP_KERNEL::Exception& e)
- {
- std::ostringstream oss; oss << "Error when computing residual of cell #" << i << " in source/m1 mesh ! Maybe the neighbours of this cell in mesh are not well connected !\n" << "The deep reason is the following : " << e.what();
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
- for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
- (*it).second->decrRef();
- }
-}
-
-/**
- * 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. This method only works using nodal connectivity consideration. Coordinates of nodes are ignored here.
- * 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)
- *
- * \sa DataArrayInt::sortEachPairToMakeALinkedList
- */
-DataArrayInt *MEDCouplingUMesh::orderConsecutiveCells1D() const
-{
- checkFullyDefined();
- if(getMeshDimension()!=1)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::orderConsecutiveCells1D works on unstructured mesh with meshdim = 1 !");
-
- // Check that this is a line (and not a more complex 1D mesh) - each point is used at most by 2 segments:
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> _d(DataArrayInt::New()),_dI(DataArrayInt::New());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> _rD(DataArrayInt::New()),_rDI(DataArrayInt::New());
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m_points(buildDescendingConnectivity(_d, _dI, _rD, _rDI));
- const int *d(_d->getConstPointer()), *dI(_dI->getConstPointer());
- const int *rD(_rD->getConstPointer()), *rDI(_rDI->getConstPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> _dsi(_rDI->deltaShiftIndex());
- const int * dsi(_dsi->getConstPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> 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<DataArrayInt> result(DataArrayInt::New());
- result->alloc(nc,1);
-
- // set of edges not used so far
- std::set<int> edgeSet;
- for (int i=0; i<nc; edgeSet.insert(i), i++);
-
- int startSeg=0;
- int newIdx=0;
- // while we have points with only one neighbor segments
- do
- {
- std::list<int> 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<MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
-{
- MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node> nTmp(n); nTmp->incrRef();
- std::map<MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int>::const_iterator it(m.find(nTmp));
- if(it==m.end())
- throw INTERP_KERNEL::Exception("Internal error in remapping !");
- int 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<MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
-{
- int sz(c.size());
- if(sz<=1)
- return false;
- bool presenceOfOn(false);
- for(int i=0;i<sz;i++)
- {
- INTERP_KERNEL::ElementaryEdge *e(c[i]);
- if(e->getLoc()!=INTERP_KERNEL::FULL_ON_1)
- continue ;
- IKGeo2DInternalMapper2(e->getStartNode(),m,forbVal0,forbVal1,isect);
- IKGeo2DInternalMapper2(e->getEndNode(),m,forbVal0,forbVal1,isect);
- }
- 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 \ref numbering-indirect 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.
- * So this method try to reduce at most the number of new nodes. The only case that can lead this method to add
- * nodes if a SEG3 is split without information of middle.
- * \b WARNING : is returned value is different from 0 a call to MEDCouplingUMesh::mergeNodes is necessary to
- * avoid to have a non conform mesh.
- *
- * \return int - the number of new nodes created (in most of cases 0).
- *
- * \throw If \a this is not coherent.
- * \throw If \a this has not spaceDim equal to 2.
- * \throw If \a this has not meshDim equal to 2.
- * \throw If some subcells needed to be split are orphan.
- * \sa MEDCouplingUMesh::conformize2D
- */
-int MEDCouplingUMesh::split2DCells(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *midOpt, const DataArrayInt *midOptI)
-{
- if(!desc || !descI || !subNodesInSeg || !subNodesInSegI)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCells : the 4 first arrays must be not null !");
- desc->checkAllocated(); descI->checkAllocated(); subNodesInSeg->checkAllocated(); subNodesInSegI->checkAllocated();
- if(getSpaceDimension()!=2 || getMeshDimension()!=2)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCells : This method only works for meshes with spaceDim=2 and meshDim=2 !");
- if(midOpt==0 && midOptI==0)
- {
- split2DCellsLinear(desc,descI,subNodesInSeg,subNodesInSegI);
- return 0;
- }
- else if(midOpt!=0 && midOptI!=0)
- return split2DCellsQuadratic(desc,descI,subNodesInSeg,subNodesInSegI,midOpt,midOptI);
- else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCells : middle parameters must be set to null for all or not null for all.");
-}
-
-/*!
- * \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
- * 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 invokation of this method.
- *
- * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
- * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
- *
- * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
- * 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
- * that the user is expected to deal with.
- *
- * \throw If \a this is not coherent.
- * \throw If \a this has not spaceDim equal to 2.
- * \throw If \a this has not meshDim equal to 2.
- * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
- */
-DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
-{
- static const int SPACEDIM=2;
- checkCoherency();
- if(getSpaceDimension()!=2 || getMeshDimension()!=2)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
- const int *c(mDesc->getNodalConnectivity()->getConstPointer()),*ci(mDesc->getNodalConnectivityIndex()->getConstPointer()),*rd(revDesc1->getConstPointer()),*rdi(revDescIndx1->getConstPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
- const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
- int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
- std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
- std::vector<double> addCoo;
- BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
- for(int i=0;i<nDescCell;i++)
- {
- std::vector<int> candidates;
- myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
- for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
- if(*it>i)
- {
- std::map<MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int> 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;
- INTERP_KERNEL::QuadraticPolygon c1,c2;
- e1->intersectWith(e2,merge,c1,c2);
- e1->decrRef(); e2->decrRef();
- if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
- overlapEdge[i].push_back(*it);
- if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
- overlapEdge[*it].push_back(i);
- }
- }
- // splitting done. sort intersect point in intersectEdge.
- std::vector< std::vector<int> > middle(nDescCell);
- int nbOf2DCellsToBeSplit(0);
- bool middleNeedsToBeUsed(false);
- std::vector<bool> cells2DToTreat(nDescCell,false);
- for(int i=0;i<nDescCell;i++)
- {
- std::vector<int>& isect(intersectEdge[i]);
- int sz((int)isect.size());
- if(sz>1)
- {
- std::map<MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node>,int> 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]);
- if(idx1-idx0!=1)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
- if(!cells2DToTreat[rd[idx0]])
- {
- cells2DToTreat[rd[idx0]]=true;
- nbOf2DCellsToBeSplit++;
- }
- // try to reuse at most eventual 'middle' of SEG3
- std::vector<int>& mid(middle[i]);
- mid.resize(sz+1,-1);
- if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
- {
- middleNeedsToBeUsed=true;
- const std::vector<int>& candidates(overlapEdge[i]);
- std::vector<int> trueCandidates;
- for(std::vector<int>::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<sz+1;j++)
- {
- for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
- {
- int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
- if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
- { mid[j]=*itc; break; }
- }
- stNode=endNode;
- endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
- }
- }
- }
- }
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
- if(nbOf2DCellsToBeSplit==0)
- return ret.retn();
- //
- int *retPtr(ret->getPointer());
- for(int i=0;i<nCell;i++)
- if(cells2DToTreat[i])
- *retPtr++=i;
- //
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> 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;
- if(middleNeedsToBeUsed)
- { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
- int 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;
- if(nbOfNodesCreated!=0)
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
- }
- return ret.retn();
-}
-
-/*!
- * This non const method works on 2D mesh. This method scans every cell in \a this and look if each edge constituting this cell is not mergeable with neighbors edges of that cell.
- * If yes, the cell is "repaired" to minimize at most its number of edges. So this method do not change the overall shape of cells in \a this (with eps precision).
- * This method do not take care of shared edges between cells, so this method can lead to a non conform mesh (\a this). If a conform mesh is required you're expected
- * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
- * This method works on any 2D geometric types of cell (even static one). If a cell is touched its type becomes dynamic automaticaly. For 2D "repaired" quadratic cells
- * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
- *
- * If the returned array is empty it means that nothing has changed in \a this (as if it were a const method). If the array is not empty the connectivity of \a this is modified
- * using new instance, idem for coordinates.
- *
- * 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.
- *
- * \throw If \a this is not coherent.
- * \throw If \a this has not spaceDim equal to 2.
- * \throw If \a this has not meshDim equal to 2.
- *
- * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
- */
-DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
-{
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
- checkCoherency();
- 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::QUADRATIC_PLANAR::_arc_detection_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
- const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> 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;i<nbOfCells;i++,newciptr++,ciptr++)
- {
- if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,newc,appendedCoords))
- ret->pushBackSilent(i);
- newciptr[1]=newc->getNumberOfTuples();
- }
- //
- if(ret->empty())
- return ret.retn();
- if(!appendedCoords->empty())
- {
- appendedCoords->rearrange(2);
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
- //non const part
- setCoords(newCoords);
- }
- //non const part
- setConnectivity(newc,newci,true);
- return ret.retn();
-}
-
-/*!
- * \param [out] intersectEdge1 - for each cell in \a m1Desc returns the result of the split. The result is given using pair of int given resp start and stop.
- * So for all edge \a i in \a m1Desc \a intersectEdge1[i] is of length 2*n where n is the number of sub edges.
- * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
- * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
- * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
- * \param [out] addCoo - nodes to be append at the end
- * \param [out] mergedNodes - gives all pair of nodes of \a m2Desc that have same location than some nodes in \a m1Desc. key is id in \a m2Desc offseted and value is id in \a m1Desc.
- */
-void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
- std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2, std::vector<double>& addCoo, std::map<int,int>& mergedNodes)
-{
- static const int SPACEDIM=2;
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
- const int *c1(m1Desc->getNodalConnectivity()->getConstPointer()),*ci1(m1Desc->getNodalConnectivityIndex()->getConstPointer());
- // Build BB tree of all edges in the tool mesh (second mesh)
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2Desc->getBoundingBoxForBBTree());
- const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
- int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
- intersectEdge1.resize(nDescCell1);
- colinear2.resize(nDescCell2);
- subDiv2.resize(nDescCell2);
- BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
-
- std::vector<int> candidates1(1);
- int offset1(m1Desc->getNumberOfNodes());
- int offset2(offset1+m2Desc->getNumberOfNodes());
- for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
- {
- std::vector<int> 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<INTERP_KERNEL::Node *,int> map1,map2;
- // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
- INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
- candidates1[0]=i;
- INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
- // This following part is to avoid that some removed nodes (for example due to a merge between pol1 and pol2) are replaced by a newly created one
- // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
- std::set<INTERP_KERNEL::Node *> nodes;
- pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
- std::size_t szz(nodes.size());
- std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node> > nodesSafe(szz);
- std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
- for(std::size_t iii=0;iii<szz;iii++,itt++)
- { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
- // end of protection
- // Performs egde cutting:
- pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
- delete pol2;
- delete pol1;
- }
- else
- // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
- intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
- }
-}
-
-/*!
- * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
- * It builds the descending connectivity of the two meshes, and then using a binary tree
- * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
- * 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<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
- MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
- std::vector<double>& addCoo,
- MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&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();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
- m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
- m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
- std::map<int,int> 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();
-}
-
-/*!
- * This method performs the 2nd step of Partition of 2D mesh.
- * This method has 4 inputs :
- * - a mesh 'm1' with meshDim==1 and a SpaceDim==2
- * - a mesh 'm2' with meshDim==1 and a SpaceDim==2
- * - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids randomly sorted.
- * The aim of this method is to sort the splitting nodes, if any, and to put them in 'intersectEdge' output parameter based on edges of mesh 'm2'
- * Nodes end up lying consecutively on a cutted edge.
- * \param m1 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method.
- * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
- * \param m2 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method.
- * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
- * \param[out] intersectEdge the same content as subDiv, but correclty oriented.
- */
-void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
- const std::vector<double>& addCoo,
- const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
-{
- int offset1=m1->getNumberOfNodes();
- int ncell=m2->getNumberOfCells();
- const int *c=m2->getNodalConnectivity()->getConstPointer();
- const int *cI=m2->getNodalConnectivityIndex()->getConstPointer();
- const double *coo=m2->getCoords()->getConstPointer();
- const double *cooBis=m1->getCoords()->getConstPointer();
- int offset2=offset1+m2->getNumberOfNodes();
- intersectEdge.resize(ncell);
- for(int i=0;i<ncell;i++,cI++)
- {
- const std::vector<int>& divs=subDiv[i];
- int nnode=cI[1]-cI[0]-1;
- std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;
- std::map<INTERP_KERNEL::Node *, int> mapp22;
- for(int j=0;j<nnode;j++)
- {
- INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
- int nnid=c[(*cI)+j+1];
- mapp2[nnid]=std::pair<INTERP_KERNEL::Node *,bool>(nn,true);
- mapp22[nn]=nnid+offset1;
- }
- INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
- for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
- ((*it).second.first)->decrRef();
- std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
- std::map<INTERP_KERNEL::Node *,int> mapp3;
- for(std::size_t j=0;j<divs.size();j++)
- {
- int id=divs[j];
- INTERP_KERNEL::Node *tmp=0;
- if(id<offset1)
- tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
- else if(id<offset2)
- tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
- else
- tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
- addNodes[j]=tmp;
- mapp3[tmp]=id;
- }
- e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
- for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
- (*it)->decrRef();
- e->decrRef();
- }
-}
-
-/*!
- * This method is part of the Slice3D algorithm. It is the first step of assembly process, ones coordinates have been computed (by MEDCouplingUMesh::split3DCurveWithPlane method).
- * This method allows to compute given the status of 3D curve cells and the descending connectivity 3DSurf->3DCurve to deduce the intersection of each 3D surf cells
- * with a plane. The result will be put in 'cut3DSuf' out parameter.
- * \param [in] cut3DCurve input paramter that gives for each 3DCurve cell if it owns fully to the plane or partially.
- * \param [out] nodesOnPlane, returns all the nodes that are on the plane.
- * \param [in] nodal3DSurf is the nodal connectivity of 3D surf mesh.
- * \param [in] nodalIndx3DSurf is the nodal connectivity index of 3D surf mesh.
- * \param [in] nodal3DCurve is the nodal connectivity of 3D curve mesh.
- * \param [in] nodal3DIndxCurve is the nodal connectivity index of 3D curve mesh.
- * \param [in] desc is the descending connectivity 3DSurf->3DCurve
- * \param [in] descIndx is the descending connectivity index 3DSurf->3DCurve
- * \param [out] cut3DSuf input/output param.
- */
-void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<int>& cut3DCurve, std::vector<int>& nodesOnPlane, const int *nodal3DSurf, const int *nodalIndx3DSurf,
- const int *nodal3DCurve, const int *nodalIndx3DCurve,
- const int *desc, const int *descIndx,
- std::vector< std::pair<int,int> >& cut3DSurf)
-{
- std::set<int> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
- int nbOf3DSurfCell=(int)cut3DSurf.size();
- for(int i=0;i<nbOf3DSurfCell;i++)
- {
- std::vector<int> res;
- int offset=descIndx[i];
- int nbOfSeg=descIndx[i+1]-offset;
- for(int j=0;j<nbOfSeg;j++)
- {
- int edgeId=desc[offset+j];
- int status=cut3DCurve[edgeId];
- if(status!=-2)
- {
- if(status>-1)
- res.push_back(status);
- else
- {
- res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
- res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
- }
- }
- }
- switch(res.size())
- {
- case 2:
- {
- cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
- break;
- }
- case 1:
- case 0:
- {
- std::set<int> s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]);
- std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector<int> >(res));
- if(res.size()==2)
- {
- cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
- }
- else
- {
- cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
- }
- break;
- }
- default:
- {// case when plane is on a multi colinear edge of a polyhedron
- if((int)res.size()==2*nbOfSeg)
- {
- cut3DSurf[i].first=-2; cut3DSurf[i].second=i;
- }
- else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
- }
- }
- }
-}
-
-/*!
- * \a this is expected to be a mesh with spaceDim==3 and meshDim==3. If not an exception will be thrown.
- * This method is part of the Slice3D algorithm. It is the second step of assembly process, ones coordinates have been computed (by MEDCouplingUMesh::split3DCurveWithPlane method).
- * This method allows to compute given the result of 3D surf cells with plane and the descending connectivity 3D->3DSurf to deduce the intersection of each 3D cells
- * with a plane. The result will be put in 'nodalRes' 'nodalResIndx' and 'cellIds' out parameters.
- * \param cut3DSurf input paramter that gives for each 3DSurf its intersection with plane (result of MEDCouplingUMesh::AssemblyForSplitFrom3DCurve).
- * \param desc is the descending connectivity 3D->3DSurf
- * \param descIndx is the descending connectivity index 3D->3DSurf
- */
-void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<int,int> >& cut3DSurf,
- const int *desc, const int *descIndx,
- DataArrayInt *nodalRes, DataArrayInt *nodalResIndx, DataArrayInt *cellIds) const
-{
- checkFullyDefined();
- if(getMeshDimension()!=3 || getSpaceDimension()!=3)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
- const int *nodal3D=_nodal_connec->getConstPointer();
- const int *nodalIndx3D=_nodal_connec_index->getConstPointer();
- int nbOfCells=getNumberOfCells();
- for(int i=0;i<nbOfCells;i++)
- {
- std::map<int, std::set<int> > m;
- int offset=descIndx[i];
- int nbOfFaces=descIndx[i+1]-offset;
- int start=-1;
- int end=-1;
- for(int j=0;j<nbOfFaces;j++)
- {
- const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
- if(p.first!=-1 && p.second!=-1)
- {
- if(p.first!=-2)
- {
- start=p.first; end=p.second;
- m[p.first].insert(p.second);
- m[p.second].insert(p.first);
- }
- else
- {
- 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<int> tmp=new int[sz];
- INTERP_KERNEL::NormalizedCellType cmsId;
- unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId);
- start=tmp[0]; end=tmp[nbOfNodesSon-1];
- for(unsigned k=0;k<nbOfNodesSon;k++)
- {
- m[tmp[k]].insert(tmp[(k+1)%nbOfNodesSon]);
- m[tmp[(k+1)%nbOfNodesSon]].insert(tmp[k]);
- }
- }
- }
- }
- if(m.empty())
- continue;
- std::vector<int> conn(1,(int)INTERP_KERNEL::NORM_POLYGON);
- int prev=end;
- while(end!=start)
- {
- std::map<int, std::set<int> >::const_iterator it=m.find(start);
- const std::set<int>& s=(*it).second;
- std::set<int> s2; s2.insert(prev);
- std::set<int> s3;
- std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
- if(s3.size()==1)
- {
- int val=*s3.begin();
- conn.push_back(start);
- prev=start;
- start=val;
- }
- else
- start=end;
- }
- conn.push_back(end);
- if(conn.size()>3)
- {
- nodalRes->insertAtTheEnd(conn.begin(),conn.end());
- nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples());
- cellIds->pushBackSilent(i);
- }
- }
-}
-
-/*!
- * This method compute the convex hull of a single 2D cell. This method tries to conserve at maximum the given input connectivity. In particular, if the orientation of cell is not clockwise
- * as in MED format norm. If definitely the result of Jarvis algorithm is not matchable with the input connectivity, the result will be copied into \b nodalConnecOut parameter and
- * the geometric cell type set to INTERP_KERNEL::NORM_POLYGON.
- * This method excepts that \b coords parameter is expected to be in dimension 2. [ \b nodalConnBg , \b nodalConnEnd ) is the nodal connectivity of the input
- * cell (geometric cell type included at the position 0). If the meshdimension of the input cell is not equal to 2 an INTERP_KERNEL::Exception will be thrown.
- *
- * \return false if the input connectivity represents already the convex hull, true if the input cell needs to be reordered.
- */
-bool MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis(const double *coords, const int *nodalConnBg, const int *nodalConnEnd, DataArrayInt *nodalConnecOut)
-{
- std::size_t sz=std::distance(nodalConnBg,nodalConnEnd);
- if(sz>=4)
- {
- const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)*nodalConnBg);
- if(cm.getDimension()==2)
- {
- const int *node=nodalConnBg+1;
- int startNode=*node++;
- double refX=coords[2*startNode];
- for(;node!=nodalConnEnd;node++)
- {
- if(coords[2*(*node)]<refX)
- {
- startNode=*node;
- refX=coords[2*startNode];
- }
- }
- std::vector<int> tmpOut; tmpOut.reserve(sz); tmpOut.push_back(startNode);
- refX=1e300;
- double tmp1;
- double tmp2[2];
- double angle0=-M_PI/2;
- //
- int nextNode=-1;
- int prevNode=-1;
- double resRef;
- double angleNext=0.;
- while(nextNode!=startNode)
- {
- nextNode=-1;
- resRef=1e300;
- for(node=nodalConnBg+1;node!=nodalConnEnd;node++)
- {
- if(*node!=tmpOut.back() && *node!=prevNode)
- {
- tmp2[0]=coords[2*(*node)]-coords[2*tmpOut.back()]; tmp2[1]=coords[2*(*node)+1]-coords[2*tmpOut.back()+1];
- double angleM=INTERP_KERNEL::EdgeArcCircle::GetAbsoluteAngle(tmp2,tmp1);
- double res;
- if(angleM<=angle0)
- res=angle0-angleM;
- else
- res=angle0-angleM+2.*M_PI;
- if(res<resRef)
- {
- nextNode=*node;
- resRef=res;
- angleNext=angleM;
- }
- }
- }
- if(nextNode!=startNode)
- {
- angle0=angleNext-M_PI;
- if(angle0<-M_PI)
- angle0+=2*M_PI;
- prevNode=tmpOut.back();
- tmpOut.push_back(nextNode);
- }
- }
- std::vector<int> tmp3(2*(sz-1));
- std::vector<int>::iterator it=std::copy(nodalConnBg+1,nodalConnEnd,tmp3.begin());
- std::copy(nodalConnBg+1,nodalConnEnd,it);
- if(std::search(tmp3.begin(),tmp3.end(),tmpOut.begin(),tmpOut.end())!=tmp3.end())
- {
- nodalConnecOut->insertAtTheEnd(nodalConnBg,nodalConnEnd);
- return false;
- }
- if(std::search(tmp3.rbegin(),tmp3.rend(),tmpOut.begin(),tmpOut.end())!=tmp3.rend())
- {
- nodalConnecOut->insertAtTheEnd(nodalConnBg,nodalConnEnd);
- return false;
- }
- else
- {
- nodalConnecOut->pushBackSilent((int)INTERP_KERNEL::NORM_POLYGON);
- nodalConnecOut->insertAtTheEnd(tmpOut.begin(),tmpOut.end());
- return true;
- }
- }
- else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis : invalid 2D cell connectivity !");
- }
- else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis : invalid 2D cell connectivity !");
-}
-
-/*!
- * This method works on an input pair (\b arr, \b arrIndx) where \b arr indexes is in \b arrIndx.
- * This method will not impact the size of inout parameter \b arrIndx but the size of \b arr will be modified in case of suppression.
- *
- * \param [in] idsToRemoveBg begin of set of ids to remove in \b arr (included)
- * \param [in] idsToRemoveEnd end of set of ids to remove in \b arr (excluded)
- * \param [in,out] arr array in which the remove operation will be done.
- * \param [in,out] arrIndx array in the remove operation will modify
- * \param [in] offsetForRemoval (by default 0) offset so that for each i in [0,arrIndx->getNumberOfTuples()-1) removal process will be performed in the following range [arr+arrIndx[i]+offsetForRemoval,arr+arr[i+1])
- * \return true if \b arr and \b arrIndx have been modified, false if not.
- */
-bool MEDCouplingUMesh::RemoveIdsFromIndexedArrays(const int *idsToRemoveBg, const int *idsToRemoveEnd, DataArrayInt *arr, DataArrayInt *arrIndx, int offsetForRemoval)
-{
- if(!arrIndx || !arr)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::RemoveIdsFromIndexedArrays : some input arrays are empty !");
- if(offsetForRemoval<0)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::RemoveIdsFromIndexedArrays : offsetForRemoval should be >=0 !");
- std::set<int> s(idsToRemoveBg,idsToRemoveEnd);
- int nbOfGrps=arrIndx->getNumberOfTuples()-1;
- int *arrIPtr=arrIndx->getPointer();
- *arrIPtr++=0;
- int previousArrI=0;
- const int *arrPtr=arr->getConstPointer();
- std::vector<int> arrOut;//no utility to switch to DataArrayInt because copy always needed
- for(int i=0;i<nbOfGrps;i++,arrIPtr++)
- {
- if(*arrIPtr-previousArrI>offsetForRemoval)
- {
- for(const int *work=arrPtr+previousArrI+offsetForRemoval;work!=arrPtr+*arrIPtr;work++)
- {
- if(s.find(*work)==s.end())
- arrOut.push_back(*work);
- }
- }
- previousArrI=*arrIPtr;
- *arrIPtr=(int)arrOut.size();
- }
- if(arr->getNumberOfTuples()==(int)arrOut.size())
- return false;
- arr->alloc((int)arrOut.size(),1);
- std::copy(arrOut.begin(),arrOut.end(),arr->getPointer());
- return true;
-}
-
-/*!
- * This method works on a pair input (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn
- * (\ref numbering-indirect).
- * This method returns the result of the extraction ( specified by a set of ids in [\b idsOfSelectBg , \b idsOfSelectEnd ) ).
- * The selection of extraction is done standardly in new2old format.
- * This method returns indexed arrays (\ref numbering-indirect) using 2 arrays (arrOut,arrIndexOut).
- *
- * \param [in] idsOfSelectBg begin of set of ids of the input extraction (included)
- * \param [in] idsOfSelectEnd end of set of ids of the input extraction (excluded)
- * \param [in] arrIn arr origin array from which the extraction will be done.
- * \param [in] arrIndxIn is the input index array allowing to walk into \b arrIn
- * \param [out] arrOut the resulting array
- * \param [out] arrIndexOut the index array of the resulting array \b arrOut
- * \sa MEDCouplingUMesh::ExtractFromIndexedArrays2
- */
-void MEDCouplingUMesh::ExtractFromIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
- DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut)
-{
- if(!arrIn || !arrIndxIn)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays : input pointer is NULL !");
- arrIn->checkAllocated(); arrIndxIn->checkAllocated();
- if(arrIn->getNumberOfComponents()!=1 || arrIndxIn->getNumberOfComponents()!=1)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays : input arrays must have exactly one component !");
- std::size_t sz=std::distance(idsOfSelectBg,idsOfSelectEnd);
- const int *arrInPtr=arrIn->getConstPointer();
- const int *arrIndxPtr=arrIndxIn->getConstPointer();
- int nbOfGrps=arrIndxIn->getNumberOfTuples()-1;
- if(nbOfGrps<0)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays : The format of \"arrIndxIn\" is invalid ! Its nb of tuples should be >=1 !");
- int maxSizeOfArr=arrIn->getNumberOfTuples();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arro=DataArrayInt::New();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arrIo=DataArrayInt::New();
- arrIo->alloc((int)(sz+1),1);
- const int *idsIt=idsOfSelectBg;
- int *work=arrIo->getPointer();
- *work++=0;
- int lgth=0;
- for(std::size_t i=0;i<sz;i++,work++,idsIt++)
- {
- if(*idsIt>=0 && *idsIt<nbOfGrps)
- lgth+=arrIndxPtr[*idsIt+1]-arrIndxPtr[*idsIt];
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::ExtractFromIndexedArrays : id located on pos #" << i << " value is " << *idsIt << " ! Must be in [0," << nbOfGrps << ") !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- if(lgth>=work[-1])
- *work=lgth;
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::ExtractFromIndexedArrays : id located on pos #" << i << " value is " << *idsIt << " and at this pos arrIndxIn[" << *idsIt;
- oss << "+1]-arrIndxIn[" << *idsIt << "] < 0 ! The input index array is bugged !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
- arro->alloc(lgth,1);
- work=arro->getPointer();
- idsIt=idsOfSelectBg;
- for(std::size_t i=0;i<sz;i++,idsIt++)
- {
- if(arrIndxPtr[*idsIt]>=0 && arrIndxPtr[*idsIt+1]<=maxSizeOfArr)
- work=std::copy(arrInPtr+arrIndxPtr[*idsIt],arrInPtr+arrIndxPtr[*idsIt+1],work);
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::ExtractFromIndexedArrays : id located on pos #" << i << " value is " << *idsIt << " arrIndx[" << *idsIt << "] must be >= 0 and arrIndx[";
- oss << *idsIt << "+1] <= " << maxSizeOfArr << " (the size of arrIn)!";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
- arrOut=arro.retn();
- arrIndexOut=arrIo.retn();
-}
-
-/*!
- * This method works on a pair input (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn
- * (\ref numbering-indirect).
- * This method returns the result of the extraction ( specified by a set of ids with a slice given by \a idsOfSelectStart, \a idsOfSelectStop and \a idsOfSelectStep ).
- * The selection of extraction is done standardly in new2old format.
- * This method returns indexed arrays (\ref numbering-indirect) using 2 arrays (arrOut,arrIndexOut).
- *
- * \param [in] idsOfSelectStart begin of set of ids of the input extraction (included)
- * \param [in] idsOfSelectStop end of set of ids of the input extraction (excluded)
- * \param [in] idsOfSelectStep
- * \param [in] arrIn arr origin array from which the extraction will be done.
- * \param [in] arrIndxIn is the input index array allowing to walk into \b arrIn
- * \param [out] arrOut the resulting array
- * \param [out] arrIndexOut the index array of the resulting array \b arrOut
- * \sa MEDCouplingUMesh::ExtractFromIndexedArrays
- */
-void MEDCouplingUMesh::ExtractFromIndexedArrays2(int idsOfSelectStart, int idsOfSelectStop, int idsOfSelectStep, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
- DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut)
-{
- if(!arrIn || !arrIndxIn)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays2 : input pointer is NULL !");
- arrIn->checkAllocated(); arrIndxIn->checkAllocated();
- if(arrIn->getNumberOfComponents()!=1 || arrIndxIn->getNumberOfComponents()!=1)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays2 : input arrays must have exactly one component !");
- int sz=DataArrayInt::GetNumberOfItemGivenBESRelative(idsOfSelectStart,idsOfSelectStop,idsOfSelectStep,"MEDCouplingUMesh::ExtractFromIndexedArrays2 : Input slice ");
- const int *arrInPtr=arrIn->getConstPointer();
- const int *arrIndxPtr=arrIndxIn->getConstPointer();
- int nbOfGrps=arrIndxIn->getNumberOfTuples()-1;
- if(nbOfGrps<0)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays2 : The format of \"arrIndxIn\" is invalid ! Its nb of tuples should be >=1 !");
- int maxSizeOfArr=arrIn->getNumberOfTuples();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arro=DataArrayInt::New();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arrIo=DataArrayInt::New();
- arrIo->alloc((int)(sz+1),1);
- int idsIt=idsOfSelectStart;
- int *work=arrIo->getPointer();
- *work++=0;
- int lgth=0;
- for(int i=0;i<sz;i++,work++,idsIt+=idsOfSelectStep)
- {
- if(idsIt>=0 && idsIt<nbOfGrps)
- lgth+=arrIndxPtr[idsIt+1]-arrIndxPtr[idsIt];
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::ExtractFromIndexedArrays2 : id located on pos #" << i << " value is " << idsIt << " ! Must be in [0," << nbOfGrps << ") !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- if(lgth>=work[-1])
- *work=lgth;
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::ExtractFromIndexedArrays2 : id located on pos #" << i << " value is " << idsIt << " and at this pos arrIndxIn[" << idsIt;
- oss << "+1]-arrIndxIn[" << idsIt << "] < 0 ! The input index array is bugged !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
- arro->alloc(lgth,1);
- work=arro->getPointer();
- idsIt=idsOfSelectStart;
- for(int i=0;i<sz;i++,idsIt+=idsOfSelectStep)
- {
- if(arrIndxPtr[idsIt]>=0 && arrIndxPtr[idsIt+1]<=maxSizeOfArr)
- work=std::copy(arrInPtr+arrIndxPtr[idsIt],arrInPtr+arrIndxPtr[idsIt+1],work);
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::ExtractFromIndexedArrays2 : id located on pos #" << i << " value is " << idsIt << " arrIndx[" << idsIt << "] must be >= 0 and arrIndx[";
- oss << idsIt << "+1] <= " << maxSizeOfArr << " (the size of arrIn)!";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
- arrOut=arro.retn();
- arrIndexOut=arrIo.retn();
-}
-
-/*!
- * This method works on an input pair (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn.
- * This method builds an output pair (\b arrOut,\b arrIndexOut) that is a copy from \b arrIn for all cell ids \b not \b in [ \b idsOfSelectBg , \b idsOfSelectEnd ) and for
- * cellIds \b in [ \b idsOfSelectBg , \b idsOfSelectEnd ) a copy coming from the corresponding values in input pair (\b srcArr, \b srcArrIndex).
- * This method is an generalization of MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx that performs the same thing but by without building explicitely a result output arrays.
- *
- * \param [in] idsOfSelectBg begin of set of ids of the input extraction (included)
- * \param [in] idsOfSelectEnd end of set of ids of the input extraction (excluded)
- * \param [in] arrIn arr origin array from which the extraction will be done.
- * \param [in] arrIndxIn is the input index array allowing to walk into \b arrIn
- * \param [in] srcArr input array that will be used as source of copy for ids in [ \b idsOfSelectBg, \b idsOfSelectEnd )
- * \param [in] srcArrIndex index array of \b srcArr
- * \param [out] arrOut the resulting array
- * \param [out] arrIndexOut the index array of the resulting array \b arrOut
- *
- * \sa MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx
- */
-void MEDCouplingUMesh::SetPartOfIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
- const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex,
- DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut)
-{
- if(arrIn==0 || arrIndxIn==0 || srcArr==0 || srcArrIndex==0)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::SetPartOfIndexedArrays : presence of null pointer in input parameter !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arro=DataArrayInt::New();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> arrIo=DataArrayInt::New();
- int nbOfTuples=arrIndxIn->getNumberOfTuples()-1;
- std::vector<bool> v(nbOfTuples,true);
- int offset=0;
- const int *arrIndxInPtr=arrIndxIn->getConstPointer();
- const int *srcArrIndexPtr=srcArrIndex->getConstPointer();
- for(const int *it=idsOfSelectBg;it!=idsOfSelectEnd;it++,srcArrIndexPtr++)
- {
- if(*it>=0 && *it<nbOfTuples)
- {
- v[*it]=false;
- offset+=(srcArrIndexPtr[1]-srcArrIndexPtr[0])-(arrIndxInPtr[*it+1]-arrIndxInPtr[*it]);
- }
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::SetPartOfIndexedArrays : On pos #" << std::distance(idsOfSelectBg,it) << " value is " << *it << " not in [0," << nbOfTuples << ") !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
- srcArrIndexPtr=srcArrIndex->getConstPointer();
- arrIo->alloc(nbOfTuples+1,1);
- arro->alloc(arrIn->getNumberOfTuples()+offset,1);
- const int *arrInPtr=arrIn->getConstPointer();
- const int *srcArrPtr=srcArr->getConstPointer();
- int *arrIoPtr=arrIo->getPointer(); *arrIoPtr++=0;
- int *arroPtr=arro->getPointer();
- for(int ii=0;ii<nbOfTuples;ii++,arrIoPtr++)
- {
- if(v[ii])
- {
- arroPtr=std::copy(arrInPtr+arrIndxInPtr[ii],arrInPtr+arrIndxInPtr[ii+1],arroPtr);
- *arrIoPtr=arrIoPtr[-1]+(arrIndxInPtr[ii+1]-arrIndxInPtr[ii]);
- }
- else
- {
- std::size_t pos=std::distance(idsOfSelectBg,std::find(idsOfSelectBg,idsOfSelectEnd,ii));
- arroPtr=std::copy(srcArrPtr+srcArrIndexPtr[pos],srcArrPtr+srcArrIndexPtr[pos+1],arroPtr);
- *arrIoPtr=arrIoPtr[-1]+(srcArrIndexPtr[pos+1]-srcArrIndexPtr[pos]);
- }
- }
- arrOut=arro.retn();
- arrIndexOut=arrIo.retn();
-}
-
-/*!
- * This method works on an input pair (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn.
- * This method is an specialization of MEDCouplingUMesh::SetPartOfIndexedArrays in the case of assignement do not modify the index in \b arrIndxIn.
- *
- * \param [in] idsOfSelectBg begin of set of ids of the input extraction (included)
- * \param [in] idsOfSelectEnd end of set of ids of the input extraction (excluded)
- * \param [in,out] arrInOut arr origin array from which the extraction will be done.
- * \param [in] arrIndxIn is the input index array allowing to walk into \b arrIn
- * \param [in] srcArr input array that will be used as source of copy for ids in [ \b idsOfSelectBg , \b idsOfSelectEnd )
- * \param [in] srcArrIndex index array of \b srcArr
- *
- * \sa MEDCouplingUMesh::SetPartOfIndexedArrays
- */
-void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx(const int *idsOfSelectBg, const int *idsOfSelectEnd, DataArrayInt *arrInOut, const DataArrayInt *arrIndxIn,
- const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex)
-{
- if(arrInOut==0 || arrIndxIn==0 || srcArr==0 || srcArrIndex==0)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx : presence of null pointer in input parameter !");
- int nbOfTuples=arrIndxIn->getNumberOfTuples()-1;
- const int *arrIndxInPtr=arrIndxIn->getConstPointer();
- const int *srcArrIndexPtr=srcArrIndex->getConstPointer();
- int *arrInOutPtr=arrInOut->getPointer();
- const int *srcArrPtr=srcArr->getConstPointer();
- for(const int *it=idsOfSelectBg;it!=idsOfSelectEnd;it++,srcArrIndexPtr++)
- {
- if(*it>=0 && *it<nbOfTuples)
- {
- if(srcArrIndexPtr[1]-srcArrIndexPtr[0]==arrIndxInPtr[*it+1]-arrIndxInPtr[*it])
- std::copy(srcArrPtr+srcArrIndexPtr[0],srcArrPtr+srcArrIndexPtr[1],arrInOutPtr+arrIndxInPtr[*it]);
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx : On pos #" << std::distance(idsOfSelectBg,it) << " id (idsOfSelectBg[" << std::distance(idsOfSelectBg,it)<< "]) is " << *it << " arrIndxIn[id+1]-arrIndxIn[id]!=srcArrIndex[pos+1]-srcArrIndex[pos] !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
- else
- {
- std::ostringstream oss; oss << "MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx : On pos #" << std::distance(idsOfSelectBg,it) << " value is " << *it << " not in [0," << nbOfTuples << ") !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
- }
-}
-
-/*!
- * This method works on a pair input (\b arrIn, \b arrIndxIn) where \b arr indexes is in \b arrIndxIn.
- * This method expects that these two input arrays come from the output of MEDCouplingUMesh::computeNeighborsOfCells method.
- * This method start from id 0 that will be contained in output DataArrayInt. It searches then all neighbors of id0 looking at arrIn[arrIndxIn[0]:arrIndxIn[0+1]].
- * Then it is repeated recursively until either all ids are fetched or no more ids are reachable step by step.
- * A negative value in \b arrIn means that it is ignored.
- * This method is useful to see if a mesh is contiguous regarding its connectivity. If it is not the case the size of returned array is different from arrIndxIn->getNumberOfTuples()-1.
- *