-/*!
- * \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<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);
- for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it2=m.begin();it2!=m.end();it2++)
- (*it2).first->decrRef();
- }
- }
- // 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<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();
- for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it2=m.begin();it2!=m.end();it2++)
- (*it2).first->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();
-}
-
-/*!
- * 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)
- throw(INTERP_KERNEL::Exception)
-{
- static const int SPACEDIM=2;
- // 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);
- const int *c1=m1Desc->getNodalConnectivity()->getConstPointer();
- const int *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();
- int 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=m1->getNumberOfNodes();
- int offset2=offset1+m2->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);
- delete pol2;
- delete pol1;
- }
- else
- intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i+1]);
- }
- 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) throw(INTERP_KERNEL::Exception)
-{
- 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 throw(INTERP_KERNEL::Exception)
-{
- 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);
- }
- }
-}
-