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 "MEDCouplingFieldDouble.hxx"
  24 #include "CellModel.hxx"
  25 #include "VolSurfUser.txx"
  26 #include "InterpolationUtils.hxx"
  27 #include "PointLocatorAlgos.txx"
  28 #include "BBTree.txx"
  29 #include "BBTreeDst.txx"
  30 #include "DirectedBoundingBox.hxx"
  31 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  32 #include "InterpKernelAutoPtr.hxx"
  33 #include "InterpKernelGeo2DNode.hxx"
  34 #include "InterpKernelGeo2DEdgeLin.hxx"
  35 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  36 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
  37 #include "TranslationRotationMatrix.hxx"
  38 #include "VectorUtils.hxx"
  39 #include "MEDCouplingSkyLineArray.hxx"
  40
  41 #include <sstream>
  42 #include <fstream>
  43 #include <numeric>
  44 #include <cstring>
  45 #include <limits>
  46 #include <list>
  47 #include <assert.h>
  48
  49 using namespace MEDCoupling;
  50
  51 /// @cond INTERNAL
  52
  53 int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
  54 {
  55   if(id!=-1)
  56     return id;
  57   else
  58     {
  59       int ret(nodesCnter++);
  60       double newPt[2];
  61       e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt);
  62       addCoo.insertAtTheEnd(newPt,newPt+2);
  63       return ret;
  64     }
  65 }
  66
  67 int InternalAddPointOriented(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
  68 {
  69   if(id!=-1)
  70     return id;
  71   else
  72     {
  73       int ret(nodesCnter++);
  74       double newPt[2];
  75       e->getMiddleOfPointsOriented(coo+2*startId,coo+2*endId,newPt);
  76       addCoo.insertAtTheEnd(newPt,newPt+2);
  77       return ret;
  78     }
  79 }
  80
  81
  82 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)
  83 {
  84   int tmp[3];
  85   int trueStart(start>=0?start:nbOfEdges+start);
  86   tmp[0]=linOrArc?(int)INTERP_KERNEL::NORM_QPOLYG:(int)INTERP_KERNEL::NORM_POLYGON; tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp];
  87   newConnOfCell->insertAtTheEnd(tmp,tmp+3);
  88   if(linOrArc)
  89     {
  90       if(stp-start>1)
  91         {
  92           int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
  93           InternalAddPointOriented(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2);
  94           middles.push_back(tmp3+offset);
  95         }
  96       else
  97         middles.push_back(connBg[trueStart+nbOfEdges]);
  98     }
  99 }
 100
 101 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)
 102 {
 103   int tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]);
 104   newConnOfCell->pushBackSilent(tmpEnd);
 105   if(linOrArc)
 106     {
 107       if(stp-start>1)
 108         {
 109           int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
 110           InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
 111           middles.push_back(tmp3+offset);
 112         }
 113       else
 114         middles.push_back(connBg[start+nbOfEdges]);
 115     }
 116 }
 117
 118 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)
 119 {
 120   // only the quadratic point to deal with:
 121   if(linOrArc)
 122     {
 123       if(stp-start>1)  // if we are covering more than one segment we need to create a new mid point
 124         {
 125           int tmpSrt(connBg[start]),tmpEnd(connBg[stp % nbOfEdges]);  // % to handle last seg.
 126           int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
 127           InternalAddPointOriented(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
 128           middles.push_back(tmp3+offset);
 129         }
 130       else
 131         middles.push_back(connBg[start+nbOfEdges]);
 132     }
 133 }
 134
 135 void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
 136 {
 137   MCAuto<INTERP_KERNEL::Node> nTmp(n); nTmp->incrRef();
 138   std::map<MCAuto<INTERP_KERNEL::Node>,int>::const_iterator it(m.find(nTmp));
 139   if(it==m.end())
 140     throw INTERP_KERNEL::Exception("Internal error in remapping !");
 141   int v((*it).second);
 142   if(v==forbVal0 || v==forbVal1)
 143     return ;
 144   if(std::find(isect.begin(),isect.end(),v)==isect.end())
 145     isect.push_back(v);
 146 }
 147
 148 bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map<MCAuto<INTERP_KERNEL::Node>,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
 149 {
 150   int sz(c.size());
 151   if(sz<=1)
 152     return false;
 153   bool presenceOfOn(false);
 154   for(int i=0;i<sz;i++)
 155     {
 156       INTERP_KERNEL::ElementaryEdge *e(c[i]);
 157       if(e->getLoc()!=INTERP_KERNEL::FULL_ON_1)
 158         continue ;
 159       IKGeo2DInternalMapper2(e->getStartNode(),m,forbVal0,forbVal1,isect);
 160       IKGeo2DInternalMapper2(e->getEndNode(),m,forbVal0,forbVal1,isect);
 161     }
 162   return presenceOfOn;
 163 }
 164
 165
 166 namespace MEDCoupling
 167 {
 168
 169   INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge2(INTERP_KERNEL::NormalizedCellType typ, const int *bg, const double *coords2D, std::map< MCAuto<INTERP_KERNEL::Node>,int>& m)
 170   {
 171     INTERP_KERNEL::Edge *ret(0);
 172     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]));
 173     m[n0]=bg[0]; m[n1]=bg[1];
 174     switch(typ)
 175     {
 176       case INTERP_KERNEL::NORM_SEG2:
 177         {
 178           ret=new INTERP_KERNEL::EdgeLin(n0,n1);
 179           break;
 180         }
 181       case INTERP_KERNEL::NORM_SEG3:
 182         {
 183           INTERP_KERNEL::Node *n2(new INTERP_KERNEL::Node(coords2D[2*bg[2]],coords2D[2*bg[2]+1])); m[n2]=bg[2];
 184           INTERP_KERNEL::EdgeLin *e1(new INTERP_KERNEL::EdgeLin(n0,n2)),*e2(new INTERP_KERNEL::EdgeLin(n2,n1));
 185           INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
 186           // is the SEG3 degenerated, and thus can be reduced to a SEG2?
 187           bool colinearity(inters.areColinears());
 188           delete e1; delete e2;
 189           if(colinearity)
 190             { ret=new INTERP_KERNEL::EdgeLin(n0,n1); }
 191           else
 192             { ret=new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1); }
 193           break;
 194         }
 195       default:
 196         throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge2 : Expecting a mesh with spaceDim==2 and meshDim==1 !");
 197     }
 198     return ret;
 199   }
 200
 201   INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >& mapp2, const int *bg)
 202   {
 203     INTERP_KERNEL::Edge *ret=0;
 204     switch(typ)
 205     {
 206       case INTERP_KERNEL::NORM_SEG2:
 207         {
 208           ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
 209           break;
 210         }
 211       case INTERP_KERNEL::NORM_SEG3:
 212         {
 213           INTERP_KERNEL::EdgeLin *e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[2]].first);
 214           INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]].first,mapp2[bg[1]].first);
 215           INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
 216           // is the SEG3 degenerated, and thus can be reduced to a SEG2?
 217           bool colinearity=inters.areColinears();
 218           delete e1; delete e2;
 219           if(colinearity)
 220             ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
 221           else
 222             ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]].first,mapp2[bg[2]].first,mapp2[bg[1]].first);
 223           mapp2[bg[2]].second=false;
 224           break;
 225         }
 226       default:
 227         throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !");
 228     }
 229     return ret;
 230   }
 231
 232   /*!
 233    * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed by the sub set of cells in 'candidates' taken from
 234    * the global mesh 'mDesc'.
 235    * The input mesh 'mDesc' must be so that mDim==1 and spaceDim==2.
 236    * 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'.
 237    */
 238   INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates,
 239                                                                    std::map<INTERP_KERNEL::Node *,int>& mapp)
 240   {
 241     mapp.clear();
 242     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.
 243     const double *coo=mDesc->getCoords()->getConstPointer();
 244     const int *c=mDesc->getNodalConnectivity()->getConstPointer();
 245     const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
 246     std::set<int> s;
 247     for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
 248       s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
 249     for(std::set<int>::const_iterator it2=s.begin();it2!=s.end();it2++)
 250       {
 251         INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]);
 252         mapp2[*it2]=std::pair<INTERP_KERNEL::Node *,bool>(n,true);
 253       }
 254     INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon;
 255     for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
 256       {
 257         INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]];
 258         ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1));
 259       }
 260     for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++)
 261       {
 262         if((*it2).second.second)
 263           mapp[(*it2).second.first]=(*it2).first;
 264         ((*it2).second.first)->decrRef();
 265       }
 266     return ret;
 267   }
 268
 269   INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo)
 270   {
 271     if(nodeId>=offset2)
 272       {
 273         int locId=nodeId-offset2;
 274         return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
 275       }
 276     if(nodeId>=offset1)
 277       {
 278         int locId=nodeId-offset1;
 279         return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]);
 280       }
 281     return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]);
 282   }
 283
 284   /**
 285    * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI).
 286    */
 287   void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
 288                                         const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
 289                                         /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
 290   {
 291     for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
 292       {
 293         int eltId1=abs(*desc1)-1;
 294         for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
 295           {
 296             std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.find(*it1);
 297             if(it==mappRev.end())
 298               {
 299                 INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo);
 300                 mapp[node]=*it1;
 301                 mappRev[*it1]=node;
 302               }
 303           }
 304       }
 305   }
 306 }
 307
 308
 309
 310 /*!
 311  * 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 ) .
 312  * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
 313  * \param forbiddenPoints the list of points that should not be removed in the process
 314  */
 315 bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset,
 316                                          const std::map<int, bool>& forbiddenPoints,
 317                                          DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
 318 {
 319   std::size_t sz(std::distance(connBg,connEnd));
 320   if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
 321     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
 322   sz--;
 323   INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
 324   INTERP_KERNEL::AutoPtr<int> tmpConn2(new int[sz]);
 325   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
 326   unsigned nbs(cm.getNumberOfSons2(connBg+1,sz));
 327   unsigned nbOfHit(0); // number of fusions operated
 328   int posBaseElt(0),posEndElt(0),nbOfTurn(0);
 329   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
 330   INTERP_KERNEL::NormalizedCellType typeOfSon;
 331   std::vector<int> middles;
 332   bool ret(false);
 333   for(;(nbOfTurn+nbOfHit)<nbs;nbOfTurn++)
 334     {
 335       cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
 336       std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
 337       INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
 338       posEndElt = posBaseElt+1;
 339
 340       // Look backward first: are the final edges of the cells colinear with the first ones?
 341       // This initializes posBaseElt.
 342       if(nbOfTurn==0)
 343         {
 344           for(unsigned i=1;i<nbs && nbOfHit<maxNbOfHit;i++) // 2nd condition is to avoid ending with a cell with one single edge
 345             {
 346               cm.fillSonCellNodalConnectivity2(nbs-i,connBg+1,sz,tmpConn2,typeOfSon);
 347               // Identify common point:
 348               int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
 349               auto itE(forbiddenPoints.end());
 350               if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
 351                 break;
 352               INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
 353               INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
 354               bool isColinear=eint->areColinears();
 355               if(isColinear)
 356                 {
 357                   nbOfHit++;
 358                   posBaseElt--;
 359                   ret=true;
 360                 }
 361               delete eint;
 362               eCand->decrRef();
 363               if(!isColinear)
 364                 break;
 365               // Update last connectivity
 366               std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
 367             }
 368         }
 369       // Now move forward:
 370       const unsigned fwdStart = (nbOfTurn == 0 ? 0 : posBaseElt);  // the first element to be inspected going forward
 371       for(unsigned j=fwdStart+1;j<nbs && nbOfHit<maxNbOfHit;j++)  // 2nd condition is to avoid ending with a cell with one single edge
 372         {
 373           cm.fillSonCellNodalConnectivity2((int)j,connBg+1,sz,tmpConn2,typeOfSon); // get edge #j's connectivity
 374           // Identify common point:
 375           int commPoint = std::find((int *)tmpConn, tmpConn+2, tmpConn2[0]) != tmpConn+2 ? tmpConn2[0] : tmpConn2[1];
 376           auto itE(forbiddenPoints.end());
 377           if (forbiddenPoints.find(commPoint) != itE) // is the junction point in the list of points we can not remove?
 378             break;
 379           INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn2,coords,m));
 380           INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
 381           bool isColinear(eint->areColinears());
 382           if(isColinear)
 383             {
 384               nbOfHit++;
 385               posEndElt++;
 386               ret=true;
 387             }
 388           delete eint;
 389           eCand->decrRef();
 390           if(!isColinear)
 391               break;
 392           // Update last connectivity
 393           std::copy((int *)tmpConn2, tmpConn2+sz, (int *)tmpConn);
 394         }
 395       //push [posBaseElt,posEndElt) in newConnOfCell using e
 396       // The if clauses below are (voluntary) not mutually exclusive: on a quad cell with 2 edges, the end of the connectivity is also its beginning!
 397       if(nbOfTurn==0)
 398         // at the beginning of the connectivity (insert type)
 399         EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
 400       else if((nbOfHit+nbOfTurn) != (nbs-1))
 401         // in the middle
 402         EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
 403       if ((nbOfHit+nbOfTurn) == (nbs-1))
 404         // at the end (only quad points to deal with)
 405         EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
 406       posBaseElt=posEndElt;
 407       e->decrRef();
 408     }
 409   if(!middles.empty())
 410     newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
 411   return ret;
 412 }
 413
 414
 415
 416 bool IsColinearOfACellOf(const std::vector< std::vector<int> >& intersectEdge1, const std::vector<int>& candidates, int start, int stop, int& retVal)
 417 {
 418   if(candidates.empty())
 419     return false;
 420   for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
 421     {
 422       const std::vector<int>& pool(intersectEdge1[*it]);
 423       int tmp[2]; tmp[0]=start; tmp[1]=stop;
 424       if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
 425         {
 426           retVal=*it+1;
 427           return true;
 428         }
 429       tmp[0]=stop; tmp[1]=start;
 430       if(std::search(pool.begin(),pool.end(),tmp,tmp+2)!=pool.end())
 431         {
 432           retVal=-*it-1;
 433           return true;
 434         }
 435     }
 436   return false;
 437 }
 438
 439 /*!
 440  * This method performs the 2nd step of Partition of 2D mesh.
 441  * This method has 4 inputs :
 442  *  - a mesh 'm1' with meshDim==1 and a SpaceDim==2
 443  *  - a mesh 'm2' with meshDim==1 and a SpaceDim==2
 444  *  - subDiv of size 'm2->getNumberOfCells()' that lists for each seg cell in 'm2' the splitting node ids randomly sorted.
 445  * 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'
 446  * Nodes end up lying consecutively on a cutted edge.
 447  * \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.
 448  * (Only present for its coords in case of 'subDiv' shares some nodes of 'm1')
 449  * \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.
 450  * \param addCoo input parameter with additional nodes linked to intersection of the 2 meshes.
 451  * \param[out] intersectEdge the same content as subDiv, but correctly oriented.
 452  */
 453 void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
 454                                            const std::vector<double>& addCoo,
 455                                            const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
 456 {
 457   int offset1=m1->getNumberOfNodes();
 458   int ncell=m2->getNumberOfCells();
 459   const int *c=m2->getNodalConnectivity()->begin();
 460   const int *cI=m2->getNodalConnectivityIndex()->begin();
 461   const double *coo=m2->getCoords()->begin();
 462   const double *cooBis=m1->getCoords()->begin();
 463   int offset2=offset1+m2->getNumberOfNodes();
 464   intersectEdge.resize(ncell);
 465   for(int i=0;i<ncell;i++,cI++)
 466     {
 467       const std::vector<int>& divs=subDiv[i];
 468       int nnode=cI[1]-cI[0]-1;
 469       std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;
 470       std::map<INTERP_KERNEL::Node *, int> mapp22;
 471       for(int j=0;j<nnode;j++)
 472         {
 473           INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
 474           int nnid=c[(*cI)+j+1];
 475           mapp2[nnid]=std::pair<INTERP_KERNEL::Node *,bool>(nn,true);
 476           mapp22[nn]=nnid+offset1;
 477         }
 478       INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
 479       for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
 480         ((*it).second.first)->decrRef();
 481       std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
 482       std::map<INTERP_KERNEL::Node *,int> mapp3;
 483       for(std::size_t j=0;j<divs.size();j++)
 484         {
 485           int id=divs[j];
 486           INTERP_KERNEL::Node *tmp=0;
 487           if(id<offset1)
 488             tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
 489           else if(id<offset2)
 490             tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
 491           else
 492             tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
 493           addNodes[j]=tmp;
 494           mapp3[tmp]=id;
 495         }
 496       e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
 497       for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
 498         (*it)->decrRef();
 499       e->decrRef();
 500     }
 501 }
 502
 503 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,
 504                                      MCAuto<DataArrayInt>& idsInRetColinear, MCAuto<DataArrayInt>& idsInMesh1DForIdsInRetColinear)
 505 {
 506   idsInRetColinear=DataArrayInt::New(); idsInRetColinear->alloc(0,1);
 507   idsInMesh1DForIdsInRetColinear=DataArrayInt::New(); idsInMesh1DForIdsInRetColinear->alloc(0,1);
 508   int nCells(mesh1D->getNumberOfCells());
 509   if(nCells!=(int)intersectEdge2.size())
 510     throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
 511   const DataArrayDouble *coo2(mesh1D->getCoords());
 512   const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
 513   const double *coo2Ptr(coo2->begin());
 514   int offset1(coords1->getNumberOfTuples());
 515   int offset2(offset1+coo2->getNumberOfTuples());
 516   int offset3(offset2+addCoo.size()/2);
 517   std::vector<double> addCooQuad;
 518   MCAuto<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
 519   int tmp[4],cicnt(0),kk(0);
 520   for(int i=0;i<nCells;i++)
 521     {
 522       std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
 523       INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
 524       const std::vector<int>& subEdges(intersectEdge2[i]);
 525       int nbSubEdge(subEdges.size()/2);
 526       for(int j=0;j<nbSubEdge;j++,kk++)
 527         {
 528           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));
 529           MCAuto<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
 530           INTERP_KERNEL::Edge *e2Ptr(e2);
 531           std::map<int,int>::const_iterator itm;
 532           if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
 533             {
 534               tmp[0]=INTERP_KERNEL::NORM_SEG3;
 535               itm=mergedNodes.find(subEdges[2*j]);
 536               tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
 537               itm=mergedNodes.find(subEdges[2*j+1]);
 538               tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
 539               tmp[3]=offset3+(int)addCooQuad.size()/2;
 540               double tmp2[2];
 541               e2->getBarycenter(tmp2); addCooQuad.insert(addCooQuad.end(),tmp2,tmp2+2);
 542               cicnt+=4;
 543               cOut->insertAtTheEnd(tmp,tmp+4);
 544               ciOut->pushBackSilent(cicnt);
 545             }
 546           else
 547             {
 548               tmp[0]=INTERP_KERNEL::NORM_SEG2;
 549               itm=mergedNodes.find(subEdges[2*j]);
 550               tmp[1]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j];
 551               itm=mergedNodes.find(subEdges[2*j+1]);
 552               tmp[2]=itm!=mergedNodes.end()?(*itm).second:subEdges[2*j+1];
 553               cicnt+=3;
 554               cOut->insertAtTheEnd(tmp,tmp+3);
 555               ciOut->pushBackSilent(cicnt);
 556             }
 557           int tmp00;
 558           if(IsColinearOfACellOf(intersectEdge1,colinear2[i],tmp[1],tmp[2],tmp00))
 559             {
 560               idsInRetColinear->pushBackSilent(kk);
 561               idsInMesh1DForIdsInRetColinear->pushBackSilent(tmp00);
 562             }
 563         }
 564       e->decrRef();
 565     }
 566   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
 567   ret->setConnectivity(cOut,ciOut,true);
 568   MCAuto<DataArrayDouble> arr3(DataArrayDouble::New());
 569   arr3->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
 570   MCAuto<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,C_DEALLOC,(int)addCooQuad.size()/2,2);
 571   std::vector<const DataArrayDouble *> coordss(4);
 572   coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
 573   MCAuto<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
 574   ret->setCoords(arr);
 575   return ret.retn();
 576 }
 577
 578 MEDCouplingUMesh *BuildRefined2DCellLinear(const DataArrayDouble *coords, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
 579 {
 580   std::vector<int> allEdges;
 581   for(const int *it2(descBg);it2!=descEnd;it2++)
 582     {
 583       const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
 584       if(*it2>0)
 585         allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
 586       else
 587         allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
 588     }
 589   std::size_t nb(allEdges.size());
 590   if(nb%2!=0)
 591     throw INTERP_KERNEL::Exception("BuildRefined2DCellLinear : internal error 1 !");
 592   std::size_t nbOfEdgesOf2DCellSplit(nb/2);
 593   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
 594   ret->setCoords(coords);
 595   ret->allocateCells(1);
 596   std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
 597   for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
 598     connOut[kk]=allEdges[2*kk];
 599   ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,connOut.size(),&connOut[0]);
 600   return ret.retn();
 601 }
 602
 603 MEDCouplingUMesh *BuildRefined2DCellQuadratic(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
 604 {
 605   const int *c(mesh2D->getNodalConnectivity()->begin()),*ci(mesh2D->getNodalConnectivityIndex()->begin());
 606   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[cellIdInMesh2D]]));
 607   std::size_t ii(0);
 608   unsigned sz(cm.getNumberOfSons2(c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1));
 609   if(sz!=std::distance(descBg,descEnd))
 610     throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 1 !");
 611   INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]]);
 612   std::vector<int> allEdges,centers;
 613   const double *coordsPtr(coords->begin());
 614   MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
 615   int offset(coords->getNumberOfTuples());
 616   for(const int *it2(descBg);it2!=descEnd;it2++,ii++)
 617     {
 618       INTERP_KERNEL::NormalizedCellType typeOfSon;
 619       cm.fillSonCellNodalConnectivity2(ii,c+ci[cellIdInMesh2D]+1,ci[cellIdInMesh2D+1]-ci[cellIdInMesh2D]-1,tmpPtr,typeOfSon);
 620       const std::vector<int>& edge1(intersectEdge1[std::abs(*it2)-1]);
 621       if(*it2>0)
 622         allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
 623       else
 624         allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
 625       if(edge1.size()==2)
 626         centers.push_back(tmpPtr[2]);//special case where no subsplit of edge -> reuse the original center.
 627       else
 628         {//the current edge has been subsplit -> create corresponding centers.
 629           std::size_t nbOfCentersToAppend(edge1.size()/2);
 630           std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
 631           MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpPtr,coordsPtr,m));
 632           std::vector<int>::const_iterator it3(allEdges.end()-edge1.size());
 633           for(std::size_t k=0;k<nbOfCentersToAppend;k++)
 634             {
 635               double tmpp[2];
 636               const double *aa(coordsPtr+2*(*it3++));
 637               const double *bb(coordsPtr+2*(*it3++));
 638               ee->getMiddleOfPoints(aa,bb,tmpp);
 639               addCoo->insertAtTheEnd(tmpp,tmpp+2);
 640               centers.push_back(offset+k);
 641             }
 642         }
 643     }
 644   std::size_t nb(allEdges.size());
 645   if(nb%2!=0)
 646     throw INTERP_KERNEL::Exception("BuildRefined2DCellQuadratic : internal error 2 !");
 647   std::size_t nbOfEdgesOf2DCellSplit(nb/2);
 648   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("",2));
 649   if(addCoo->empty())
 650     ret->setCoords(coords);
 651   else
 652     {
 653       addCoo->rearrange(2);
 654       addCoo=DataArrayDouble::Aggregate(coords,addCoo);
 655       ret->setCoords(addCoo);
 656     }
 657   ret->allocateCells(1);
 658   std::vector<int> connOut(nbOfEdgesOf2DCellSplit);
 659   for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
 660     connOut[kk]=allEdges[2*kk];
 661   connOut.insert(connOut.end(),centers.begin(),centers.end());
 662   ret->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,connOut.size(),&connOut[0]);
 663   return ret.retn();
 664 }
 665
 666 /*!
 667  * This method creates a refinement of a cell in \a mesh2D. Those cell is defined by descending connectivity and the sorted subdivided nodal connectivity
 668  * of those edges.
 669  *
 670  * \param [in] mesh2D - The origin 2D mesh. \b Warning \b coords are not those of \a mesh2D. But mesh2D->getCoords()==coords[:mesh2D->getNumberOfNodes()]
 671  */
 672 MEDCouplingUMesh *BuildRefined2DCell(const DataArrayDouble *coords, const MEDCouplingUMesh *mesh2D, int cellIdInMesh2D, const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1)
 673 {
 674   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(mesh2D->getTypeOfCell(cellIdInMesh2D)));
 675   if(!cm.isQuadratic())
 676     return BuildRefined2DCellLinear(coords,descBg,descEnd,intersectEdge1);
 677   else
 678     return BuildRefined2DCellQuadratic(coords,mesh2D,cellIdInMesh2D,descBg,descEnd,intersectEdge1);
 679 }
 680
 681 void AddCellInMesh2D(MEDCouplingUMesh *mesh2D, const std::vector<int>& conn, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edges)
 682 {
 683   bool isQuad(false);
 684   for(std::vector< MCAuto<INTERP_KERNEL::Edge> >::const_iterator it=edges.begin();it!=edges.end();it++)
 685     {
 686       const INTERP_KERNEL::Edge *ee(*it);
 687       if(dynamic_cast<const INTERP_KERNEL::EdgeArcCircle *>(ee))
 688         isQuad=true;
 689     }
 690   if(!isQuad)
 691     mesh2D->insertNextCell(INTERP_KERNEL::NORM_POLYGON,conn.size(),&conn[0]);
 692   else
 693     {
 694       const double *coo(mesh2D->getCoords()->begin());
 695       std::size_t sz(conn.size());
 696       std::vector<double> addCoo;
 697       std::vector<int> conn2(conn);
 698       int offset(mesh2D->getNumberOfNodes());
 699       for(std::size_t i=0;i<sz;i++)
 700         {
 701           double tmp[2];
 702           edges[(i+1)%sz]->getMiddleOfPoints(coo+2*conn[i],coo+2*conn[(i+1)%sz],tmp);// tony a chier i+1 -> i
 703           addCoo.insert(addCoo.end(),tmp,tmp+2);
 704           conn2.push_back(offset+(int)i);
 705         }
 706       mesh2D->getCoords()->rearrange(1);
 707       mesh2D->getCoords()->pushBackValsSilent(&addCoo[0],&addCoo[0]+addCoo.size());
 708       mesh2D->getCoords()->rearrange(2);
 709       mesh2D->insertNextCell(INTERP_KERNEL::NORM_QPOLYG,conn2.size(),&conn2[0]);
 710     }
 711 }
 712
 713 /*!
 714  * \b WARNING edges in out1 coming from \a splitMesh1D are \b NOT oriented because only used for equation of curve.
 715  *
 716  * This method cuts in 2 parts the input 2D cell given using boundaries description (\a edge1Bis and \a edge1BisPtr) using
 717  * a set of edges defined in \a splitMesh1D.
 718  */
 719 void BuildMesh2DCutInternal2(const MEDCouplingUMesh *splitMesh1D, const std::vector<int>& edge1Bis, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edge1BisPtr,
 720                              std::vector< std::vector<int> >& out0, std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& out1)
 721 {
 722   std::size_t nb(edge1Bis.size()/2);
 723   std::size_t nbOfEdgesOf2DCellSplit(nb/2);
 724   int iEnd(splitMesh1D->getNumberOfCells());
 725   if(iEnd==0)
 726     throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal2 : internal error ! input 1D mesh must have at least one cell !");
 727   std::size_t ii,jj;
 728   const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
 729   for(ii=0;ii<nb && edge1Bis[2*ii]!=cSplitPtr[ciSplitPtr[0]+1];ii++);
 730   for(jj=ii;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd-1]+2];jj++);
 731   //
 732   if(jj==nb)
 733     {//the edges splitMesh1D[iStart:iEnd] does not fully cut the current 2D cell -> single output cell
 734       out0.resize(1); out1.resize(1);
 735       std::vector<int>& connOut(out0[0]);
 736       connOut.resize(nbOfEdgesOf2DCellSplit);
 737       std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr(out1[0]);
 738       edgesPtr.resize(nbOfEdgesOf2DCellSplit);
 739       for(std::size_t kk=0;kk<nbOfEdgesOf2DCellSplit;kk++)
 740         {
 741           connOut[kk]=edge1Bis[2*kk];
 742           edgesPtr[kk]=edge1BisPtr[2*kk];
 743         }
 744     }
 745   else
 746     {
 747       // [i,iEnd[ contains the
 748       out0.resize(2); out1.resize(2);
 749       std::vector<int>& connOutLeft(out0[0]);
 750       std::vector<int>& connOutRight(out0[1]);//connOutLeft should end with edge1Bis[2*ii] and connOutRight should end with edge1Bis[2*jj+1]
 751       std::vector< MCAuto<INTERP_KERNEL::Edge> >& eleft(out1[0]);
 752       std::vector< MCAuto<INTERP_KERNEL::Edge> >& eright(out1[1]);
 753       for(std::size_t k=ii;k<jj+1;k++)
 754         { connOutLeft.push_back(edge1Bis[2*k+1]); eleft.push_back(edge1BisPtr[2*k+1]); }
 755       std::vector< MCAuto<INTERP_KERNEL::Edge> > ees(iEnd);
 756       for(int ik=0;ik<iEnd;ik++)
 757         {
 758           std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
 759           MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cSplitPtr[ciSplitPtr[ik]],cSplitPtr+ciSplitPtr[ik]+1,splitMesh1D->getCoords()->begin(),m));
 760           ees[ik]=ee;
 761         }
 762       for(int ik=iEnd-1;ik>=0;ik--)
 763         connOutLeft.push_back(cSplitPtr[ciSplitPtr[ik]+1]);
 764       for(std::size_t k=jj+1;k<nbOfEdgesOf2DCellSplit+ii;k++)
 765         { connOutRight.push_back(edge1Bis[2*k+1]); eright.push_back(edge1BisPtr[2*k+1]); }
 766       eleft.insert(eleft.end(),ees.rbegin(),ees.rend());
 767       for(int ik=0;ik<iEnd;ik++)
 768         connOutRight.push_back(cSplitPtr[ciSplitPtr[ik]+2]);
 769       eright.insert(eright.end(),ees.begin(),ees.end());
 770     }
 771 }
 772
 773 struct CellInfo
 774 {
 775 public:
 776   CellInfo() { }
 777   CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
 778 public:
 779   std::vector<int> _edges;
 780   std::vector< MCAuto<INTERP_KERNEL::Edge> > _edges_ptr;
 781 };
 782
 783 CellInfo::CellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr)
 784 {
 785   std::size_t nbe(edges.size());
 786   std::vector<int> edges2(2*nbe); std::vector< MCAuto<INTERP_KERNEL::Edge> > edgesPtr2(2*nbe);
 787   for(std::size_t i=0;i<nbe;i++)
 788     {
 789       edges2[2*i]=edges[i]; edges2[2*i+1]=edges[(i+1)%nbe];
 790       edgesPtr2[2*i]=edgesPtr[(i+1)%nbe]; edgesPtr2[2*i+1]=edgesPtr[(i+1)%nbe];//tony a chier
 791     }
 792   _edges.resize(4*nbe); _edges_ptr.resize(4*nbe);
 793   std::copy(edges2.begin(),edges2.end(),_edges.begin()); std::copy(edges2.begin(),edges2.end(),_edges.begin()+2*nbe);
 794   std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()); std::copy(edgesPtr2.begin(),edgesPtr2.end(),_edges_ptr.begin()+2*nbe);
 795 }
 796
 797 class EdgeInfo
 798 {
 799 public:
 800   EdgeInfo(int istart, int iend, const MCAuto<MEDCouplingUMesh>& mesh):_istart(istart),_iend(iend),_mesh(mesh),_left(-7),_right(-7) { }
 801   EdgeInfo(int istart, int iend, int pos, const MCAuto<INTERP_KERNEL::Edge>& edge):_istart(istart),_iend(iend),_edge(edge),_left(pos),_right(pos+1) { }
 802   bool isInMyRange(int pos) const { return pos>=_istart && pos<_iend; }
 803   void somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
 804   void feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const;
 805 private:
 806   int _istart;
 807   int _iend;
 808   MCAuto<MEDCouplingUMesh> _mesh;
 809   MCAuto<INTERP_KERNEL::Edge> _edge;
 810   int _left;
 811   int _right;
 812 };
 813
 814 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
 815 {
 816   const MEDCouplingUMesh *mesh(_mesh);
 817   if(mesh)
 818     return ;
 819   if(_right<pos)
 820     return ;
 821   if(_left>pos)
 822     { _left++; _right++; return ; }
 823   if(_right==pos)
 824     {
 825       bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
 826       if((isLeft && isRight) || (!isLeft && !isRight))
 827         throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
 828       if(isLeft)
 829         return ;
 830       if(isRight)
 831         {
 832           _right++;
 833           return ;
 834         }
 835     }
 836   if(_left==pos)
 837     {
 838       bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
 839       if((isLeft && isRight) || (!isLeft && !isRight))
 840         throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
 841       if(isLeft)
 842         {
 843           _right++;
 844           return ;
 845         }
 846       if(isRight)
 847         {
 848           _left++;
 849           _right++;
 850           return ;
 851         }
 852     }
 853 }
 854
 855 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
 856 {
 857   const MEDCouplingUMesh *mesh(_mesh);
 858   if(!mesh)
 859     {
 860       neighbors[0]=offset+_left; neighbors[1]=offset+_right;
 861     }
 862   else
 863     {// not fully splitting cell case
 864       if(mesh2D->getNumberOfCells()==1)
 865         {//little optimization. 1 cell no need to find in which cell mesh is !
 866           neighbors[0]=offset; neighbors[1]=offset;
 867           return;
 868         }
 869       else
 870         {
 871           MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
 872           int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
 873           if(cellId==-1)
 874             throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
 875           neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
 876         }
 877     }
 878 }
 879
 880 class VectorOfCellInfo
 881 {
 882 public:
 883   VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
 884   std::size_t size() const { return _pool.size(); }
 885   int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
 886   void setMeshAt(std::size_t 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);
 887   const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
 888   const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
 889   MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
 890   void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
 891 private:
 892   int getZePosOfEdgeGivenItsGlobalId(int pos) const;
 893   void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
 894   const CellInfo& get(int pos) const;
 895   CellInfo& get(int pos);
 896 private:
 897   std::vector<CellInfo> _pool;
 898   MCAuto<MEDCouplingUMesh> _ze_mesh;
 899   std::vector<EdgeInfo> _edge_info;
 900 };
 901
 902 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
 903 {
 904   _pool[0]._edges=edges;
 905   _pool[0]._edges_ptr=edgesPtr;
 906 }
 907
 908 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
 909 {
 910   if(_pool.empty())
 911     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
 912   if(_pool.size()==1)
 913     return 0;
 914   const MEDCouplingUMesh *zeMesh(_ze_mesh);
 915   if(!zeMesh)
 916     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
 917   MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
 918   return zeMesh->getCellContainingPoint(barys->begin(),eps);
 919 }
 920
 921 void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend,
 922                                  const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges,
 923                                  const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
 924 {
 925   get(pos);//to check pos
 926   bool isFast(pos==0 && _pool.size()==1);
 927   std::size_t sz(edges.size());
 928   // dealing with edges
 929   if(sz==1)
 930     _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
 931   else
 932     _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
 933   //
 934   std::vector<CellInfo> pool(_pool.size()-1+sz);
 935   for(std::size_t i=0;i<pos;i++)
 936     pool[i]=_pool[i];
 937   for(std::size_t j=0;j<sz;j++)
 938     pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
 939   for(int i=pos+1;i<(int)_pool.size();i++)
 940     pool[i+sz-1]=_pool[i];
 941   _pool=pool;
 942   //
 943   if(sz==2)
 944     updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
 945   //
 946   if(isFast)
 947     {
 948       _ze_mesh=mesh;
 949       return ;
 950     }
 951   //
 952   std::vector< MCAuto<MEDCouplingUMesh> > ms;
 953   if(pos>0)
 954     {
 955       MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
 956       ms.push_back(elt);
 957     }
 958   ms.push_back(mesh);
 959   if(pos<_ze_mesh->getNumberOfCells()-1)
 960   {
 961     MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
 962     ms.push_back(elt);
 963   }
 964   std::vector< const MEDCouplingUMesh *> ms2(ms.size());
 965   for(std::size_t j=0;j<ms2.size();j++)
 966     ms2[j]=ms[j];
 967   _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
 968 }
 969
 970 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
 971 {
 972   _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
 973 }
 974
 975 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
 976 {
 977   if(pos<0)
 978     throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
 979   int ret(0);
 980   for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
 981     {
 982       if((*it).isInMyRange(pos))
 983         return ret;
 984     }
 985   throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
 986 }
 987
 988 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
 989 {
 990   get(pos);//to perform the sanity check;
 991   if(_edge_info.empty())
 992     return ;
 993   std::size_t sz(_edge_info.size()-1);
 994   for(std::size_t i=0;i<sz;i++)
 995     _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
 996 }
 997
 998 const CellInfo& VectorOfCellInfo::get(int pos) const
 999 {
1000   if(pos<0 || pos>=(int)_pool.size())
1001     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
1002   return _pool[pos];
1003 }
1004
1005 CellInfo& VectorOfCellInfo::get(int pos)
1006 {
1007   if(pos<0 || pos>=(int)_pool.size())
1008     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
1009   return _pool[pos];
1010 }
1011
1012 /*!
1013  * Given :
1014  * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
1015  * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
1016  *
1017  * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1018  *
1019  * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1020  *
1021  * \param [in] allEdges a list of pairs (beginNode, endNode). Represents all edges (already cut) in the single 2D cell being handled here. Linked with \a allEdgesPtr to get the equation of edge.
1022  */
1023 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
1024                                          MCAuto<DataArrayInt>& idsLeftRight)
1025 {
1026   int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1027   if(nbCellsInSplitMesh1D==0)
1028     throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1029   const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1030   std::size_t nb(allEdges.size()),jj;
1031   if(nb%2!=0)
1032     throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1033   std::vector<int> edge1Bis(nb*2);
1034   std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1035   std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1036   std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1037   std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1038   std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1039   //
1040   idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1041   int *idsLeftRightPtr(idsLeftRight->getPointer());
1042   VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1043
1044   // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1045   // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1046   // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1047   // of the connectivity.
1048   MCAuto <DataArrayInt> renumb(DataArrayInt::New());
1049   renumb->alloc(nbCellsInSplitMesh1D,1);
1050   const int * renumbP(renumb->begin());
1051
1052   int i, first=cSplitPtr[1];
1053   // Follow 1D line backward as long as it is connected:
1054   for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1055     first=cSplitPtr[ciSplitPtr[i]+1];
1056   if (i < nbCellsInSplitMesh1D-1)
1057     {
1058       // Build circular permutation to shift consecutive edges together
1059       renumb->iota(i+1);
1060       renumb->applyModulus(nbCellsInSplitMesh1D);
1061       splitMesh1D->renumberCells(renumbP, false);
1062       cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1063       ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1064     }
1065   else
1066     renumb->iota();
1067   //
1068
1069   for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1070     {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1071       int iEnd(iStart);
1072       for(;iEnd<nbCellsInSplitMesh1D;)
1073         {
1074           for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1075           if(jj!=nb)
1076             break;
1077           else
1078             iEnd++;
1079         }
1080       if(iEnd<nbCellsInSplitMesh1D)
1081         iEnd++;
1082
1083       MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1084       int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1085       //
1086       MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1087       retTmp->setCoords(splitMesh1D->getCoords());
1088       retTmp->allocateCells();
1089
1090       std::vector< std::vector<int> > out0;
1091       std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1092
1093       BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1094       for(std::size_t cnt=0;cnt<out0.size();cnt++)
1095         AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1096       pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1097       //
1098       iStart=iEnd;
1099     }
1100   for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1101     pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1102
1103   return pool.getZeMesh().retn();
1104 }
1105
1106 /*
1107  * splitMesh1D is an input parameter but might have its cells renumbered.
1108  */
1109 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1110                                      const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1111                                      MCAuto<DataArrayInt>& idsLeftRight)
1112 {
1113   const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1114   //
1115   std::vector<int> allEdges;
1116   std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1117   for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1118     {
1119       int edgeId(std::abs(*it)-1);
1120       std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1121       MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1122       const std::vector<int>& edge1(intersectEdge1[edgeId]);
1123       if(*it>0)
1124         allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1125       else
1126         allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1127       std::size_t sz(edge1.size());
1128       for(std::size_t cnt=0;cnt<sz;cnt++)
1129         allEdgesPtr.push_back(ee);
1130     }
1131   //
1132   return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1133 }
1134
1135 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1136 {
1137   if(!typ1.isQuadratic() && !typ2.isQuadratic())
1138     {//easy case comparison not
1139       return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1140     }
1141   else if(typ1.isQuadratic() && typ2.isQuadratic())
1142     {
1143       bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1144       if(!status0)
1145         return false;
1146       if(conn1[2]==conn2[2])
1147         return true;
1148       const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1149       double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1150       return dist<eps;
1151     }
1152   else
1153     {//only one is quadratic
1154       bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1155       if(!status0)
1156         return false;
1157       const double *a(0),*bb(0),*be(0);
1158       if(typ1.isQuadratic())
1159         {
1160           a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1161         }
1162       else
1163         {
1164           a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1165         }
1166       double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1167       double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1168       return dist<eps;
1169     }
1170 }
1171
1172 /*!
1173  * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1174  * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1175  *
1176  * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1177  */
1178 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1179 {
1180   if(candidatesIn2DEnd==candidatesIn2DBg)
1181     throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1182   const double *coo(mesh2DSplit->getCoords()->begin());
1183   if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1184     return *candidatesIn2DBg;
1185   int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1186   MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1187   if(cellIdInMesh1DSplitRelative<0)
1188     cur1D->changeOrientationOfCells();
1189   const int *c1D(cur1D->getNodalConnectivity()->begin());
1190   const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1191   for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1192     {
1193       MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1194       const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1195       const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1196       unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1197       INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1198       for(unsigned it2=0;it2<sz;it2++)
1199         {
1200           INTERP_KERNEL::NormalizedCellType typeOfSon;
1201           cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1202           const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1203           if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1204             return *it;
1205         }
1206     }
1207   throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1208 }
1209
1210 /*!
1211  * \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.
1212  *                               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.
1213  *                               And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1214  * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1215  * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1216  * \param [out] addCoo - nodes to be append at the end
1217  * \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 offsetted and value is id in \a m1Desc.
1218  */
1219 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1220                                          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)
1221 {
1222   static const int SPACEDIM=2;
1223   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1224   const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1225   // Build BB tree of all edges in the tool mesh (second mesh)
1226   MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1227   const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1228   int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1229   intersectEdge1.resize(nDescCell1);
1230   colinear2.resize(nDescCell2);
1231   subDiv2.resize(nDescCell2);
1232   BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1233
1234   std::vector<int> candidates1(1);
1235   int offset1(m1Desc->getNumberOfNodes());
1236   int offset2(offset1+m2Desc->getNumberOfNodes());
1237   for(int i=0;i<nDescCell1;i++)  // for all edges in the first mesh
1238     {
1239       std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1240       myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1241       if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1242         {
1243           std::map<INTERP_KERNEL::Node *,int> map1,map2;
1244           // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1245           INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1246           candidates1[0]=i;
1247           INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
1248           // 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
1249           // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1250           std::set<INTERP_KERNEL::Node *> nodes;
1251           pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1252           std::size_t szz(nodes.size());
1253           std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1254           std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1255           for(std::size_t iii=0;iii<szz;iii++,itt++)
1256             { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1257           // end of protection
1258           // Performs edge cutting:
1259           pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1260           delete pol2;
1261           delete pol1;
1262         }
1263       else
1264         // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1265         intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1266     }
1267 }
1268
1269
1270 /*!
1271  * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1272  * It builds the descending connectivity of the two meshes, and then using a binary tree
1273  * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1274  * Documentation about parameters  colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1275  */
1276 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1277                                                    std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1278                                                    MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1279                                                    std::vector<double>& addCoo,
1280                                                    MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1281 {
1282   // Build desc connectivity
1283   desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1284   desc2=DataArrayInt::New();
1285   descIndx2=DataArrayInt::New();
1286   revDesc2=DataArrayInt::New();
1287   revDescIndx2=DataArrayInt::New();
1288   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1289   MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1290   m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1291   m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1292   MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1293   std::map<int,int> notUsedMap;
1294   Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1295   m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1296   m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1297 }
1298
1299 /**
1300  * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1301  * (newly created) nodes corresponding to the edge intersections.
1302  * Output params:
1303  * @param[out] cr, crI connectivity of the resulting mesh
1304  * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1305  * TODO: describe input parameters
1306  */
1307 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1308                                                          const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1309                                                          const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1310                                                          const std::vector<double>& addCoords,
1311                                                          std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1312 {
1313   static const int SPACEDIM=2;
1314   const double *coo1(m1->getCoords()->begin());
1315   const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1316   int offset1(m1->getNumberOfNodes());
1317   const double *coo2(m2->getCoords()->begin());
1318   const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1319   int offset2(offset1+m2->getNumberOfNodes());
1320   int offset3(offset2+((int)addCoords.size())/2);
1321   MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1322   const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1323   // Here a BBTree on 2D-cells, not on segments:
1324   BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1325   int ncell1(m1->getNumberOfCells());
1326   crI.push_back(0);
1327   for(int i=0;i<ncell1;i++)
1328     {
1329       std::vector<int> candidates2;
1330       myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1331       std::map<INTERP_KERNEL::Node *,int> mapp;
1332       std::map<int,INTERP_KERNEL::Node *> mappRev;
1333       INTERP_KERNEL::QuadraticPolygon pol1;
1334       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1335       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1336       // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1337       MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1338       // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
1339       pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1340           desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1341       //
1342       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
1343       std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1344       INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1345       for(it1.first();!it1.finished();it1.next())
1346         edges1.insert(it1.current()->getPtr());
1347       //
1348       std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1349       std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1350       int ii=0;
1351       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1352         {
1353           INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1354           const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1355           // Complete mapping with elements coming from the current cell it2 in mesh2:
1356           MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1357           // pol2 is the new QP in the final merged result.
1358           pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1359               pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1360         }
1361       ii=0;
1362       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1363         {
1364           INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1365           pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1366           //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1367           pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1368         }
1369       // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1370       // by m2 but that we still want to keep in the final result.
1371       if(!edges1.empty())
1372         {
1373           try
1374           {
1375               INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1376           }
1377           catch(INTERP_KERNEL::Exception& e)
1378           {
1379               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();
1380               throw INTERP_KERNEL::Exception(oss.str());
1381           }
1382         }
1383       for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1384         (*it).second->decrRef();
1385     }
1386 }
1387
1388 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1389 {
1390   std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1391   if(it!=conn.end())
1392     return ;
1393   std::size_t sz(conn.size());
1394   std::size_t found(std::numeric_limits<std::size_t>::max());
1395   for(std::size_t i=0;i<sz;i++)
1396     {
1397       int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1398       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]};
1399       double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1400       std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1401       std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1402       double v3[3];
1403       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];
1404       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]);
1405       if(normm2<eps)
1406         if(dotTest>eps && dotTest<1.-eps)
1407           {
1408             found=i;
1409             break;
1410           }
1411     }
1412   if(found==std::numeric_limits<std::size_t>::max())
1413     throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1414   conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1415 }
1416
1417 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1418 {
1419   std::size_t sz(conn.size());
1420   std::vector<int> *curPart(&part0);
1421   for(std::size_t i=0;i<sz;i++)
1422     {
1423       int nextt(conn[(i+1)%sz]);
1424       (*curPart).push_back(nextt);
1425       if(nextt==pt0 || nextt==pt1)
1426         {
1427           if(curPart==&part0)
1428             curPart=&part1;
1429           else
1430             curPart=&part0;
1431           (*curPart).push_back(nextt);
1432         }
1433     }
1434 }
1435
1436 /*!
1437  * 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.
1438  */
1439 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1440                                             const int *desc, const int *descIndx, const double *coords, double eps,
1441                                             std::vector<std::vector<int> >& res) const
1442 {
1443   checkFullyDefined();
1444   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1445     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1446   const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1447   int nbOfCells(getNumberOfCells());
1448   if(nbOfCells!=1)
1449     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1450   for(int i=0;i<nbOfCells;i++)
1451     {
1452       int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1453       for(int j=0;j<nbOfFaces;j++)
1454         {
1455           const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1456           const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1457           int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1458           INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1459           INTERP_KERNEL::NormalizedCellType cmsId;
1460           unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1461           std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1462           if(p.first!=-1 && p.second!=-1)
1463             {
1464               if(p.first!=-2)
1465                 {
1466                   InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1467                   InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1468                   std::vector<int> elt1,elt2;
1469                   SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1470                   res.push_back(elt1);
1471                   res.push_back(elt2);
1472                 }
1473               else
1474                 res.push_back(elt);
1475             }
1476           else
1477             res.push_back(elt);
1478         }
1479     }
1480 }
1481
1482 /*!
1483  * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additional nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
1484  *
1485  * \sa MEDCouplingUMesh::split2DCells
1486  */
1487 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1488 {
1489   checkConnectivityFullyDefined();
1490   int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1491   MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1492   const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1493   int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1494   int prevPosOfCi(ciPtr[0]);
1495   for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1496     {
1497       int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1498       *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1499       for(int j=0;j<sz;j++)
1500         {
1501           int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1502           for(int k=0;k<sz2;k++)
1503             *cPtr++=subPtr[offset2+k];
1504           if(j!=sz-1)
1505             *cPtr++=oldConn[prevPosOfCi+j+2];
1506           deltaSz+=sz2;
1507         }
1508       prevPosOfCi=ciPtr[1];
1509       ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1510     }
1511   if(c->end()!=cPtr)
1512     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1513   _nodal_connec->decrRef();
1514   _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1515 }
1516
1517
1518 /*!
1519  * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1520  *
1521  * \return  int - the number of new nodes created.
1522  * \sa MEDCouplingUMesh::split2DCells
1523  */
1524 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1525 {
1526   checkConsistencyLight();
1527   int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1528   MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1529   MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1530   const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1531   const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1532   const double *oldCoordsPtr(getCoords()->begin());
1533   int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1534   int prevPosOfCi(ciPtr[0]);
1535   for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1536     {
1537       int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1538       for(int j=0;j<sz;j++)
1539         { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1540       *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1541       for(int j=0;j<sz;j++)//loop over subedges of oldConn
1542         {
1543           int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1544           if(sz2==0)
1545             {
1546               if(j<sz-1)
1547                 cPtr[1]=oldConn[prevPosOfCi+2+j];
1548               cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1549               continue;
1550             }
1551           std::vector<INTERP_KERNEL::Node *> ns(3);
1552           ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1553           ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1554           ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1555           MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1556           for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1557             {
1558               cPtr[1]=subPtr[offset2+k];
1559               cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1560             }
1561           int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1562           if(j!=sz-1)
1563             { cPtr[1]=tmpEnd; }
1564           cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1565         }
1566       prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1567       ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1568     }
1569   if(c->end()!=cPtr)
1570     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1571   _nodal_connec->decrRef();
1572   _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1573   addCoo->rearrange(2);
1574   MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1575   setCoords(coo);
1576   return addCoo->getNumberOfTuples();
1577 }
1578
1579
1580 /// @endcond
1581
1582 /*!
1583  * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1584  * returns a result mesh constituted by polygons.
1585  * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1586  * all nodes from m2.
1587  * The meshes should be in 2D space. In
1588  * addition, returns two arrays mapping cells of the result mesh to cells of the input
1589  * meshes.
1590  *  \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
1591  *                      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)
1592  *  \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
1593  *                      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)
1594  *  \param [in] eps - precision used to detect coincident mesh entities.
1595  *  \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1596  *         cell an id of the cell of \a m1 it comes from. The caller is to delete
1597  *         this array using decrRef() as it is no more needed.
1598  *  \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1599  *         cell an id of the cell of \a m2 it comes from. -1 value means that a
1600  *         result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1601  *         any cell of \a m2. The caller is to delete this array using decrRef() as
1602  *         it is no more needed.
1603  *  \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1604  *         MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1605  *         is no more needed.
1606  *  \throw If the coordinates array is not set in any of the meshes.
1607  *  \throw If the nodal connectivity of cells is not defined in any of the meshes.
1608  *  \throw If any of the meshes is not a 2D mesh in 2D space.
1609  *
1610  *  \sa conformize2D, mergeNodes
1611  */
1612 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1613                                                       double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1614 {
1615   if(!m1 || !m2)
1616     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1617   m1->checkFullyDefined();
1618   m2->checkFullyDefined();
1619   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1620   if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1621     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2  with meshdim equal to 2 and spaceDim equal to 2 too!");
1622
1623   // Step 1: compute all edge intersections (new nodes)
1624   std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1625   MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1626   DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1627   std::vector<double> addCoo,addCoordsQuadratic;  // coordinates of newly created nodes
1628   IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1629                               m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1630                               addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1631   revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1632   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1633   MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1634
1635   // Step 2: re-order newly created nodes according to the ordering found in m2
1636   std::vector< std::vector<int> > intersectEdge2;
1637   BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1638   subDiv2.clear(); dd5=0; dd6=0;
1639
1640   // Step 3:
1641   std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1642   std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1643   BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1644                                     /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1645
1646   // Step 4: Prepare final result:
1647   MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1648   addCooDa->alloc((int)(addCoo.size())/2,2);
1649   std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1650   MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1651   addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1652   std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1653   std::vector<const DataArrayDouble *> coordss(4);
1654   coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1655   MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1656   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1657   MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1658   MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1659   MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1660   MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1661   ret->setConnectivity(conn,connI,true);
1662   ret->setCoords(coo);
1663   cellNb1=c1.retn(); cellNb2=c2.retn();
1664   return ret.retn();
1665 }
1666
1667 /*!
1668  * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1669  * 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
1670  * and finally, in case of quadratic polygon the centers of edges new nodes.
1671  * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1672  *
1673  * \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
1674  *                      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)
1675  * \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
1676  *                      you can invoke orderConsecutiveCells1D on \a mesh1D.
1677  * \param [in] eps - precision used to perform intersections and localization operations.
1678  * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1679  * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1680  * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1681  *                               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.
1682  * \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
1683  *                               and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1684  *                               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.
1685  *
1686  * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1687  */
1688 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1689 {
1690   if(!mesh2D || !mesh1D)
1691     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1692   mesh2D->checkFullyDefined();
1693   mesh1D->checkFullyDefined();
1694   const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1695   if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1696     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1697   // Step 1: compute all edge intersections (new nodes)
1698   std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1699   std::vector<double> addCoo,addCoordsQuadratic;  // coordinates of newly created nodes
1700   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1701   //
1702   // Build desc connectivity
1703   DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1704   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1705   MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1706   std::map<int,int> mergedNodes;
1707   Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1708   // use mergeNodes to fix intersectEdge1
1709   for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1710     {
1711       std::size_t n((*it0).size()/2);
1712       int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1713       std::map<int,int>::const_iterator it1;
1714       it1=mergedNodes.find(eltStart);
1715       if(it1!=mergedNodes.end())
1716         (*it0)[0]=(*it1).second;
1717       it1=mergedNodes.find(eltEnd);
1718       if(it1!=mergedNodes.end())
1719         (*it0)[2*n-1]=(*it1).second;
1720     }
1721   //
1722   MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1723   addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
1724   // Step 2: re-order newly created nodes according to the ordering found in m2
1725   std::vector< std::vector<int> > intersectEdge2;
1726   BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1727   subDiv2.clear();
1728   // Step 3: compute splitMesh1D
1729   MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1730   MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1731   MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1732       idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1733   MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1734   MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1735   // deal with cells in mesh2D that are not cut but only some of their edges are
1736   MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1737   idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1738   idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1739   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
1740   if(!idsInDesc2DToBeRefined->empty())
1741     {
1742       DataArrayInt *out0(0),*outi0(0);
1743       MEDCouplingUMesh::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1744       MCAuto<DataArrayInt> outi0s(outi0);
1745       out0s=out0;
1746       out0s=out0s->buildUnique();
1747       out0s->sort(true);
1748     }
1749   //
1750   MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1751   MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1752   MCAuto<DataArrayInt> elts,eltsIndex;
1753   mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1754   MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1755   if (eltsIndex->getNumberOfTuples() > 1)
1756     eltsIndex2 = eltsIndex->deltaShiftIndex();
1757   MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1758   if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1759     throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1760   MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1761   MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1762   if((DataArrayInt *)out0s)
1763     untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1764   std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1765   // OK all is ready to insert in ret2 mesh
1766   if(!untouchedCells->empty())
1767     {// the most easy part, cells in mesh2D not impacted at all
1768       outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1769       outMesh2DSplit.back()->setCoords(ret1->getCoords());
1770       ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1771     }
1772   if((DataArrayInt *)out0s)
1773     {// here dealing with cells in out0s but not in cellsToBeModified
1774       MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1775       const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1776       for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1777         {
1778           outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1779           ret1->setCoords(outMesh2DSplit.back()->getCoords());
1780         }
1781       int offset(ret2->getNumberOfTuples());
1782       ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1783       MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1784       partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1785       int kk(0),*ret3ptr(partOfRet3->getPointer());
1786       for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1787         {
1788           int faceId(std::abs(*it)-1);
1789           for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1790             {
1791               int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1792               if(tmp!=-1)
1793                 {
1794                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1795                     ret3ptr[2*kk]=tmp+offset;
1796                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1797                     ret3ptr[2*kk+1]=tmp+offset;
1798                 }
1799               else
1800                 {//the current edge is shared by a 2D cell that will be split just after
1801                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1802                     ret3ptr[2*kk]=-(*it2+1);
1803                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1804                     ret3ptr[2*kk+1]=-(*it2+1);
1805                 }
1806             }
1807         }
1808       m1Desc->setCoords(ret1->getCoords());
1809       ret1NonCol->setCoords(ret1->getCoords());
1810       ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1811       if(!outMesh2DSplit.empty())
1812         {
1813           DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1814           for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1815             (*itt)->setCoords(da);
1816         }
1817     }
1818   cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1819   for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1820     {
1821       MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1822       idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1823       MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1824       MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1825       MCAuto<DataArrayInt> partOfRet3;
1826       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));
1827       ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1828       outMesh2DSplit.push_back(splitOfOneCell);
1829       for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1830         ret2->pushBackSilent(*it);
1831     }
1832   //
1833   std::size_t nbOfMeshes(outMesh2DSplit.size());
1834   std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1835   for(std::size_t i=0;i<nbOfMeshes;i++)
1836     tmp[i]=outMesh2DSplit[i];
1837   //
1838   ret1->getCoords()->setInfoOnComponents(compNames);
1839   MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1840   // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1841   ret3->rearrange(1);
1842   MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1843   for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1844     {
1845       int old2DCellId(-ret3->getIJ(*it,0)-1);
1846       MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1847       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
1848     }
1849   ret3->changeValue(std::numeric_limits<int>::max(),-1);
1850   ret3->rearrange(2);
1851   //
1852   splitMesh1D=ret1.retn();
1853   splitMesh2D=ret2D.retn();
1854   cellIdInMesh2D=ret2.retn();
1855   cellIdInMesh1D=ret3.retn();
1856 }
1857
1858 /*!
1859  * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1860  * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1861  * This method performs a conformization of \b this. So if a edge in \a this can be split into entire edges in \a this method
1862  * 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).
1863  * 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 invocation of this method.
1864  *
1865  * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1866  * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1867  *
1868  * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1869  * This method expects that all nodes in \a this are not closer than \a eps.
1870  * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1871  *
1872  * \param [in] eps the relative error to detect merged edges.
1873  * \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
1874  *                           that the user is expected to deal with.
1875  *
1876  * \throw If \a this is not coherent.
1877  * \throw If \a this has not spaceDim equal to 2.
1878  * \throw If \a this has not meshDim equal to 2.
1879  * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1880  */
1881 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1882 {
1883   static const int SPACEDIM=2;
1884   checkConsistencyLight();
1885   if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1886     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1887   MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1888   MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1889   const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1890   MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1891   const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1892   int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1893   std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1894   std::vector<double> addCoo;
1895   BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1896   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1897   for(int i=0;i<nDescCell;i++)
1898     {
1899       std::vector<int> candidates;
1900       myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1901       for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1902         if(*it>i)  // we're dealing with pair of edges, no need to treat the same pair twice
1903           {
1904             std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1905             INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1906                 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1907             INTERP_KERNEL::MergePoints merge;
1908             INTERP_KERNEL::QuadraticPolygon c1,c2;
1909             e1->intersectWith(e2,merge,c1,c2);
1910             e1->decrRef(); e2->decrRef();
1911             if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1912               overlapEdge[i].push_back(*it);
1913             if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1914               overlapEdge[*it].push_back(i);
1915           }
1916     }
1917   // splitting done. sort intersect point in intersectEdge.
1918   std::vector< std::vector<int> > middle(nDescCell);
1919   int nbOf2DCellsToBeSplit(0);
1920   bool middleNeedsToBeUsed(false);
1921   std::vector<bool> cells2DToTreat(nDescCell,false);
1922   for(int i=0;i<nDescCell;i++)
1923     {
1924       std::vector<int>& isect(intersectEdge[i]);
1925       int sz((int)isect.size());
1926       if(sz>1)
1927         {
1928           std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1929           INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
1930           e->sortSubNodesAbs(coords,isect);
1931           e->decrRef();
1932         }
1933       if(sz!=0)
1934         {
1935           int idx0(rdi[i]),idx1(rdi[i+1]);
1936           if(idx1-idx0!=1)
1937             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
1938           if(!cells2DToTreat[rd[idx0]])
1939             {
1940               cells2DToTreat[rd[idx0]]=true;
1941               nbOf2DCellsToBeSplit++;
1942             }
1943           // try to reuse at most eventual 'middle' of SEG3
1944           std::vector<int>& mid(middle[i]);
1945           mid.resize(sz+1,-1);
1946           if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
1947             {
1948               middleNeedsToBeUsed=true;
1949               const std::vector<int>& candidates(overlapEdge[i]);
1950               std::vector<int> trueCandidates;
1951               for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
1952                 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
1953                   trueCandidates.push_back(*itc);
1954               int stNode(c[ci[i]+1]),endNode(isect[0]);
1955               for(int j=0;j<sz+1;j++)
1956                 {
1957                   for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
1958                     {
1959                       int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
1960                       if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
1961                         { mid[j]=*itc; break; }
1962                     }
1963                   stNode=endNode;
1964                   endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
1965                 }
1966             }
1967         }
1968     }
1969   MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
1970   if(nbOf2DCellsToBeSplit==0)
1971     return ret.retn();
1972   //
1973   int *retPtr(ret->getPointer());
1974   for(int i=0;i<nCell;i++)
1975     if(cells2DToTreat[i])
1976       *retPtr++=i;
1977   //
1978   MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
1979   DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
1980   MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
1981   DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
1982   if(middleNeedsToBeUsed)
1983     { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
1984   MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
1985   int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
1986   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.
1987   setPartOfMySelf(ret->begin(),ret->end(),*modif);
1988   {
1989     bool areNodesMerged; int newNbOfNodes;
1990     if(nbOfNodesCreated!=0)
1991       MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
1992   }
1993   return ret.retn();
1994 }
1995
1996 /*!
1997  * 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.
1998  * 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).
1999  * 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
2000  * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2001  * This method works on any 2D geometric types of cell (even static one). If a cell is touched its type becomes dynamic automatically. For 2D "repaired" quadratic cells
2002  * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2003  *
2004  * 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
2005  * using new instance, idem for coordinates.
2006  *
2007  * 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.
2008  *
2009  * \return DataArrayInt  * - The list of cellIds in \a this that have at least one edge colinearized.
2010  *
2011  * \throw If \a this is not coherent.
2012  * \throw If \a this has not spaceDim equal to 2.
2013  * \throw If \a this has not meshDim equal to 2.
2014  *
2015  * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2016  */
2017 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2018 {
2019   internalColinearize2D(eps, false);
2020 }
2021
2022 /*!
2023  * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2024  * In a given 2D cell, if two edges are colinear and the junction point is used by a third edge, the two edges will not be
2025  * merged, contrary to colinearize2D().
2026  *
2027  * \sa MEDCouplingUMesh::colinearize2D
2028  */
2029 DataArrayInt *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2030 {
2031   internalColinearize2D(eps, true);
2032 }
2033
2034
2035 /*!
2036  * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2037  */
2038 DataArrayInt *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2039 {
2040   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2041   checkConsistencyLight();
2042   if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2043     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2044   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2045   int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2046   const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2047   MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2048   std::map<int, bool> forbiddenPoints;  // list of points that can not be removed (or it will break conformity)
2049   if(stayConform) //
2050     {
2051       // A point that is used by more than 2 edges can not be removed without breaking conformity:
2052       MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descI1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescI1(DataArrayInt::New());
2053       MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2054       MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2055       MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2056       MCAuto<DataArrayInt> dsi(revDescI2->deltaShiftIndex());
2057       MCAuto<DataArrayInt> ids(dsi->findIdsGreaterThan(2));
2058       const int * cPtr(mDesc0D->getNodalConnectivity()->begin());
2059       for(auto it = ids->begin(); it != ids->end(); it++)
2060          forbiddenPoints[cPtr[2*(*it)+1]] = true;  // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2061     }
2062
2063   MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2064   const double *coords(_coords->begin());
2065   int *newciptr(newci->getPointer());
2066   for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2067     {
2068       if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2069         ret->pushBackSilent(i);
2070       newciptr[1]=newc->getNumberOfTuples();
2071     }
2072   //
2073   if(ret->empty())
2074     return ret.retn();
2075   if(!appendedCoords->empty())
2076     {
2077       appendedCoords->rearrange(2);
2078       MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2079       //non const part
2080       setCoords(newCoords);
2081     }
2082   //non const part
2083   setConnectivity(newc,newci,true);
2084   return ret.retn();
2085 }
2086
2087 ///@cond INTERNAL
2088 /**
2089  * c, cI describe a wire mesh in 3D space, in local numbering
2090  * startNode, endNode in global numbering
2091  *\return true if the segment is indeed split
2092  */
2093 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2094                                             const int * c, const int * cI, const int *idsBg, const int *endBg,
2095                                             std::vector<int> & pointIds, std::vector<int> & hitSegs)
2096 {
2097   using namespace std;
2098
2099   const int SPACEDIM=3;
2100   typedef pair<double, int> PairDI;
2101   set< PairDI > x;
2102   for (const int * it = idsBg; it != endBg; ++it)
2103     {
2104       assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2105       int start = c[cI[*it]+1], end = c[cI[*it]+2];
2106       x.insert(make_pair(coo[start*SPACEDIM], start));  // take only X coords
2107       x.insert(make_pair(coo[end*SPACEDIM], end));
2108     }
2109
2110   vector<PairDI> xx(x.begin(), x.end());
2111   sort(xx.begin(),xx.end());
2112   pointIds.reserve(xx.size());
2113
2114   // Keep what is inside [startNode, endNode]:
2115   int go = 0;
2116   for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2117     {
2118       const int idx = (*it).second;
2119       if (!go)
2120         {
2121           if (idx == startNode)   go = 1;
2122           if (idx == endNode)     go = 2;
2123           if (go)                 pointIds.push_back(idx);
2124           continue;
2125         }
2126       pointIds.push_back(idx);
2127       if (idx == endNode || idx == startNode)
2128         break;
2129     }
2130
2131 //  vector<int> pointIds2(pointIds.size()+2);
2132 //  copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2133 //  pointIds2[0] = startNode;
2134 //  pointIds2[pointIds2.size()-1] = endNode;
2135
2136   if (go == 2)
2137     reverse(pointIds.begin(), pointIds.end());
2138
2139   // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2140   for (const int * it = idsBg; it != endBg; ++it)
2141     {
2142       int start = c[cI[*it]+1], end = c[cI[*it]+2];
2143       vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2144       if (itStart == pointIds.end()) continue;
2145       vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2146       if (itEnd == pointIds.end())               continue;
2147       if (abs(distance(itEnd, itStart)) != 1)    continue;
2148       hitSegs.push_back(*it);   // segment is undivided.
2149     }
2150
2151   return (pointIds.size() > 2); // something else apart start and end node
2152 }
2153
2154 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2155                                           const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2156 {
2157   using namespace std;
2158   int dst = distance(sIdxConn, sIdxConnE);
2159   modifiedFace.reserve(dst + insidePoints.size()-2);
2160   modifiedFace.resize(dst);
2161   copy(sIdxConn, sIdxConnE, modifiedFace.data());
2162
2163   vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2164   vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2165   if (startPos == shortEnd)
2166     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2167   vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2168   if (endPos == shortEnd)
2169     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2170   int d = distance(startPos, endPos);
2171   if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2172     modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end());  // insidePoints also contains start and end node. Those don't need to be inserted.
2173   else
2174     modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2175 }
2176
2177 ///@endcond
2178
2179
2180 /*!
2181  * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2182  * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2183  * This method performs a conformization of \b this.
2184  *
2185  * Only polyhedron cells are supported. You can call convertAllToPoly()
2186  *
2187  * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2188  * This method expects that all nodes in \a this are not closer than \a eps.
2189  * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2190  *
2191  * \param [in] eps the relative error to detect merged edges.
2192  * \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
2193  *                           that the user is expected to deal with.
2194  *
2195  * \throw If \a this is not coherent.
2196  * \throw If \a this has not spaceDim equal to 3.
2197  * \throw If \a this has not meshDim equal to 3.
2198  * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2199  */
2200 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2201 {
2202   using namespace std;
2203
2204   static const int SPACEDIM=3;
2205   checkConsistencyLight();
2206   if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2207     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2208   if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2209     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2210
2211   MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2212   const double * coo(_coords->begin());
2213   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2214
2215   {
2216     /*************************
2217      *  STEP 1  -- faces
2218      *************************/
2219     MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2220     MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2221     const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2222     const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2223     MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2224
2225     // Build BBTree
2226     MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2227     const double *bbox(bboxArr->begin()); getCoords()->begin();
2228     int nDescCell(mDesc->getNumberOfCells());
2229     BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2230     // Surfaces - handle biggest first
2231     MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2232     DataArrayDouble * surfs = surfF->getArray();
2233     // Normal field
2234     MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2235     DataArrayDouble * normals = normalsF->getArray();
2236     const double * normalsP = normals->getConstPointer();
2237
2238     // Sort faces by decreasing surface:
2239     vector< pair<double,int> > S;
2240     for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2241       {
2242         pair<double,int> p = make_pair(surfs->begin()[i], i);
2243         S.push_back(p);
2244       }
2245     sort(S.rbegin(),S.rend()); // reverse sort
2246     vector<bool> hit(nDescCell);
2247     fill(hit.begin(), hit.end(), false);
2248     vector<int> hitPoly; // the final result: which 3D cells have been modified.
2249
2250     for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2251       {
2252         int faceIdx = (*it).second;
2253         if (hit[faceIdx]) continue;
2254
2255         vector<int> candidates, cands2;
2256         myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2257         // Keep only candidates whose normal matches the normal of current face
2258         for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2259           {
2260             bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2261             if (*it2 != faceIdx && col)
2262               cands2.push_back(*it2);
2263           }
2264         if (!cands2.size())
2265           continue;
2266
2267         // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2268         INTERP_KERNEL::TranslationRotationMatrix rotation;
2269         INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2270                                                                    coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2271                                                                    coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2272
2273         MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false));  // false=zipCoords is called
2274         MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2275         double * cooPartRef(mPartRef->_coords->getPointer());
2276         double * cooPartCand(mPartCand->_coords->getPointer());
2277         for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2278           rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2279         for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2280           rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2281
2282         // Localize faces in 2D thanks to barycenters
2283         MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2284         vector<int> compo; compo.push_back(2);
2285         MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2286         MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2287         if (!idsGoodPlane->getNumberOfTuples())
2288           continue;
2289
2290         baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2291
2292         compo[0] = 0; compo.push_back(1);
2293         MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2294         mPartRef->changeSpaceDimension(2,0.0);
2295         MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2296         mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2297
2298         if (!cc->getNumberOfTuples())
2299           continue;
2300         MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2301
2302         {
2303           MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2304           if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2305             {
2306               ostringstream oss;
2307               oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2308               throw INTERP_KERNEL::Exception(oss.str());
2309             }
2310         }
2311
2312         MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2313         // Boundary face:
2314         if (!ids->getNumberOfTuples())
2315           continue;
2316
2317         double checkSurf=0.0;
2318         const int * idsGoodPlaneP(idsGoodPlane->begin());
2319         for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2320           {
2321             int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2322             hit[faceIdx2] = true;
2323             checkSurf += surfs->begin()[faceIdx2];
2324           }
2325         if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2326           {
2327             ostringstream oss;
2328             oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2329             throw INTERP_KERNEL::Exception(oss.str());
2330           }
2331
2332         // For all polyhedrons using this face, replace connectivity:
2333         vector<int> polyIndices, packsIds, facePack;
2334         for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2335           polyIndices.push_back(revDescP[ii]);
2336         ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2337
2338         // Current face connectivity
2339         const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2340         const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2341         connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2342         // Deletion of old faces
2343         int jj=0;
2344         for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2345           {
2346             if (packsIds[jj] == -1)
2347               // The below should never happen - if a face is used several times, with a different layout of the nodes
2348               // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2349               // faces which are actually used only once, by a single cell. This is different for edges below.
2350               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2351             else
2352               connSla->deletePack(*it2, packsIds[jj]);
2353           }
2354         // Insertion of new faces:
2355         for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2356           {
2357             // Build pack from the face to insert:
2358             int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2359             int facePack2Sz;
2360             const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2361             // Insert it in all hit polyhedrons:
2362             for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2363               connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz);  // without the type
2364           }
2365       }
2366   }  // end step1
2367
2368   // Set back modified connectivity
2369   MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2370   MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2371   connSla->convertToPolyhedronConn(cAuto, cIAuto);
2372
2373   {
2374     /************************
2375      *  STEP 2 -- edges
2376      ************************/
2377     // Now we have a face-conform mesh.
2378
2379     // Recompute descending
2380     MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2381     // Rebuild desc connectivity with orientation this time!!
2382     MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2383     const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2384     const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2385     const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2386     MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2387     MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2388     MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2389     MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2390     MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2391 //    std::cout << "writing!\n";
2392 //    mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2393 //    mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2394     const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2395     const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2396     MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2397     const double *bbox2(bboxArr->begin());
2398     int nDesc2Cell=mDesc2->getNumberOfCells();
2399     BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2400
2401     // Edges - handle longest first
2402     MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2403     DataArrayDouble * lens = lenF->getArray();
2404
2405     // Sort edges by decreasing length:
2406     vector<pair<double,int> > S;
2407     for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2408       {
2409         pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2410         S.push_back(p);
2411       }
2412     sort(S.rbegin(),S.rend()); // reverse sort
2413
2414     vector<bool> hit(nDesc2Cell);
2415     fill(hit.begin(), hit.end(), false);
2416
2417     for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2418       {
2419         int eIdx = (*it).second;
2420         if (hit[eIdx])
2421           continue;
2422
2423         vector<int> candidates, cands2;
2424         myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2425         // Keep only candidates colinear with current edge
2426         double vCurr[3];
2427         unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2428         for (int i3=0; i3 < 3; i3++)  // TODO: use fillSonCellNodalConnectivity2 or similar?
2429           vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2430         for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2431           {
2432             double vOther[3];
2433             unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2434             for (int i3=0; i3 < 3; i3++)
2435               vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2436             bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2437             // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2438             // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2439             if (col)
2440               cands2.push_back(*it2);
2441           }
2442         if (cands2.size() == 1 && cands2[0] == eIdx)  // see warning above
2443           continue;
2444
2445         // Now rotate edges to bring them on Ox
2446         int startNode = cDesc2[cIDesc2[eIdx]+1];
2447         int endNode = cDesc2[cIDesc2[eIdx]+2];
2448         INTERP_KERNEL::TranslationRotationMatrix rotation;
2449         INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2450         MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false));  // false=zipCoords is called
2451         MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2452         MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2453         int nbElemsNotM1;
2454         {
2455           MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2456           nbElemsNotM1 = tmp->getNbOfElems();
2457         }
2458         MCAuto<DataArrayInt>  nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2459         double * cooPartRef(mPartRef->_coords->getPointer());
2460         double * cooPartCand(mPartCand->_coords->getPointer());
2461         for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2462           rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2463         for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2464           rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2465
2466
2467         // Eliminate all edges for which y or z is not null
2468         MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2469         vector<int> compo; compo.push_back(1);
2470         MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2471         compo[0] = 2;
2472         MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2473         MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2474         MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2475         MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2476         if (!idsGoodLine->getNumberOfTuples())
2477           continue;
2478
2479         // Now the ordering along the Ox axis:
2480         std::vector<int> insidePoints, hitSegs;
2481         bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2482             mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2483             idsGoodLine->begin(), idsGoodLine->end(),
2484             /*out*/insidePoints, hitSegs);
2485         // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2486         for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2487           hit[cands2[*its]] = true;
2488
2489         if (!isSplit)  // current segment remains in one piece
2490           continue;
2491
2492         // Get original node IDs in global coords array
2493         for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2494           *iit = nodeMapInv->begin()[*iit];
2495
2496         vector<int> polyIndices, packsIds, facePack;
2497         // For each face implying this edge
2498         for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2499           {
2500             int faceIdx = revDescP2[ii];
2501             // each cell where this face is involved connectivity will be modified:
2502             ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2503
2504             // Current face connectivity
2505             const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2506             const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2507
2508             vector<int> modifiedFace;
2509             ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2510             modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2511             connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2512           }
2513       }
2514
2515     // Rebuild 3D connectivity from descending:
2516     MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2517     MCAuto<DataArrayInt> superIdx(DataArrayInt::New());  superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2518     MCAuto<DataArrayInt> idx(DataArrayInt::New());       idx->alloc(1);                         idx->fillWithValue(0);
2519     MCAuto<DataArrayInt> vals(DataArrayInt::New());      vals->alloc(0);
2520     newConn->set3(superIdx, idx, vals);
2521     for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2522       for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2523         {
2524           int sz, faceIdx = abs(descP[jj])-1;
2525           bool orient = descP[jj]>0;
2526           const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2527           if (orient)
2528             newConn->pushBackPack(ii, p+1, p+sz);  // +1 to skip type
2529           else
2530             {
2531               vector<int> rev(sz-1);
2532               for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2533               newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2534             }
2535         }
2536     // And finally:
2537     newConn->convertToPolyhedronConn(cAuto, cIAuto);
2538   } // end step2
2539
2540   ret = ret->buildUniqueNotSorted();
2541   return ret.retn();
2542 }
2543
2544