-/**
- * 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();
- const int *connI1=m1->getNodalConnectivityIndex()->getConstPointer();
- int offset1=m1->getNumberOfNodes();
- const double *coo2=m2->getCoords()->getConstPointer();
- const int *conn2=m2->getNodalConnectivity()->getConstPointer();
- const int *connI2=m2->getNodalConnectivityIndex()->getConstPointer();
- int offset2=offset1+m2->getNumberOfNodes();
- int offset3=offset2+((int)addCoords.size())/2;
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
- const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
- // Here a BBTree on 2D-cells, not on segments:
- BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
- int ncell1=m1->getNumberOfCells();
- crI.push_back(0);
- for(int i=0;i<ncell1;i++)
- {
- std::vector<int> candidates2;
- myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
- std::map<INTERP_KERNEL::Node *,int> mapp;
- std::map<int,INTERP_KERNEL::Node *> mappRev;
- INTERP_KERNEL::QuadraticPolygon pol1;
- INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
- const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
- // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
- MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
- // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
- pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
- desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
- //
- std::set<INTERP_KERNEL::Edge *> edges1;// store all edges of pol1 that are NOT consumed by intersect cells. If any after iteration over candidates2 -> a part of pol1 should appear in result
- std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
- INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
- for(it1.first();!it1.finished();it1.next())
- edges1.insert(it1.current()->getPtr());
- //
- std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
- std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
- int ii=0;
- for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
- {
- INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
- const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
- // Complete mapping with elements coming from the current cell it2 in mesh2:
- MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
- // pol2 is the new QP in the final merged result.
- pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
- pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
- }
- ii=0;
- for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
- {
- pol1.initLocationsWithOther(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();
- }
-}
-
-/*!
- * 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 !");
- }
- }
- }
-}