-/*!
- * \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 have same nature each other (Linear,Quadratic).
- * 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
- *
- * 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 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);
- }
- }
-}
-