src/MEDCoupling/MEDCouplingUMesh_intersection.cxx

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