-class VectorOfCellInfo
-{
-public:
- VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<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 MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs);
- const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
- const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
- MCAuto<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< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
- const CellInfo& get(int pos) const;
- CellInfo& get(int pos);
-private:
- std::vector<CellInfo> _pool;
- MCAuto<MEDCouplingUMesh> _ze_mesh;
- std::vector<EdgeInfo> _edge_info;
-};
-
-VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<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 !");
- MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
- return zeMesh->getCellContainingPoint(barys->begin(),eps);
-}
-
-void VectorOfCellInfo::setMeshAt(int pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges, const std::vector< std::vector< MCAuto<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< MCAuto<MEDCouplingUMesh> > ms;
- if(pos>0)
- {
- MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
- ms.push_back(elt);
- }
- ms.push_back(mesh);
- if(pos<_ze_mesh->getNumberOfCells()-1)
- {
- MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(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< 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()->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);
- 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().c_str());
- }
- }
- for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
- (*it).second->decrRef();
- }
-}