-void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<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< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
- MCAuto<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< MCAuto<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++;
- //
- MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
- int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
- //
- MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
- retTmp->setCoords(splitMesh1D->getCoords());
- retTmp->allocateCells();
-
- std::vector< std::vector<int> > out0;
- std::vector< std::vector< MCAuto<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,
- MCAuto<DataArrayInt>& idsLeftRight)
-{
- const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
- //
- std::vector<int> allEdges;
- std::vector< MCAuto<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< MCAuto<INTERP_KERNEL::Node>,int> m;
- MCAuto<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);
- MCAuto<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++)
- {
- MCAuto<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());
- MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
- MCAuto<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;
- }
- //
- MCAuto<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
- MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
- MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
- MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
- idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
- MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
- MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
- // deal with cells in mesh2D that are not cut but only some of their edges are
- MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
- idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
- idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
- MCAuto<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);
- MCAuto<DataArrayInt> outi0s(outi0);
- out0s=out0;
- out0s=out0s->buildUnique();
- out0s->sort(true);
- }
- //
- MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
- MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
- MCAuto<DataArrayInt> elts,eltsIndex;
- mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
- MCAuto<DataArrayInt> eltsIndex2(eltsIndex->deltaShiftIndex());
- MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
- if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
- throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
- MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
- MCAuto<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< MCAuto<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
- MCAuto<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());
- MCAuto<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->findIdFirstEqual(*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< MCAuto<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++)
- {
- MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
- idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
- MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
- MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
- MCAuto<DataArrayInt> partOfRet3;
- MCAuto<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);
- MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
- // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
- ret3->rearrange(1);
- MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStricltyNegative());
- for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
- {
- int old2DCellId(-ret3->getIJ(*it,0)-1);
- MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
- ret3->setIJ(*it,0,FindRightCandidateAmong(ret2D,candidates->begin(),candidates->end(),ret1,*it%2==0?-((*it)/2+1):(*it)/2+1,eps));// div by 2 because 2 components natively in ret3
- }
- ret3->changeValue(std::numeric_limits<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()->begin());
- const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
- int offset1(m1->getNumberOfNodes());
- const double *coo2(m2->getCoords()->begin());
- const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
- int offset2(offset1+m2->getNumberOfNodes());
- int offset3(offset2+((int)addCoords.size())/2);
- MCAuto<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());
- }
- }
- 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:
- MCAuto<DataArrayInt> _d(DataArrayInt::New()),_dI(DataArrayInt::New());
- MCAuto<DataArrayInt> _rD(DataArrayInt::New()),_rDI(DataArrayInt::New());
- MCAuto<MEDCouplingUMesh> m_points(buildDescendingConnectivity(_d, _dI, _rD, _rDI));
- const int *d(_d->begin()), *dI(_dI->begin());
- const int *rD(_rD->begin()), *rDI(_rDI->begin());
- MCAuto<DataArrayInt> _dsi(_rDI->deltaShiftIndex());
- const int * dsi(_dsi->begin());
- MCAuto<DataArrayInt> dsii = _dsi->findIdsNotInRange(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());
- MCAuto<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<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
-{
- MCAuto<INTERP_KERNEL::Node> nTmp(n); nTmp->incrRef();
- std::map<MCAuto<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<MCAuto<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
-