src/MEDCoupling/MEDCouplingUMesh_intersection.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingUMesh.hxx"
  22 #include "MEDCoupling1GTUMesh.hxx"
  23 #include "CellModel.hxx"
  24 #include "VolSurfUser.txx"
  25 #include "InterpolationUtils.hxx"
  26 #include "PointLocatorAlgos.txx"
  27 #include "BBTree.txx"
  28 #include "BBTreeDst.txx"
  29 #include "DirectedBoundingBox.hxx"
  30 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  31 #include "InterpKernelAutoPtr.hxx"
  32 #include "InterpKernelGeo2DNode.hxx"
  33 #include "InterpKernelGeo2DEdgeLin.hxx"
  34 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  35 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
  36
  37 #include <sstream>
  38 #include <fstream>
  39 #include <numeric>
  40 #include <cstring>
  41 #include <limits>
  42 #include <list>
  43
  44 using namespace MEDCoupling;
  45
  46 /// @cond INTERNAL
  47
  48 int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
  49 {
  50   if(id!=-1)
  51     return id;
  52   else
  53     {
  54       int ret(nodesCnter++);
  55       double newPt[2];
  56       e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt);
  57       addCoo.insertAtTheEnd(newPt,newPt+2);
  58       return ret;
  59     }
  60 }
  61
  62 int InternalAddPointOriented(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
  63 {
  64   if(id!=-1)
  65     return id;
  66   else
  67     {
  68       int ret(nodesCnter++);
  69       double newPt[2];
  70       e->getMiddleOfPointsOriented(coo+2*startId,coo+2*endId,newPt);
  71       addCoo.insertAtTheEnd(newPt,newPt+2);
  72       return ret;
  73     }
  74 }
  75
  76
  77 void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
  78 {
  79   int tmp[3];
  80   int trueStart(start>=0?start:nbOfEdges+start);
  81   tmp[0]=linOrArc?(int)INTERP_KERNEL::NORM_QPOLYG:(int)INTERP_KERNEL::NORM_POLYGON; tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp];
  82   newConnOfCell->insertAtTheEnd(tmp,tmp+3);
  83   if(linOrArc)
  84     {
  85       if(stp-start>1)
  86         {
  87           int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
  88           InternalAddPointOriented(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2);
  89           middles.push_back(tmp3+offset);
  90         }
  91       else
  92         middles.push_back(connBg[trueStart+nbOfEdges]);
  93     }
  94 }
  95
  96 void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
  97 {
  98   int tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]);
  99   newConnOfCell->pushBackSilent(tmpEnd);
 100   if(linOrArc)
 101     {
 102       if(stp-start>1)
 103         {
 104           int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
 105           InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
 106           middles.push_back(tmp3+offset);
 107         }
 108       else
 109         middles.push_back(connBg[start+nbOfEdges]);
 110     }
 111 }
 112
 113 void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
 114 {
 115   // only the quadratic point to deal with:
 116   if(linOrArc)
 117     {
 118       if(stp-start>1)
 119         {
 120           int tmpSrt(connBg[start]),tmpEnd(connBg[stp]);
 121           int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
 122           InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
 123           middles.push_back(tmp3+offset);
 124         }
 125       else
 126         middles.push_back(connBg[start+nbOfEdges]);
 127     }
 128 }
 129
 130 void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
 131 {
 132   MCAuto<INTERP_KERNEL::Node> nTmp(n); nTmp->incrRef();
 133   std::map<MCAuto<INTERP_KERNEL::Node>,int>::const_iterator it(m.find(nTmp));
 134   if(it==m.end())
 135     throw INTERP_KERNEL::Exception("Internal error in remapping !");
 136   int v((*it).second);
 137   if(v==forbVal0 || v==forbVal1)
 138     return ;
 139   if(std::find(isect.begin(),isect.end(),v)==isect.end())
 140     isect.push_back(v);
 141 }
 142
 143 bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
 144 {
 145   int sz(c.size());
 146   if(sz<=1)
 147     return false;
 148   bool presenceOfOn(false);
 149   for(int i=0;i<sz;i++)
 150     {
 151       INTERP_KERNEL::ElementaryEdge *e(c[i]);
 152       if(e->getLoc()!=INTERP_KERNEL::FULL_ON_1)
 153         continue ;
 154       IKGeo2DInternalMapper2(e->getStartNode(),m,forbVal0,forbVal1,isect);
 155       IKGeo2DInternalMapper2(e->getEndNode(),m,forbVal0,forbVal1,isect);
 156     }
 157   return presenceOfOn;
 158 }
 159
 160
 161 namespace MEDCoupling
 162 {
 163
 164   INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge2(INTERP_KERNEL::NormalizedCellType typ, const int *bg, const double *coords2D, std::map< MCAuto<INTERP_KERNEL::Node>,int>& m)
 165   {
 166     INTERP_KERNEL::Edge *ret(0);
 167     MCAuto<INTERP_KERNEL::Node> n0(new INTERP_KERNEL::Node(coords2D[2*bg[0]],coords2D[2*bg[0]+1])),n1(new INTERP_KERNEL::Node(coords2D[2*bg[1]],coords2D[2*bg[1]+1]));
 168     m[n0]=bg[0]; m[n1]=bg[1];
 169     switch(typ)
 170     {
 171       case INTERP_KERNEL::NORM_SEG2:
 172         {
 173           ret=new INTERP_KERNEL::EdgeLin(n0,n1);
 174           break;
 175         }
 176       case INTERP_KERNEL::NORM_SEG3:
 177         {
 178           INTERP_KERNEL::Node *n2(new INTERP_KERNEL::Node(coords2D[2*bg[2]],coords2D[2*bg[2]+1])); m[n2]=bg[2];
 179           INTERP_KERNEL::EdgeLin *e1(new INTERP_KERNEL::EdgeLin(n0,n2)),*e2(new INTERP_KERNEL::EdgeLin(n2,n1));
 180           INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
 181           // is the SEG3 degenerated, and thus can be reduced to a SEG2?
 182           bool colinearity(inters.areColinears());
 183           delete e1; delete e2;
 184           if(colinearity)
 185             { ret=new INTERP_KERNEL::EdgeLin(n0,n1); }
 186           else
 187             { ret=new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1); }
 188           break;
 189         }
 190       default:
 191         throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge2 : Expecting a mesh with spaceDim==2 and meshDim==1 !");
 192     }
 193     return ret;
 194   }
 195
 196   INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >& mapp2, const int *bg)
 197   {
 198     INTERP_KERNEL::Edge *ret=0;
 199     switch(typ)
 200     {
 201       case INTERP_KERNEL::NORM_SEG2:
 202         {
 203           ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
 204           break;
 205         }
 206       case INTERP_KERNEL::NORM_SEG3:
 207         {
 208           INTERP_KERNEL::EdgeLin *e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[2]].first);
 209           INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]].first,mapp2[bg[1]].first);
 210           INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
 211           // is the SEG3 degenerated, and thus can be reduced to a SEG2?
 212           bool colinearity=inters.areColinears();
 213           delete e1; delete e2;
 214           if(colinearity)
 215             ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
 216           else
 217             ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]].first,mapp2[bg[2]].first,mapp2[bg[1]].first);
 218           mapp2[bg[2]].second=false;
 219           break;
 220         }
 221       default:
 222         throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !");
 223     }
 224     return ret;
 225   }
 226
 227   /*!
 228    * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed by the sub set of cells in 'candidates' taken from
 229    * the global mesh 'mDesc'.
 230    * The input mesh 'mDesc' must be so that mDim==1 and spaceDim==2.
 231    * 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'.
 232    */
 233   INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates,
 234                                                                    std::map<INTERP_KERNEL::Node *,int>& mapp)
 235   {
 236     mapp.clear();
 237     std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;//bool is for a flag specifying if node is boundary (true) or only a middle for SEG3.
 238     const double *coo=mDesc->getCoords()->getConstPointer();
 239     const int *c=mDesc->getNodalConnectivity()->getConstPointer();
 240     const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
 241     std::set<int> s;
 242     for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
 243       s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
 244     for(std::set<int>::const_iterator it2=s.begin();it2!=s.end();it2++)
 245       {
 246         INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]);
 247         mapp2[*it2]=std::pair<INTERP_KERNEL::Node *,bool>(n,true);
 248       }
 249     INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon;
 250     for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
 251       {
 252         INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]];
 253         ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1));
 254       }
 255     for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++)
 256       {
 257         if((*it2).second.second)
 258           mapp[(*it2).second.first]=(*it2).first;
 259         ((*it2).second.first)->decrRef();
 260       }
 261     return ret;
 262   }
 263
 264   INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo)
 265   {
 266     if(nodeId>=offset2)
 267       {
 268         int locId=nodeId-offset2;
 269         return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
 270       }
 271     if(nodeId>=offset1)
 272       {
 273         int locId=nodeId-offset1;
 274         return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]);
 275       }
 276     return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]);
 277   }
 278
 279   /**
 280    * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI).
 281    */
 282   void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
 283                                         const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
 284                                         /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
 285   {
 286     for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
 287       {
 288         int eltId1=abs(*desc1)-1;
 289         for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
 290           {
 291             std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.find(*it1);
 292             if(it==mappRev.end())
 293               {
 294                 INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo);
 295                 mapp[node]=*it1;
 296                 mappRev[*it1]=node;
 297               }
 298           }
 299       }
 300   }
 301 }
 302
 303
 304
 305 /*!
 306  * Returns true if a colinearization has been found in the given cell. If false is returned the content pushed in \a newConnOfCell is equal to [ \a connBg , \a connEnd ) .
 307  * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
 308  */
 309 bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
 310 {
 311   std::size_t sz(std::distance(connBg,connEnd));
 312   if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
 313     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
 314   sz--;
 315   INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
 316   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
 317   unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
 318   unsigned nbOfHit(0); // number of fusions operated
 319   int posBaseElt(0),posEndElt(0),nbOfTurn(0);
 320   const unsigned int maxNbOfHit = cm.isQuadratic() ? nbs-2 : nbs-3;  // a quad cell is authorized to end up with only two edges, a linear one has to keep 3 at least
 321   INTERP_KERNEL::NormalizedCellType typeOfSon;
 322   std::vector<int> middles;
 323   bool ret(false);
 324   for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
 325     {
 326       cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
 327       std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
 328       INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
 329       posEndElt = posBaseElt+1;
 330
 331       // Look backward first: are the final edges of the cells colinear with the first ones?
 332       // This initializes posBaseElt.
 333       if(nbOfTurn==0)
 334         {
 335           for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell wih one single edge
 336             {
 337               cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn,typeOfSon);
 338               INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
 339               INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
 340               bool isColinear=eint->areColinears();
 341               if(isColinear)
 342                 {
 343                   nbOfHit++;
 344                   posBaseElt--;
 345                   ret=true;
 346                 }
 347               delete eint;
 348               eCand->decrRef();
 349               if(!isColinear)
 350                 break;
 351             }
 352         }
 353       // Now move forward:
 354       const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt);  // the first element to be inspected going forward
 355       for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++)  // 2nd condition is to avoid ending with a cell wih one single edge
 356         {
 357           cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn,typeOfSon); // get edge #j's connectivity
 358           INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
 359           INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
 360           bool isColinear(eint->areColinears());
 361           if(isColinear)
 362             {
 363               nbOfHit++;
 364               posEndElt++;
 365               ret=true;
 366             }
 367           delete eint;
 368           eCand->decrRef();
 369           if(!isColinear)
 370               break;
 371         }
 372       //push [posBaseElt,posEndElt) in newConnOfCell using e
 373       // The if clauses below are (volontary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its begining!
 374       if(nbOfTurn==0)
 375         // at the begining of the connectivity (insert type)
 376         EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
 377       else if((nbOfHit+nbOfTurn) != (nbs-1))
 378         // in the middle
 379         EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
 380       if ((nbOfHit+nbOfTurn) == (nbs-1))
 381         // at the end (only quad points to deal with)
 382         EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
 383       posBaseElt=posEndElt;
 384       e->decrRef();
 385     }
 386   if(!middles.empty())
 387     newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
 388   return ret;
 389 }
 390
 391
 392
 393 bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1, const std::vector<int>& candidates, int start, int stop, int& retVal)
 394 {
 395   if(candidates.empty())
 396     return false;
 397   for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
 398     {
 399       const std::vector<int>& pool(intersectEdge1[*it]);
 400       int tmp[2]; tmp[0]=start; tmp[1]=stop;
 401       if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
 402         {
 403           retVal=*it+1;
 404           return true;
 405         }
 406       tmp[0]=stop; tmp[1]=start;
 407       if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
 408         {
 409           retVal=-*it-1;
 410           return true;
 411         }
 412     }
 413   return false;
 414 }
 415
 416 /*!
 417  * This method performs the 2nd step of Partition of 2D mesh.
 418  * This method has 4 inputs :
 419  *  - a mesh 'm1' with meshDim==1 and a SpaceDim==2
 420  *  - a mesh 'm2' with meshDim==1 and a SpaceDim==2
 421  *  - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids randomly sorted.
 422  * 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'
 423  * Nodes end up lying consecutively on a cutted edge.
 424  * \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.
 425  * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
 426  * \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.
 427  * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
 428  * \param[out] intersectEdge the same content as subDiv, but correclty oriented.
 429  */
 430 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
 431                                            const std::vector<double>& addCoo,
 432                                            const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
 433 {
 434   int offset1=m1->getNumberOfNodes();
 435   int ncell=m2->getNumberOfCells();
 436   const int *c=m2->getNodalConnectivity()->begin();
 437   const int *cI=m2->getNodalConnectivityIndex()->begin();
 438   const double *coo=m2->getCoords()->begin();
 439   const double *cooBis=m1->getCoords()->begin();
 440   int offset2=offset1+m2->getNumberOfNodes();
 441   intersectEdge.resize(ncell);
 442   for(int i=0;i<ncell;i++,cI++)
 443     {
 444       const std::vector<int>& divs=subDiv[i];
 445       int nnode=cI[1]-cI[0]-1;
 446       std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;
 447       std::map<INTERP_KERNEL::Node *, int> mapp22;
 448       for(int j=0;j<nnode;j++)
 449         {
 450           INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
 451           int nnid=c[(*cI)+j+1];
 452           mapp2[nnid]=std::pair<INTERP_KERNEL::Node *,bool>(nn,true);
 453           mapp22[nn]=nnid+offset1;
 454         }
 455       INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
 456       for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
 457         ((*it).second.first)->decrRef();
 458       std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
 459       std::map<INTERP_KERNEL::Node *,int> mapp3;
 460       for(std::size_t j=0;j<divs.size();j++)
 461         {
 462           int id=divs[j];
 463           INTERP_KERNEL::Node *tmp=0;
 464           if(id<offset1)
 465             tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
 466           else if(id<offset2)
 467             tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
 468           else
 469             tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
 470           addNodes[j]=tmp;
 471           mapp3[tmp]=id;
 472         }
 473       e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
 474       for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
 475         (*it)->decrRef();
 476       e->decrRef();
 477     }
 478 }
 479
 480 MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<int> >& intersectEdge2, const DataArrayDouble *coords1, const std::vector<double>& addCoo, const std::map<int,int>& mergedNodes, const std::vector< std::vector<int> >& colinear2, const std::vector< std::vector<int> >& intersectEdge1,
 481                                      MCAuto<DataArrayInt>& idsInRetColinear, MCAuto<DataArrayInt>& idsInMesh1DForIdsInRetColinear)
 482 {
 483   idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1);
 484   idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
 485   int nCells(mesh1D->getNumberOfCells());
 486   if(nCells!=(int)intersectEdge2.size())
 487     throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
 488   const DataArrayDouble *coo2(mesh1D->getCoords());
 489   const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
 490   const double *coo2Ptr(coo2->begin());
 491   int offset1(coords1->getNumberOfTuples());
 492   int offset2(offset1+coo2->getNumberOfTuples());
 493   int offset3(offset2+addCoo.size()/2);
 494   std::vector<double> addCooQuad;
 495   MCAuto<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
 496   int tmp[4],cicnt(0),kk(0);
 497   for(int i=0;i<nCells;i++)
 498     {
 499       std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
 500       INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
 501       const std::vector<int>& subEdges(intersectEdge2[i]);
 502       int nbSubEdge(subEdges.size()/2);
 503       for(int j=0;j<nbSubEdge;j++,kk++)
 504         {
 505           MCAuto<INTERP_KERNEL::Node> n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo)),n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
 506           MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
 507           INTERP_KERNEL::Edge *e2Ptr(e2);
 508           std::map<int,int>::const_iterator itm;
 509           if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
 510             {
 511               tmp[0]=INTERP_KERNEL::NORM_SEG3;
 512               itm=mergedNodes.find(subEdges[2*j]);
 513               tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
 514               itm=mergedNodes.find(subEdges[2*j+1]);
 515               tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
 516               tmp[3]=offset3+(int)addCooQuad.size()/2;
 517               double tmp2[2];
 518               e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
 519               cicnt+=4;
 520               cOut->insertAtTheEnd(tmp,tmp+4);
 521               ciOut->pushBackSilent(cicnt);
 522             }
 523           else
 524             {
 525               tmp[0]=INTERP_KERNEL::NORM_SEG2;
 526               itm=mergedNodes.find(subEdges[2*j]);
 527               tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
 528               itm=mergedNodes.find(subEdges[2*j+1]);
 529               tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
 530               cicnt+=3;
 531               cOut->insertAtTheEnd(tmp,tmp+3);
 532               ciOut->pushBackSilent(cicnt);
 533             }
 534           int tmp00;
 535           if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
 536             {
 537               idsInRetColinear->pushBackSilent(kk);
 538               idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
 539             }
 540         }
 541       e->decrRef();
 542     }
 543   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
 544   ret->setConnectivity(cOut,ciOut,true);
 545   MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
 546   arr3->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
 547   MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,C_DEALLOC,(int)addCooQuad.size()/2,2);
 548   std::vector<const DataArrayDouble *> coordss(4);
 549   coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
 550   MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
 551   ret->setCoords(arr);
 552   return ret.retn();
 553 }
 554
 555 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
 556 {
 557   std::vector<int> allEdges;
 558   for(const int *it2(descBg);it2!=descEnd;it2++)
 559     {
 560       const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
 561       if(*it2>0)
 562         allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
 563       else
 564         allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
 565     }
 566   std::size_t nb(allEdges.size());
 567   if(nb%2!=0)
 568     throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
 569   std::size_t nbOfEdgesOf2DCellSplit(nb/2);
 570   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
 571   ret->setCoords(coords);
 572   ret->allocateCells(1);
 573   std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
 574   for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
 575     connOut[kk]=allEdges[2*kk];
 576   ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]);
 577   return ret.retn();
 578 }
 579
 580 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
 581 {
 582   const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
 583   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
 584   std::size_t ii(0);
 585   unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
 586   if(sz!=std::distance(descBg,descEnd))
 587     throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
 588   INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
 589   std::vector<int> allEdges,centers;
 590   const double *coordsPtr(coords->begin());
 591   MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
 592   int offset(coords->getNumberOfTuples());
 593   for(const int *it2(descBg);it2!=descEnd;it2++,ii++)
 594     {
 595       INTERP_KERNEL::NormalizedCellType typeOfSon;
 596       cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
 597       const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
 598       if(*it2>0)
 599         allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
 600       else
 601         allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
 602       if(edge1.size()==2)
 603         centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
 604       else
 605         {//the current edge has been subsplit -> create corresponding centers.
 606           std::size_t nbOfCentersToAppend(edge1.size()/2);
 607           std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
 608           MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
 609           std::vector<int>::const_iterator it3(allEdges.end()-edge1.size());
 610           for(std::size_t k=0;k<nbOfCentersToAppend;k++)
 611             {
 612               double tmpp[2];
 613               const double *aa(coordsPtr+2*(*it3++));
 614               const double *bb(coordsPtr+2*(*it3++));
 615               ee->getMiddleOfPoints(aa,bb,tmpp);
 616               addCoo->insertAtTheEnd(tmpp,tmpp+2);
 617               centers.push_back(offset+k);
 618             }
 619         }
 620     }
 621   std::size_t nb(allEdges.size());
 622   if(nb%2!=0)
 623     throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
 624   std::size_t nbOfEdgesOf2DCellSplit(nb/2);
 625   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
 626   if(addCoo->empty())
 627     ret->setCoords(coords);
 628   else
 629     {
 630       addCoo->rearrange(2);
 631       addCoo=DataArrayDouble::Aggregate(coords,addCoo);
 632       ret->setCoords(addCoo);
 633     }
 634   ret->allocateCells(1);
 635   std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
 636   for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
 637     connOut[kk]=allEdges[2*kk];
 638   connOut.insert(connOut.end(),centers.begin(),centers.end());
 639   ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]);
 640   return ret.retn();
 641 }
 642
 643 /*!
 644  * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
 645  * of those edges.
 646  *
 647  * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
 648  */
 649 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
 650 {
 651   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
 652   if(!cm.isQuadratic())
 653     return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
 654   else
 655     return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
 656 }
 657
 658 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<int>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
 659 {
 660   bool isQuad(false);
 661   for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
 662     {
 663       const INTERP_KERNEL::Edge *ee(*it);
 664       if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
 665         isQuad=true;
 666     }
 667   if(!isQuad)
 668     mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]);
 669   else
 670     {
 671       const double *coo(mesh2D->getCoords()->begin());
 672       std::size_t sz(conn.size());
 673       std::vector<double> addCoo;
 674       std::vector<int> conn2(conn);
 675       int offset(mesh2D->getNumberOfNodes());
 676       for(std::size_t i=0;i<sz;i++)
 677         {
 678           double tmp[2];
 679           edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
 680           addCoo.insert(addCoo.end(),tmp,tmp+2);
 681           conn2.push_back(offset+(int)i);
 682         }
 683       mesh2D->getCoords()->rearrange(1);
 684       mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
 685       mesh2D->getCoords()->rearrange(2);
 686       mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]);
 687     }
 688 }
 689
 690 /*!
 691  * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
 692  *
 693  * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
 694  * a set of edges defined in \a splitMesh1D.
 695  */
 696 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
 697                              std::vector< std::vector<int> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
 698 {
 699   std::size_t nb(edge1Bis.size()/2);
 700   std::size_t nbOfEdgesOf2DCellSplit(nb/2);
 701   int iEnd(splitMesh1D->getNumberOfCells());
 702   if(iEnd==0)
 703     throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
 704   std::size_t ii,jj;
 705   const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
 706   for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
 707   for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
 708   //
 709   if(jj==nb)
 710     {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
 711       out0.resize(1); out1.resize(1);
 712       std::vector<int>& connOut(out0[0]);
 713       connOut.resize(nbOfEdgesOf2DCellSplit);
 714       std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
 715       edgesPtr.resize(nbOfEdgesOf2DCellSplit);
 716       for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
 717         {
 718           connOut[kk]=edge1Bis[2*kk];
 719           edgesPtr[kk]=edge1BisPtr[2*kk];
 720         }
 721     }
 722   else
 723     {
 724       // [i,iEnd[ contains the
 725       out0.resize(2); out1.resize(2);
 726       std::vector<int>& connOutLeft(out0[0]);
 727       std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
 728       std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
 729       std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
 730       for(std::size_t k=ii;k<jj+1;k++)
 731         { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
 732       std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
 733       for(int ik=0;ik<iEnd;ik++)
 734         {
 735           std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
 736           MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
 737           ees[ik]=ee;
 738         }
 739       for(int ik=iEnd-1;ik>=0;ik--)
 740         connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
 741       for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
 742         { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
 743       eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
 744       for(int ik=0;ik<iEnd;ik++)
 745         connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
 746       eright.insert(eright.end(),ees.begin(),ees.end());
 747     }
 748 }
 749
 750 struct CellInfo
 751 {
 752 public:
 753   CellInfo() { }
 754   CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
 755 public:
 756   std::vector<int> _edges;
 757   std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
 758 };
 759
 760 CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
 761 {
 762   std::size_t nbe(edges.size());
 763   std::vector<int> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
 764   for(std::size_t i=0;i<nbe;i++)
 765     {
 766       edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
 767       edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
 768     }
 769   _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
 770   std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
 771   std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
 772 }
 773
 774 class EdgeInfo
 775 {
 776 public:
 777   EdgeInfo(int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
 778   EdgeInfo(int istart, int iend, int pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
 779   bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; }
 780   void somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
 781   void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
 782 private:
 783   int _istart;
 784   int _iend;
 785   MCAuto<MEDCouplingUMesh> _mesh;
 786   MCAuto<INTERP_KERNEL::Edge> _edge;
 787   int _left;
 788   int _right;
 789 };
 790
 791 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
 792 {
 793   const MEDCouplingUMesh *mesh(_mesh);
 794   if(mesh)
 795     return ;
 796   if(_right<pos)
 797     return ;
 798   if(_left>pos)
 799     { _left++; _right++; return ; }
 800   if(_right==pos)
 801     {
 802       bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
 803       if((isLeft && isRight) || (!isLeft && !isRight))
 804         throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
 805       if(isLeft)
 806         return ;
 807       if(isRight)
 808         {
 809           _right++;
 810           return ;
 811         }
 812     }
 813   if(_left==pos)
 814     {
 815       bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
 816       if((isLeft && isRight) || (!isLeft && !isRight))
 817         throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
 818       if(isLeft)
 819         {
 820           _right++;
 821           return ;
 822         }
 823       if(isRight)
 824         {
 825           _left++;
 826           _right++;
 827           return ;
 828         }
 829     }
 830 }
 831
 832 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
 833 {
 834   const MEDCouplingUMesh *mesh(_mesh);
 835   if(!mesh)
 836     {
 837       neighbors[0]=offset+_left; neighbors[1]=offset+_right;
 838     }
 839   else
 840     {// not fully splitting cell case
 841       if(mesh2D->getNumberOfCells()==1)
 842         {//little optimization. 1 cell no need to find in which cell mesh is !
 843           neighbors[0]=offset; neighbors[1]=offset;
 844           return;
 845         }
 846       else
 847         {
 848           MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
 849           int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
 850           if(cellId==-1)
 851             throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
 852           neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
 853         }
 854     }
 855 }
 856
 857 class VectorOfCellInfo
 858 {
 859 public:
 860   VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
 861   std::size_t size() const { return _pool.size(); }
 862   int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
 863   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);
 864   const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
 865   const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
 866   MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
 867   void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
 868 private:
 869   int getZePosOfEdgeGivenItsGlobalId(int pos) const;
 870   void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
 871   const CellInfo& get(int pos) const;
 872   CellInfo& get(int pos);
 873 private:
 874   std::vector<CellInfo> _pool;
 875   MCAuto<MEDCouplingUMesh> _ze_mesh;
 876   std::vector<EdgeInfo> _edge_info;
 877 };
 878
 879 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
 880 {
 881   _pool[0]._edges=edges;
 882   _pool[0]._edges_ptr=edgesPtr;
 883 }
 884
 885 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
 886 {
 887   if(_pool.empty())
 888     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
 889   if(_pool.size()==1)
 890     return 0;
 891   const MEDCouplingUMesh *zeMesh(_ze_mesh);
 892   if(!zeMesh)
 893     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
 894   MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
 895   return zeMesh->getCellContainingPoint(barys->begin(),eps);
 896 }
 897
 898 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)
 899 {
 900   get(pos);//to check pos
 901   bool isFast(pos==0 && _pool.size()==1);
 902   std::size_t sz(edges.size());
 903   // dealing with edges
 904   if(sz==1)
 905     _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
 906   else
 907     _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
 908   //
 909   std::vector<CellInfo> pool(_pool.size()-1+sz);
 910   for(int i=0;i<pos;i++)
 911     pool[i]=_pool[i];
 912   for(std::size_t j=0;j<sz;j++)
 913     pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
 914   for(int i=pos+1;i<(int)_pool.size();i++)
 915     pool[i+sz-1]=_pool[i];
 916   _pool=pool;
 917   //
 918   if(sz==2)
 919     updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
 920   //
 921   if(isFast)
 922     {
 923       _ze_mesh=mesh;
 924       return ;
 925     }
 926   //
 927   std::vector< MCAuto<MEDCouplingUMesh> > ms;
 928   if(pos>0)
 929     {
 930       MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
 931       ms.push_back(elt);
 932     }
 933   ms.push_back(mesh);
 934   if(pos<_ze_mesh->getNumberOfCells()-1)
 935   {
 936     MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
 937     ms.push_back(elt);
 938   }
 939   std::vector< const MEDCouplingUMesh *> ms2(ms.size());
 940   for(std::size_t j=0;j<ms2.size();j++)
 941     ms2[j]=ms[j];
 942   _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
 943 }
 944
 945 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
 946 {
 947   _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
 948 }
 949
 950 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
 951 {
 952   if(pos<0)
 953     throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
 954   int ret(0);
 955   for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
 956     {
 957       if((*it).isInMyRange(pos))
 958         return ret;
 959     }
 960   throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
 961 }
 962
 963 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
 964 {
 965   get(pos);//to check;
 966   if(_edge_info.empty())
 967     return ;
 968   std::size_t sz(_edge_info.size()-1);
 969   for(std::size_t i=0;i<sz;i++)
 970     _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
 971 }
 972
 973 const CellInfo& VectorOfCellInfo::get(int pos) const
 974 {
 975   if(pos<0 || pos>=(int)_pool.size())
 976     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
 977   return _pool[pos];
 978 }
 979
 980 CellInfo& VectorOfCellInfo::get(int pos)
 981 {
 982   if(pos<0 || pos>=(int)_pool.size())
 983     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
 984   return _pool[pos];
 985 }
 986
 987 /*!
 988  * Given :
 989  * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
 990  * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
 991  *
 992  * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
 993  *
 994  * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
 995  *
 996  * \param [in] allEdges a list of pairs (beginNode, endNode). Linked with \a allEdgesPtr to get the equation of edge.
 997  */
 998 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
 999                                          MCAuto<DataArrayInt>& idsLeftRight)
1000 {
1001   int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1002   if(nbCellsInSplitMesh1D==0)
1003     throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1004   const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1005   std::size_t nb(allEdges.size()),jj;
1006   if(nb%2!=0)
1007     throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1008   std::vector<int> edge1Bis(nb*2);
1009   std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1010   std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1011   std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1012   std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1013   std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1014   //
1015   idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1016   int *idsLeftRightPtr(idsLeftRight->getPointer());
1017   VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1018   for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1019     {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1020       int iEnd(iStart);
1021       for(;iEnd<nbCellsInSplitMesh1D;)
1022         {
1023           for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1024           if(jj!=nb)
1025             break;
1026           else
1027             iEnd++;
1028         }
1029       if(iEnd<nbCellsInSplitMesh1D)
1030         iEnd++;
1031       //
1032       MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1033       int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1034       //
1035       MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1036       retTmp->setCoords(splitMesh1D->getCoords());
1037       retTmp->allocateCells();
1038
1039       std::vector< std::vector<int> > out0;
1040       std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1041
1042       BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1043       for(std::size_t cnt=0;cnt<out0.size();cnt++)
1044         AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1045       pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1046       //
1047       iStart=iEnd;
1048     }
1049   for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1050     pool.feedEdgeInfoAt(eps,mm,offset,idsLeftRightPtr+2*mm);
1051   return pool.getZeMesh().retn();
1052 }
1053
1054 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, const MEDCouplingUMesh *splitMesh1D,
1055                                      const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1056                                      MCAuto<DataArrayInt>& idsLeftRight)
1057 {
1058   const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1059   //
1060   std::vector<int> allEdges;
1061   std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1062   for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1063     {
1064       int edgeId(std::abs(*it)-1);
1065       std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1066       MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1067       const std::vector<int>& edge1(intersectEdge1[edgeId]);
1068       if(*it>0)
1069         allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1070       else
1071         allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1072       std::size_t sz(edge1.size());
1073       for(std::size_t cnt=0;cnt<sz;cnt++)
1074         allEdgesPtr.push_back(ee);
1075     }
1076   //
1077   return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1078 }
1079
1080 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1081 {
1082   if(!typ1.isQuadratic() && !typ2.isQuadratic())
1083     {//easy case comparison not
1084       return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1085     }
1086   else if(typ1.isQuadratic() && typ2.isQuadratic())
1087     {
1088       bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1089       if(!status0)
1090         return false;
1091       if(conn1[2]==conn2[2])
1092         return true;
1093       const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1094       double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1095       return dist<eps;
1096     }
1097   else
1098     {//only one is quadratic
1099       bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1100       if(!status0)
1101         return false;
1102       const double *a(0),*bb(0),*be(0);
1103       if(typ1.isQuadratic())
1104         {
1105           a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1106         }
1107       else
1108         {
1109           a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1110         }
1111       double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1112       double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1113       return dist<eps;
1114     }
1115 }
1116
1117 /*!
1118  * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1119  * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1120  *
1121  * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1122  */
1123 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1124 {
1125   if(candidatesIn2DEnd==candidatesIn2DBg)
1126     throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1127   const double *coo(mesh2DSplit->getCoords()->begin());
1128   if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1129     return *candidatesIn2DBg;
1130   int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1131   MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1132   if(cellIdInMesh1DSplitRelative<0)
1133     cur1D->changeOrientationOfCells();
1134   const int *c1D(cur1D->getNodalConnectivity()->begin());
1135   const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1136   for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1137     {
1138       MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1139       const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1140       const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1141       unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1142       INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1143       for(unsigned it2=0;it2<sz;it2++)
1144         {
1145           INTERP_KERNEL::NormalizedCellType typeOfSon;
1146           cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1147           const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1148           if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1149             return *it;
1150         }
1151     }
1152   throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1153 }
1154
1155 /*!
1156  * \param [out] intersectEdge1 - for each cell in \a m1Desc returns the result of the split. The result is given using pair of int given resp start and stop.
1157  *                               So for all edge \a i in \a m1Desc \a  intersectEdge1[i] is of length 2*n where n is the number of sub edges.
1158  *                               And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1159  * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1160  * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1161  * \param [out] addCoo - nodes to be append at the end
1162  * \param [out] mergedNodes - gives all pair of nodes of \a m2Desc that have same location than some nodes in \a m1Desc. key is id in \a m2Desc offseted and value is id in \a m1Desc.
1163  */
1164 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1165                                          std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2, std::vector<double>& addCoo, std::map<int,int>& mergedNodes)
1166 {
1167   static const int SPACEDIM=2;
1168   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1169   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1170   const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1171   // Build BB tree of all edges in the tool mesh (second mesh)
1172   MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2Desc->getBoundingBoxForBBTree());
1173   const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1174   int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1175   intersectEdge1.resize(nDescCell1);
1176   colinear2.resize(nDescCell2);
1177   subDiv2.resize(nDescCell2);
1178   BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1179
1180   std::vector<int> candidates1(1);
1181   int offset1(m1Desc->getNumberOfNodes());
1182   int offset2(offset1+m2Desc->getNumberOfNodes());
1183   for(int i=0;i<nDescCell1;i++)  // for all edges in the first mesh
1184     {
1185       std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1186       myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1187       if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1188         {
1189           std::map<INTERP_KERNEL::Node *,int> map1,map2;
1190           // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1191           INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1192           candidates1[0]=i;
1193           INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
1194           // 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
1195           // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1196           std::set<INTERP_KERNEL::Node *> nodes;
1197           pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1198           std::size_t szz(nodes.size());
1199           std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1200           std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1201           for(std::size_t iii=0;iii<szz;iii++,itt++)
1202             { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1203           // end of protection
1204           // Performs egde cutting:
1205           pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1206           delete pol2;
1207           delete pol1;
1208         }
1209       else
1210         // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1211         intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1212     }
1213 }
1214
1215
1216 /*!
1217  * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1218  * It builds the descending connectivity of the two meshes, and then using a binary tree
1219  * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1220  * Documentation about parameters  colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1221  */
1222 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1223                                                    std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1224                                                    MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1225                                                    std::vector<double>& addCoo,
1226                                                    MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1227 {
1228   // Build desc connectivity
1229   desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1230   desc2=DataArrayInt::New();
1231   descIndx2=DataArrayInt::New();
1232   revDesc2=DataArrayInt::New();
1233   revDescIndx2=DataArrayInt::New();
1234   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1235   MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1236   m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1237   m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1238   MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1239   std::map<int,int> notUsedMap;
1240   Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1241   m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1242   m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1243 }
1244
1245 /**
1246  * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1247  * (newly created) nodes corresponding to the edge intersections.
1248  * Output params:
1249  * @param[out] cr, crI connectivity of the resulting mesh
1250  * @param[out] cNb1, cNb2 correspondance arrays giving for the merged mesh the initial cells IDs in m1 / m2
1251  * TODO: describe input parameters
1252  */
1253 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1254                                                          const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1255                                                          const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1256                                                          const std::vector<double>& addCoords,
1257                                                          std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1258 {
1259   static const int SPACEDIM=2;
1260   const double *coo1(m1->getCoords()->begin());
1261   const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1262   int offset1(m1->getNumberOfNodes());
1263   const double *coo2(m2->getCoords()->begin());
1264   const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1265   int offset2(offset1+m2->getNumberOfNodes());
1266   int offset3(offset2+((int)addCoords.size())/2);
1267   MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
1268   const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1269   // Here a BBTree on 2D-cells, not on segments:
1270   BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1271   int ncell1(m1->getNumberOfCells());
1272   crI.push_back(0);
1273   for(int i=0;i<ncell1;i++)
1274     {
1275       std::vector<int> candidates2;
1276       myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1277       std::map<INTERP_KERNEL::Node *,int> mapp;
1278       std::map<int,INTERP_KERNEL::Node *> mappRev;
1279       INTERP_KERNEL::QuadraticPolygon pol1;
1280       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1281       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1282       // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1283       MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1284       // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
1285       pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1286           desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1287       //
1288       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
1289       std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1290       INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1291       for(it1.first();!it1.finished();it1.next())
1292         edges1.insert(it1.current()->getPtr());
1293       //
1294       std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1295       std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1296       int ii=0;
1297       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1298         {
1299           INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1300           const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1301           // Complete mapping with elements coming from the current cell it2 in mesh2:
1302           MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1303           // pol2 is the new QP in the final merged result.
1304           pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1305               pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1306         }
1307       ii=0;
1308       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1309         {
1310           INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1311           pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1312           //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1313           pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1314         }
1315       // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1316       // by m2 but that we still want to keep in the final result.
1317       if(!edges1.empty())
1318         {
1319           try
1320           {
1321               INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1322           }
1323           catch(INTERP_KERNEL::Exception& e)
1324           {
1325               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();
1326               throw INTERP_KERNEL::Exception(oss.str());
1327           }
1328         }
1329       for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1330         (*it).second->decrRef();
1331     }
1332 }
1333
1334 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1335 {
1336   std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1337   if(it!=conn.end())
1338     return ;
1339   std::size_t sz(conn.size());
1340   std::size_t found(std::numeric_limits<std::size_t>::max());
1341   for(std::size_t i=0;i<sz;i++)
1342     {
1343       int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1344       double v1[3]={coords[3*pt1+0]-coords[3*pt0+0],coords[3*pt1+1]-coords[3*pt0+1],coords[3*pt1+2]-coords[3*pt0+2]},v2[3]={coords[3*nodeIdToInsert+0]-coords[3*pt0+0],coords[3*nodeIdToInsert+1]-coords[3*pt0+1],coords[3*nodeIdToInsert+2]-coords[3*pt0+2]};
1345       double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1346       std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1347       std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1348       double v3[3];
1349       v3[0]=v1[1]*v2[2]-v1[2]*v2[1]; v3[1]=v1[2]*v2[0]-v1[0]*v2[2]; v3[2]=v1[0]*v2[1]-v1[1]*v2[0];
1350       double normm2(sqrt(v3[0]*v3[0]+v3[1]*v3[1]+v3[2]*v3[2])),dotTest(v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]);
1351       if(normm2<eps)
1352         if(dotTest>eps && dotTest<1.-eps)
1353           {
1354             found=i;
1355             break;
1356           }
1357     }
1358   if(found==std::numeric_limits<std::size_t>::max())
1359     throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1360   conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1361 }
1362
1363 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1364 {
1365   std::size_t sz(conn.size());
1366   std::vector<int> *curPart(&part0);
1367   for(std::size_t i=0;i<sz;i++)
1368     {
1369       int nextt(conn[(i+1)%sz]);
1370       (*curPart).push_back(nextt);
1371       if(nextt==pt0 || nextt==pt1)
1372         {
1373           if(curPart==&part0)
1374             curPart=&part1;
1375           else
1376             curPart=&part0;
1377           (*curPart).push_back(nextt);
1378         }
1379     }
1380 }
1381
1382 /*!
1383  * this method method splits cur cells 3D Surf in sub cells 3DSurf using the previous subsplit. This method is the last one used to clip.
1384  */
1385 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1386                                             const int *desc, const int *descIndx, const double *coords, double eps,
1387                                             std::vector<std::vector<int> >& res) const
1388 {
1389   checkFullyDefined();
1390   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1391     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1392   const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1393   int nbOfCells(getNumberOfCells());
1394   if(nbOfCells!=1)
1395     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1396   for(int i=0;i<nbOfCells;i++)
1397     {
1398       int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset),start(-1),end(-1);
1399       for(int j=0;j<nbOfFaces;j++)
1400         {
1401           const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1402           const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1403           int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1404           INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1405           INTERP_KERNEL::NormalizedCellType cmsId;
1406           unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1407           std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1408           if(p.first!=-1 && p.second!=-1)
1409             {
1410               if(p.first!=-2)
1411                 {
1412                   InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1413                   InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1414                   std::vector<int> elt1,elt2;
1415                   SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1416                   res.push_back(elt1);
1417                   res.push_back(elt2);
1418                 }
1419               else
1420                 res.push_back(elt);
1421             }
1422           else
1423             res.push_back(elt);
1424         }
1425     }
1426 }
1427
1428 /*!
1429  * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additionnal nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
1430  *
1431  * \sa MEDCouplingUMesh::split2DCells
1432  */
1433 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1434 {
1435   checkConnectivityFullyDefined();
1436   int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1437   MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1438   const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1439   int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1440   int prevPosOfCi(ciPtr[0]);
1441   for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1442     {
1443       int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1444       *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1445       for(int j=0;j<sz;j++)
1446         {
1447           int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1448           for(int k=0;k<sz2;k++)
1449             *cPtr++=subPtr[offset2+k];
1450           if(j!=sz-1)
1451             *cPtr++=oldConn[prevPosOfCi+j+2];
1452           deltaSz+=sz2;
1453         }
1454       prevPosOfCi=ciPtr[1];
1455       ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1456     }
1457   if(c->end()!=cPtr)
1458     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1459   _nodal_connec->decrRef();
1460   _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1461 }
1462
1463
1464 /*!
1465  * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additionnal nodes can be added at the end of coordinates array object.
1466  *
1467  * \return  int - the number of new nodes created.
1468  * \sa MEDCouplingUMesh::split2DCells
1469  */
1470 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1471 {
1472   checkConsistencyLight();
1473   int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1474   MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1475   MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1476   const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1477   const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1478   const double *oldCoordsPtr(getCoords()->begin());
1479   int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1480   int prevPosOfCi(ciPtr[0]);
1481   for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1482     {
1483       int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1484       for(int j=0;j<sz;j++)
1485         { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1486       *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1487       for(int j=0;j<sz;j++)//loop over subedges of oldConn
1488         {
1489           int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1490           if(sz2==0)
1491             {
1492               if(j<sz-1)
1493                 cPtr[1]=oldConn[prevPosOfCi+2+j];
1494               cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1495               continue;
1496             }
1497           std::vector<INTERP_KERNEL::Node *> ns(3);
1498           ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1499           ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1500           ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1501           MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1502           for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1503             {
1504               cPtr[1]=subPtr[offset2+k];
1505               cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1506             }
1507           int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1508           if(j!=sz-1)
1509             { cPtr[1]=tmpEnd; }
1510           cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1511         }
1512       prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1513       ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1514     }
1515   if(c->end()!=cPtr)
1516     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1517   _nodal_connec->decrRef();
1518   _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1519   addCoo->rearrange(2);
1520   MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1521   setCoords(coo);
1522   return addCoo->getNumberOfTuples();
1523 }
1524
1525
1526 /// @endcond
1527
1528 /*!
1529  * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1530  * returns a result mesh constituted by polygons.
1531  * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1532  * all nodes from m2.
1533  * The meshes should be in 2D space. In
1534  * addition, returns two arrays mapping cells of the result mesh to cells of the input
1535  * meshes.
1536  *  \param [in] m1 - the first input mesh which is a partitioned object. The mesh must be so that each point in the space covered by \a m1
1537  *                      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)
1538  *  \param [in] m2 - the second input mesh which is a partition tool. The mesh must be so that each point in the space covered by \a m2
1539  *                      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)
1540  *  \param [in] eps - precision used to detect coincident mesh entities.
1541  *  \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1542  *         cell an id of the cell of \a m1 it comes from. The caller is to delete
1543  *         this array using decrRef() as it is no more needed.
1544  *  \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1545  *         cell an id of the cell of \a m2 it comes from. -1 value means that a
1546  *         result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1547  *         any cell of \a m2. The caller is to delete this array using decrRef() as
1548  *         it is no more needed.
1549  *  \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1550  *         MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1551  *         is no more needed.
1552  *  \throw If the coordinates array is not set in any of the meshes.
1553  *  \throw If the nodal connectivity of cells is not defined in any of the meshes.
1554  *  \throw If any of the meshes is not a 2D mesh in 2D space.
1555  *
1556  *  \sa conformize2D, mergeNodes
1557  */
1558 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1559                                                       double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1560 {
1561   if(!m1 || !m2)
1562     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1563   m1->checkFullyDefined();
1564   m2->checkFullyDefined();
1565   if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1566     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2  with meshdim equal to 2 and spaceDim equal to 2 too!");
1567
1568   // Step 1: compute all edge intersections (new nodes)
1569   std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1570   MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1571   DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1572   std::vector<double> addCoo,addCoordsQuadratic;  // coordinates of newly created nodes
1573   IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1574                               m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1575                               addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1576   revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1577   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1578   MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1579
1580   // Step 2: re-order newly created nodes according to the ordering found in m2
1581   std::vector< std::vector<int> > intersectEdge2;
1582   BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1583   subDiv2.clear(); dd5=0; dd6=0;
1584
1585   // Step 3:
1586   std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1587   std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1588   BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1589                                     /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1590
1591   // Step 4: Prepare final result:
1592   MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1593   addCooDa->alloc((int)(addCoo.size())/2,2);
1594   std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1595   MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1596   addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1597   std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1598   std::vector<const DataArrayDouble *> coordss(4);
1599   coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1600   MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1601   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1602   MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1603   MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1604   MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1605   MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1606   ret->setConnectivity(conn,connI,true);
1607   ret->setCoords(coo);
1608   cellNb1=c1.retn(); cellNb2=c2.retn();
1609   return ret.retn();
1610 }
1611
1612 /*!
1613  * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1614  * 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
1615  * and finaly, in case of quadratic polygon the centers of edges new nodes.
1616  * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1617  *
1618  * \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
1619  *                      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)
1620  * \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
1621  *                      you can invoke orderConsecutiveCells1D on \a mesh1D.
1622  * \param [in] eps - precision used to perform intersections and localization operations.
1623  * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1624  * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1625  * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1626  *                               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.
1627  * \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
1628  *                               and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1629  *                               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.
1630  *
1631  * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1632  */
1633 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1634 {
1635   if(!mesh2D || !mesh1D)
1636     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1637   mesh2D->checkFullyDefined();
1638   mesh1D->checkFullyDefined();
1639   const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1640   if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1641     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1642   // Step 1: compute all edge intersections (new nodes)
1643   std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1644   std::vector<double> addCoo,addCoordsQuadratic;  // coordinates of newly created nodes
1645   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1646   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1647   //
1648   // Build desc connectivity
1649   DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1650   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1651   MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1652   std::map<int,int> mergedNodes;
1653   Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1654   // use mergeNodes to fix intersectEdge1
1655   for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1656     {
1657       std::size_t n((*it0).size()/2);
1658       int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1659       std::map<int,int>::const_iterator it1;
1660       it1=mergedNodes.find(eltStart);
1661       if(it1!=mergedNodes.end())
1662         (*it0)[0]=(*it1).second;
1663       it1=mergedNodes.find(eltEnd);
1664       if(it1!=mergedNodes.end())
1665         (*it0)[2*n-1]=(*it1).second;
1666     }
1667   //
1668   MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1669   addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
1670   // Step 2: re-order newly created nodes according to the ordering found in m2
1671   std::vector< std::vector<int> > intersectEdge2;
1672   BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1673   subDiv2.clear();
1674   // Step 3: compute splitMesh1D
1675   MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1676   MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1677   MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1678       idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1679   MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1680   MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1681   // deal with cells in mesh2D that are not cut but only some of their edges are
1682   MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1683   idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1684   idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1685   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
1686   if(!idsInDesc2DToBeRefined->empty())
1687     {
1688       DataArrayInt *out0(0),*outi0(0);
1689       MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1690       MCAuto<DataArrayInt> outi0s(outi0);
1691       out0s=out0;
1692       out0s=out0s->buildUnique();
1693       out0s->sort(true);
1694     }
1695   //
1696   MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1697   MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1698   MCAuto<DataArrayInt> elts,eltsIndex;
1699   mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1700   MCAuto<DataArrayInt> eltsIndex2(eltsIndex->deltaShiftIndex());
1701   MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1702   if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1703     throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1704   MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1705   MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1706   if((DataArrayInt *)out0s)
1707     untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1708   std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1709   // OK all is ready to insert in ret2 mesh
1710   if(!untouchedCells->empty())
1711     {// the most easy part, cells in mesh2D not impacted at all
1712       outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1713       outMesh2DSplit.back()->setCoords(ret1->getCoords());
1714       ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1715     }
1716   if((DataArrayInt *)out0s)
1717     {// here dealing with cells in out0s but not in cellsToBeModified
1718       MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1719       const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1720       for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1721         {
1722           outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1723           ret1->setCoords(outMesh2DSplit.back()->getCoords());
1724         }
1725       int offset(ret2->getNumberOfTuples());
1726       ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1727       MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1728       partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1729       int kk(0),*ret3ptr(partOfRet3->getPointer());
1730       for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1731         {
1732           int faceId(std::abs(*it)-1);
1733           for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1734             {
1735               int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1736               if(tmp!=-1)
1737                 {
1738                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1739                     ret3ptr[2*kk]=tmp+offset;
1740                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1741                     ret3ptr[2*kk+1]=tmp+offset;
1742                 }
1743               else
1744                 {//the current edge is shared by a 2D cell that will be split just after
1745                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1746                     ret3ptr[2*kk]=-(*it2+1);
1747                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1748                     ret3ptr[2*kk+1]=-(*it2+1);
1749                 }
1750             }
1751         }
1752       m1Desc->setCoords(ret1->getCoords());
1753       ret1NonCol->setCoords(ret1->getCoords());
1754       ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1755       if(!outMesh2DSplit.empty())
1756         {
1757           DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1758           for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1759             (*itt)->setCoords(da);
1760         }
1761     }
1762   cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1763   for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1764     {
1765       MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1766       idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1767       MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1768       MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1769       MCAuto<DataArrayInt> partOfRet3;
1770       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));
1771       ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1772       outMesh2DSplit.push_back(splitOfOneCell);
1773       for(int i=0;i<splitOfOneCell->getNumberOfCells();i++)
1774         ret2->pushBackSilent(*it);
1775     }
1776   //
1777   std::size_t nbOfMeshes(outMesh2DSplit.size());
1778   std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1779   for(std::size_t i=0;i<nbOfMeshes;i++)
1780     tmp[i]=outMesh2DSplit[i];
1781   //
1782   ret1->getCoords()->setInfoOnComponents(compNames);
1783   MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1784   // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1785   ret3->rearrange(1);
1786   MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStricltyNegative());
1787   for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1788     {
1789       int old2DCellId(-ret3->getIJ(*it,0)-1);
1790       MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1791       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
1792     }
1793   ret3->changeValue(std::numeric_limits<int>::max(),-1);
1794   ret3->rearrange(2);
1795   //
1796   splitMesh1D=ret1.retn();
1797   splitMesh2D=ret2D.retn();
1798   cellIdInMesh2D=ret2.retn();
1799   cellIdInMesh1D=ret3.retn();
1800 }
1801
1802 /*!
1803  * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1804  * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1805  * 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
1806  * 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).
1807  * 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.
1808  *
1809  * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1810  * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1811  *
1812  * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1813  * This method expects that all nodes in \a this are not closer than \a eps.
1814  * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1815  *
1816  * \param [in] eps the relative error to detect merged edges.
1817  * \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
1818  *                           that the user is expected to deal with.
1819  *
1820  * \throw If \a this is not coherent.
1821  * \throw If \a this has not spaceDim equal to 2.
1822  * \throw If \a this has not meshDim equal to 2.
1823  * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1824  */
1825 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1826 {
1827   static const int SPACEDIM=2;
1828   checkConsistencyLight();
1829   if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1830     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1831   MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1832   MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1833   const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1834   MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
1835   const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1836   int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1837   std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1838   std::vector<double> addCoo;
1839   BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1840   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1841   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1842   for(int i=0;i<nDescCell;i++)
1843     {
1844       std::vector<int> candidates;
1845       myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1846       for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1847         if(*it>i)
1848           {
1849             std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1850             INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1851                 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1852             INTERP_KERNEL::MergePoints merge;
1853             INTERP_KERNEL::QuadraticPolygon c1,c2;
1854             e1->intersectWith(e2,merge,c1,c2);
1855             e1->decrRef(); e2->decrRef();
1856             if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1857               overlapEdge[i].push_back(*it);
1858             if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1859               overlapEdge[*it].push_back(i);
1860           }
1861     }
1862   // splitting done. sort intersect point in intersectEdge.
1863   std::vector< std::vector<int> > middle(nDescCell);
1864   int nbOf2DCellsToBeSplit(0);
1865   bool middleNeedsToBeUsed(false);
1866   std::vector<bool> cells2DToTreat(nDescCell,false);
1867   for(int i=0;i<nDescCell;i++)
1868     {
1869       std::vector<int>& isect(intersectEdge[i]);
1870       int sz((int)isect.size());
1871       if(sz>1)
1872         {
1873           std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1874           INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1875           e->sortSubNodesAbs(coords,isect);
1876           e->decrRef();
1877         }
1878       if(sz!=0)
1879         {
1880           int idx0(rdi[i]),idx1(rdi[i+1]);
1881           if(idx1-idx0!=1)
1882             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1883           if(!cells2DToTreat[rd[idx0]])
1884             {
1885               cells2DToTreat[rd[idx0]]=true;
1886               nbOf2DCellsToBeSplit++;
1887             }
1888           // try to reuse at most eventual 'middle' of SEG3
1889           std::vector<int>& mid(middle[i]);
1890           mid.resize(sz+1,-1);
1891           if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1892             {
1893               middleNeedsToBeUsed=true;
1894               const std::vector<int>& candidates(overlapEdge[i]);
1895               std::vector<int> trueCandidates;
1896               for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1897                 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1898                   trueCandidates.push_back(*itc);
1899               int stNode(c[ci[i]+1]),endNode(isect[0]);
1900               for(int j=0;j<sz+1;j++)
1901                 {
1902                   for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1903                     {
1904                       int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1905                       if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1906                         { mid[j]=*itc; break; }
1907                     }
1908                   stNode=endNode;
1909                   endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1910                 }
1911             }
1912         }
1913     }
1914   MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1915   if(nbOf2DCellsToBeSplit==0)
1916     return ret.retn();
1917   //
1918   int *retPtr(ret->getPointer());
1919   for(int i=0;i<nCell;i++)
1920     if(cells2DToTreat[i])
1921       *retPtr++=i;
1922   //
1923   MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1924   DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1925   MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1926   DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1927   if(middleNeedsToBeUsed)
1928     { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1929   MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1930   int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1931   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.
1932   setPartOfMySelf(ret->begin(),ret->end(),*modif);
1933   {
1934     bool areNodesMerged; int newNbOfNodes;
1935     if(nbOfNodesCreated!=0)
1936       MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
1937   }
1938   return ret.retn();
1939 }
1940
1941 /*!
1942  * 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.
1943  * 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).
1944  * 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
1945  * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
1946  * 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
1947  * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
1948  *
1949  * 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
1950  * using new instance, idem for coordinates.
1951  *
1952  * 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.
1953  *
1954  * \return DataArrayInt  * - The list of cellIds in \a this that have at least one edge colinearized.
1955  *
1956  * \throw If \a this is not coherent.
1957  * \throw If \a this has not spaceDim equal to 2.
1958  * \throw If \a this has not meshDim equal to 2.
1959  *
1960  * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
1961  */
1962 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
1963 {
1964   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
1965   checkConsistencyLight();
1966   if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1967     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1968   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
1969   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
1970   int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
1971   const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
1972   MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
1973   MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
1974   const double *coords(_coords->begin());
1975   int *newciptr(newci->getPointer());
1976   for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
1977     {
1978       if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,newc,appendedCoords))
1979         ret->pushBackSilent(i);
1980       newciptr[1]=newc->getNumberOfTuples();
1981     }
1982   //
1983   if(ret->empty())
1984     return ret.retn();
1985   if(!appendedCoords->empty())
1986     {
1987       appendedCoords->rearrange(2);
1988       MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
1989       //non const part
1990       setCoords(newCoords);
1991     }
1992   //non const part
1993   setConnectivity(newc,newci,true);
1994   return ret.retn();
1995 }
1996
1997