+//tony to put in private of MEDCouplingUMesh
+MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<int> >& intersectEdge2, const DataArrayDouble *coords1, const std::vector<double>& addCoo)
+{
+ int nCells(mesh1D->getNumberOfCells());
+ if(nCells!=(int)intersectEdge2.size())
+ throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
+ const DataArrayDouble *coo2(mesh1D->getCoords());
+ const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
+ const double *coo2Ptr(coo2->begin());
+ int offset1(coords1->getNumberOfTuples());
+ int offset2(offset1+coo2->getNumberOfTuples());
+ int offset3(offset2+addCoo.size()/2);
+ std::vector<double> addCooQuad;
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
+ int tmp[4],cicnt(0);
+ for(int i=0;i<nCells;i++)
+ {
+ std::map<INTERP_KERNEL::Node *,int> m;
+ INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
+ const std::vector<int>& subEdges(intersectEdge2[i]);
+ int nbSubEdge(subEdges.size()/2);
+ for(int j=0;j<nbSubEdge;j++)
+ {
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node> n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node> n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
+ INTERP_KERNEL::Edge *e2Ptr(e2);
+ if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
+ {
+ tmp[0]=INTERP_KERNEL::NORM_SEG3;
+ tmp[1]=subEdges[2*j]; tmp[2]=subEdges[2*j+1];
+ cicnt+=4;
+ cOut->insertAtTheEnd(tmp,tmp+4);
+ ciOut->pushBackSilent(cicnt);
+ }
+ else
+ {
+ tmp[0]=INTERP_KERNEL::NORM_SEG2;
+ tmp[1]=subEdges[2*j]; tmp[2]=subEdges[2*j+1]; tmp[3]=offset3+(int)addCooQuad.size()/2;
+ cicnt+=3;
+ cOut->insertAtTheEnd(tmp,tmp+3);
+ ciOut->pushBackSilent(cicnt);
+ }
+ }
+ //INTERP_KERNEL::Edge *e2(e->buildEdgeLyingOnMe());
+ for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it2=m.begin();it2!=m.end();it2++)
+ (*it2).first->decrRef();
+ e->decrRef();
+ }
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
+ ret->setConnectivity(cOut,ciOut,true);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr3(DataArrayDouble::New());
+ arr3->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,C_DEALLOC,(int)addCooQuad.size()/2,2);
+ std::vector<const DataArrayDouble *> coordss(4);
+ coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
+ ret->setCoords(arr);
+ return ret.retn();
+}
+
+/*!
+ * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
+ * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
+ * all nodes from \a mesh1D.
+ * 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.
+ * \param [in] mesh1D - the 1D mesh (spacedim=2 meshdim=1) the is the tool that will be used to intersect \a mesh2D.
+ * \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 that gives for each cell id \a i in \a splitMesh1D the 1 or 2 id(s) in \a splitMesh2D that \a i shares.
+ * 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.
+ */
+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();
+ 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());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
+ Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo);
+ MEDCouplingAutoRefCountObjectPtr<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();
+ //
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo));
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> baryRet1(ret1->getBarycenterAndOwner());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> elts,eltsIndex;
+ mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
+ splitMesh1D=ret1.retn();
+}