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;
 867   int _right;
 868 };
 869
 870 void EdgeInfo::somethingHappendAt(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
 871 {
 872   const MEDCouplingUMesh *mesh(_mesh);
 873   if(mesh)
 874     return ;
 875   if(_right<pos)
 876     return ;
 877   if(_left>pos)
 878     { _left++; _right++; return ; }
 879   if (_right > pos && _left != pos)
 880     { _right++; return ; }
 881   if(_right==pos)
 882     {
 883       bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
 884       if((isLeft && isRight) || (!isLeft && !isRight))
 885         throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 1 !");
 886       if(isLeft)
 887         return ;
 888       if(isRight)
 889         {
 890           _right++;
 891           return ;
 892         }
 893     }
 894   if(_left==pos)
 895     {
 896       bool isLeft(std::find(newLeft.begin(),newLeft.end(),_edge)!=newLeft.end()),isRight(std::find(newRight.begin(),newRight.end(),_edge)!=newRight.end());
 897       if((isLeft && isRight) || (!isLeft && !isRight))
 898         throw INTERP_KERNEL::Exception("EdgeInfo::somethingHappendAt : internal error # 2 !");
 899       if(isLeft)
 900         {
 901           _right++;
 902           return ;
 903         }
 904       if(isRight)
 905         {
 906           _left++;
 907           _right++;
 908           return ;
 909         }
 910     }
 911 }
 912
 913 void EdgeInfo::feedEdgeInfoAt(double eps, const MEDCouplingUMesh *mesh2D, int offset, int neighbors[2]) const
 914 {
 915   const MEDCouplingUMesh *mesh(_mesh);
 916   if(!mesh)
 917     {
 918       neighbors[0]=offset+_left; neighbors[1]=offset+_right;
 919     }
 920   else
 921     {// not fully splitting cell case
 922       if(mesh2D->getNumberOfCells()==1)
 923         {//little optimization. 1 cell no need to find in which cell mesh is !
 924           neighbors[0]=offset; neighbors[1]=offset;
 925           return;
 926         }
 927       else
 928         {
 929           MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
 930           int cellId(mesh2D->getCellContainingPoint(barys->begin(),eps));
 931           if(cellId==-1)
 932             throw INTERP_KERNEL::Exception("EdgeInfo::feedEdgeInfoAt : internal error !");
 933           neighbors[0]=offset+cellId; neighbors[1]=offset+cellId;
 934         }
 935     }
 936 }
 937
 938 class VectorOfCellInfo
 939 {
 940 public:
 941   VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr);
 942   std::size_t size() const { return _pool.size(); }
 943   int getPositionOf(double eps, const MEDCouplingUMesh *mesh) const;
 944   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);
 945   const std::vector<int>& getConnOf(int pos) const { return get(pos)._edges; }
 946   const std::vector< MCAuto<INTERP_KERNEL::Edge> >& getEdgePtrOf(int pos) const { return get(pos)._edges_ptr; }
 947   MCAuto<MEDCouplingUMesh> getZeMesh() const { return _ze_mesh; }
 948   void feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const;
 949 private:
 950   int getZePosOfEdgeGivenItsGlobalId(int pos) const;
 951   void updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight);
 952   const CellInfo& get(int pos) const;
 953   CellInfo& get(int pos);
 954 private:
 955   std::vector<CellInfo> _pool;
 956   MCAuto<MEDCouplingUMesh> _ze_mesh;
 957   std::vector<EdgeInfo> _edge_info;
 958 };
 959
 960 VectorOfCellInfo::VectorOfCellInfo(const std::vector<int>& edges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& edgesPtr):_pool(1)
 961 {
 962   _pool[0]._edges=edges;
 963   _pool[0]._edges_ptr=edgesPtr;
 964 }
 965
 966 int VectorOfCellInfo::getPositionOf(double eps, const MEDCouplingUMesh *mesh) const
 967 {
 968   if(_pool.empty())
 969     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : empty !");
 970   if(_pool.size()==1)
 971     return 0;
 972   const MEDCouplingUMesh *zeMesh(_ze_mesh);
 973   if(!zeMesh)
 974     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::getPositionOf : null aggregated mesh !");
 975   MCAuto<DataArrayDouble> barys(mesh->computeCellCenterOfMass());
 976   return zeMesh->getCellContainingPoint(barys->begin(),eps);
 977 }
 978
 979 void VectorOfCellInfo::setMeshAt(std::size_t pos, const MCAuto<MEDCouplingUMesh>& mesh, int istart, int iend,
 980                                  const MCAuto<MEDCouplingUMesh>& mesh1DInCase, const std::vector< std::vector<int> >& edges,
 981                                  const std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > >& edgePtrs)
 982 {
 983   get(pos);//to check pos
 984   bool isFast(pos==0 && _pool.size()==1);
 985   std::size_t sz(edges.size());
 986   // dealing with edges
 987   if(sz==1)
 988     _edge_info.push_back(EdgeInfo(istart,iend,mesh1DInCase));
 989   else
 990     _edge_info.push_back(EdgeInfo(istart,iend,pos,edgePtrs[0].back()));
 991   //
 992   std::vector<CellInfo> pool(_pool.size()-1+sz);
 993   for(std::size_t i=0;i<pos;i++)
 994     pool[i]=_pool[i];
 995   for(std::size_t j=0;j<sz;j++)
 996     pool[pos+j]=CellInfo(edges[j],edgePtrs[j]);
 997   for(int i=pos+1;i<(int)_pool.size();i++)
 998     pool[i+sz-1]=_pool[i];
 999   _pool=pool;
1000   //
1001   if(sz==2)
1002     updateEdgeInfo(pos,edgePtrs[0],edgePtrs[1]);
1003   //
1004   if(isFast)
1005     {
1006       _ze_mesh=mesh;
1007       return ;
1008     }
1009   //
1010   std::vector< MCAuto<MEDCouplingUMesh> > ms;
1011   if(pos>0)
1012     {
1013       MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(0,pos,true)));
1014       ms.push_back(elt);
1015     }
1016   ms.push_back(mesh);
1017   if(pos<_ze_mesh->getNumberOfCells()-1)
1018   {
1019     MCAuto<MEDCouplingUMesh> elt(static_cast<MEDCouplingUMesh *>(_ze_mesh->buildPartOfMySelfSlice(pos+1,_ze_mesh->getNumberOfCells(),true)));
1020     ms.push_back(elt);
1021   }
1022   std::vector< const MEDCouplingUMesh *> ms2(ms.size());
1023   for(std::size_t j=0;j<ms2.size();j++)
1024     ms2[j]=ms[j];
1025   _ze_mesh=MEDCouplingUMesh::MergeUMeshesOnSameCoords(ms2);
1026 }
1027
1028 void VectorOfCellInfo::feedEdgeInfoAt(double eps, int pos, int offset, int neighbors[2]) const
1029 {
1030   _edge_info[getZePosOfEdgeGivenItsGlobalId(pos)].feedEdgeInfoAt(eps,_ze_mesh,offset,neighbors);
1031 }
1032
1033 int VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId(int pos) const
1034 {
1035   if(pos<0)
1036     throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id ! Must be >=0 !");
1037   int ret(0);
1038   for(std::vector<EdgeInfo>::const_iterator it=_edge_info.begin();it!=_edge_info.end();it++,ret++)
1039     {
1040       if((*it).isInMyRange(pos))
1041         return ret;
1042     }
1043   throw INTERP_KERNEL::Exception("VectorOfCellInfo::getZePosOfEdgeGivenItsGlobalId : invalid id !");
1044 }
1045
1046 void VectorOfCellInfo::updateEdgeInfo(int pos, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newLeft, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& newRight)
1047 {
1048   get(pos);//to perform the sanity check;
1049   if(_edge_info.empty())
1050     return ;
1051   std::size_t sz(_edge_info.size()-1);
1052   for(std::size_t i=0;i<sz;i++)
1053     _edge_info[i].somethingHappendAt(pos,newLeft,newRight);
1054 }
1055
1056 const CellInfo& VectorOfCellInfo::get(int pos) const
1057 {
1058   if(pos<0 || pos>=(int)_pool.size())
1059     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get const : invalid pos !");
1060   return _pool[pos];
1061 }
1062
1063 CellInfo& VectorOfCellInfo::get(int pos)
1064 {
1065   if(pos<0 || pos>=(int)_pool.size())
1066     throw INTERP_KERNEL::Exception("VectorOfCellSplitter::get : invalid pos !");
1067   return _pool[pos];
1068 }
1069
1070 /*!
1071  * Given :
1072  * - a \b closed set of edges ( \a allEdges and \a allEdgesPtr ) that defines the split descending 2D cell.
1073  * - \a splitMesh1D a split 2D curve mesh contained into 2D cell defined above.
1074  *
1075  * This method returns the 2D mesh and feeds \a idsLeftRight using offset.
1076  *
1077  * Algorithm : \a splitMesh1D is cut into contiguous parts. Each contiguous parts will build incrementally the output 2D cells.
1078  *
1079  * \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.
1080  */
1081 MEDCouplingUMesh *BuildMesh2DCutInternal(double eps, MEDCouplingUMesh *splitMesh1D, const std::vector<int>& allEdges, const std::vector< MCAuto<INTERP_KERNEL::Edge> >& allEdgesPtr, int offset,
1082                                          MCAuto<DataArrayInt>& idsLeftRight)
1083 {
1084   int nbCellsInSplitMesh1D(splitMesh1D->getNumberOfCells());
1085   if(nbCellsInSplitMesh1D==0)
1086     throw INTERP_KERNEL::Exception("BuildMesh2DCutInternal : internal error ! input 1D mesh must have at least one cell !");
1087   const int *cSplitPtr(splitMesh1D->getNodalConnectivity()->begin()),*ciSplitPtr(splitMesh1D->getNodalConnectivityIndex()->begin());
1088   std::size_t nb(allEdges.size()),jj;
1089   if(nb%2!=0)
1090     throw INTERP_KERNEL::Exception("BuildMesh2DCutFrom : internal error 2 !");
1091   std::vector<int> edge1Bis(nb*2);
1092   std::vector< MCAuto<INTERP_KERNEL::Edge> > edge1BisPtr(nb*2);
1093   std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin());
1094   std::copy(allEdges.begin(),allEdges.end(),edge1Bis.begin()+nb);
1095   std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin());
1096   std::copy(allEdgesPtr.begin(),allEdgesPtr.end(),edge1BisPtr.begin()+nb);
1097   //
1098   idsLeftRight=DataArrayInt::New(); idsLeftRight->alloc(nbCellsInSplitMesh1D*2); idsLeftRight->fillWithValue(-2); idsLeftRight->rearrange(2);
1099   int *idsLeftRightPtr(idsLeftRight->getPointer());
1100   VectorOfCellInfo pool(edge1Bis,edge1BisPtr);
1101
1102   // Compute contiguous parts of splitMesh1D. We can not make the full assumption that segments are consecutive in the connectivity
1103   // (even if the user correctly called orderConsecutiveCells1D()). Indeed the tool might be a closed line whose junction point is in
1104   // splitMesh1D. There can be only one such a point, and if this happens this is necessarily at the start
1105   // of the connectivity.
1106   MCAuto <DataArrayInt> renumb(DataArrayInt::New());
1107   renumb->alloc(nbCellsInSplitMesh1D,1);
1108   const int * renumbP(renumb->begin());
1109
1110   int i, first=cSplitPtr[1];
1111   // Follow 1D line backward as long as it is connected:
1112   for (i=nbCellsInSplitMesh1D-1; cSplitPtr[ciSplitPtr[i]+2] == first; i--)
1113     first=cSplitPtr[ciSplitPtr[i]+1];
1114   if (i < nbCellsInSplitMesh1D-1)
1115     {
1116       // Build circular permutation to shift consecutive edges together
1117       renumb->iota(i+1);
1118       renumb->applyModulus(nbCellsInSplitMesh1D);
1119       splitMesh1D->renumberCells(renumbP, false);
1120       cSplitPtr = splitMesh1D->getNodalConnectivity()->begin();
1121       ciSplitPtr = splitMesh1D->getNodalConnectivityIndex()->begin();
1122     }
1123   else
1124     renumb->iota();
1125   //
1126
1127   for(int iStart=0;iStart<nbCellsInSplitMesh1D;)
1128     {// split [0:nbCellsInSplitMesh1D) in contiguous parts [iStart:iEnd)
1129       int iEnd(iStart);
1130       for(;iEnd<nbCellsInSplitMesh1D;)
1131         {
1132           for(jj=0;jj<nb && edge1Bis[2*jj+1]!=cSplitPtr[ciSplitPtr[iEnd]+2];jj++);
1133           if(jj!=nb)
1134             break;
1135           else
1136             iEnd++;
1137         }
1138       if(iEnd<nbCellsInSplitMesh1D)
1139         iEnd++;
1140
1141       MCAuto<MEDCouplingUMesh> partOfSplitMesh1D(static_cast<MEDCouplingUMesh *>(splitMesh1D->buildPartOfMySelfSlice(iStart,iEnd,1,true)));
1142       int pos(pool.getPositionOf(eps,partOfSplitMesh1D));
1143       //
1144       MCAuto<MEDCouplingUMesh>retTmp(MEDCouplingUMesh::New("",2));
1145       retTmp->setCoords(splitMesh1D->getCoords());
1146       retTmp->allocateCells();
1147
1148       std::vector< std::vector<int> > out0;
1149       std::vector< std::vector< MCAuto<INTERP_KERNEL::Edge> > > out1;
1150
1151       BuildMesh2DCutInternal2(partOfSplitMesh1D,pool.getConnOf(pos),pool.getEdgePtrOf(pos),out0,out1);
1152       for(std::size_t cnt=0;cnt<out0.size();cnt++)
1153         AddCellInMesh2D(retTmp,out0[cnt],out1[cnt]);
1154       pool.setMeshAt(pos,retTmp,iStart,iEnd,partOfSplitMesh1D,out0,out1);
1155       //
1156       iStart=iEnd;
1157     }
1158   for(int mm=0;mm<nbCellsInSplitMesh1D;mm++)
1159     pool.feedEdgeInfoAt(eps,renumbP[mm],offset,idsLeftRightPtr+2*mm);
1160
1161   return pool.getZeMesh().retn();
1162 }
1163
1164 /*
1165  * splitMesh1D is an input parameter but might have its cells renumbered.
1166  */
1167 MEDCouplingUMesh *BuildMesh2DCutFrom(double eps, int cellIdInMesh2D, const MEDCouplingUMesh *mesh2DDesc, MEDCouplingUMesh *splitMesh1D,
1168                                      const int *descBg, const int *descEnd, const std::vector< std::vector<int> >& intersectEdge1, int offset,
1169                                      MCAuto<DataArrayInt>& idsLeftRight)
1170 {
1171   const int *cdescPtr(mesh2DDesc->getNodalConnectivity()->begin()),*cidescPtr(mesh2DDesc->getNodalConnectivityIndex()->begin());
1172   //
1173   std::vector<int> allEdges;
1174   std::vector< MCAuto<INTERP_KERNEL::Edge> > allEdgesPtr; // for each sub edge in splitMesh2D the uncut Edge object of the original mesh2D
1175   for(const int *it(descBg);it!=descEnd;it++) // for all edges in the descending connectivity of the 2D mesh in relative Fortran mode
1176     {
1177       int edgeId(std::abs(*it)-1);
1178       std::map< MCAuto<INTERP_KERNEL::Node>,int> m;
1179       MCAuto<INTERP_KERNEL::Edge> ee(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)cdescPtr[cidescPtr[edgeId]],cdescPtr+cidescPtr[edgeId]+1,mesh2DDesc->getCoords()->begin(),m));
1180       const std::vector<int>& edge1(intersectEdge1[edgeId]);
1181       if(*it>0)
1182         allEdges.insert(allEdges.end(),edge1.begin(),edge1.end());
1183       else
1184         allEdges.insert(allEdges.end(),edge1.rbegin(),edge1.rend());
1185       std::size_t sz(edge1.size());
1186       for(std::size_t cnt=0;cnt<sz;cnt++)
1187         allEdgesPtr.push_back(ee);
1188     }
1189   //
1190   return BuildMesh2DCutInternal(eps,splitMesh1D,allEdges,allEdgesPtr,offset,idsLeftRight);
1191 }
1192
1193 bool AreEdgeEqual(const double *coo2D, const INTERP_KERNEL::CellModel& typ1, const int *conn1, const INTERP_KERNEL::CellModel& typ2, const int *conn2, double eps)
1194 {
1195   if(!typ1.isQuadratic() && !typ2.isQuadratic())
1196     {//easy case comparison not
1197       return conn1[0]==conn2[0] && conn1[1]==conn2[1];
1198     }
1199   else if(typ1.isQuadratic() && typ2.isQuadratic())
1200     {
1201       bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1202       if(!status0)
1203         return false;
1204       if(conn1[2]==conn2[2])
1205         return true;
1206       const double *a(coo2D+2*conn1[2]),*b(coo2D+2*conn2[2]);
1207       double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1208       return dist<eps;
1209     }
1210   else
1211     {//only one is quadratic
1212       bool status0(conn1[0]==conn2[0] && conn1[1]==conn2[1]);
1213       if(!status0)
1214         return false;
1215       const double *a(0),*bb(0),*be(0);
1216       if(typ1.isQuadratic())
1217         {
1218           a=coo2D+2*conn1[2]; bb=coo2D+2*conn2[0]; be=coo2D+2*conn2[1];
1219         }
1220       else
1221         {
1222           a=coo2D+2*conn2[2]; bb=coo2D+2*conn1[0]; be=coo2D+2*conn1[1];
1223         }
1224       double b[2]; b[0]=(be[0]+bb[0])/2.; b[1]=(be[1]+bb[1])/2.;
1225       double dist(sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])));
1226       return dist<eps;
1227     }
1228 }
1229
1230 /*!
1231  * This method returns among the cellIds [ \a candidatesIn2DBg , \a candidatesIn2DEnd ) in \a mesh2DSplit those exactly sharing \a cellIdInMesh1DSplitRelative in \a mesh1DSplit.
1232  * \a mesh2DSplit and \a mesh1DSplit are expected to share the coordinates array.
1233  *
1234  * \param [in] cellIdInMesh1DSplitRelative is in Fortran mode using sign to specify direction.
1235  */
1236 int FindRightCandidateAmong(const MEDCouplingUMesh *mesh2DSplit, const int *candidatesIn2DBg, const int *candidatesIn2DEnd, const MEDCouplingUMesh *mesh1DSplit, int cellIdInMesh1DSplitRelative, double eps)
1237 {
1238   if(candidatesIn2DEnd==candidatesIn2DBg)
1239     throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 1 !");
1240   const double *coo(mesh2DSplit->getCoords()->begin());
1241   if(std::distance(candidatesIn2DBg,candidatesIn2DEnd)==1)
1242     return *candidatesIn2DBg;
1243   int edgeId(std::abs(cellIdInMesh1DSplitRelative)-1);
1244   MCAuto<MEDCouplingUMesh> cur1D(static_cast<MEDCouplingUMesh *>(mesh1DSplit->buildPartOfMySelf(&edgeId,&edgeId+1,true)));
1245   if(cellIdInMesh1DSplitRelative<0)
1246     cur1D->changeOrientationOfCells();
1247   const int *c1D(cur1D->getNodalConnectivity()->begin());
1248   const INTERP_KERNEL::CellModel& ref1DType(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c1D[0]));
1249   for(const int *it=candidatesIn2DBg;it!=candidatesIn2DEnd;it++)
1250     {
1251       MCAuto<MEDCouplingUMesh> cur2D(static_cast<MEDCouplingUMesh *>(mesh2DSplit->buildPartOfMySelf(it,it+1,true)));
1252       const int *c(cur2D->getNodalConnectivity()->begin()),*ci(cur2D->getNodalConnectivityIndex()->begin());
1253       const INTERP_KERNEL::CellModel &cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[0]]));
1254       unsigned sz(cm.getNumberOfSons2(c+ci[0]+1,ci[1]-ci[0]-1));
1255       INTERP_KERNEL::AutoPtr<int> tmpPtr(new int[ci[1]-ci[0]]);
1256       for(unsigned it2=0;it2<sz;it2++)
1257         {
1258           INTERP_KERNEL::NormalizedCellType typeOfSon;
1259           cm.fillSonCellNodalConnectivity2(it2,c+ci[0]+1,ci[1]-ci[0]-1,tmpPtr,typeOfSon);
1260           const INTERP_KERNEL::CellModel &curCM(INTERP_KERNEL::CellModel::GetCellModel(typeOfSon));
1261           if(AreEdgeEqual(coo,ref1DType,c1D+1,curCM,tmpPtr,eps))
1262             return *it;
1263         }
1264     }
1265   throw INTERP_KERNEL::Exception("FindRightCandidateAmong : internal error 2 ! Unable to find the edge among split cell !");
1266 }
1267
1268 /*!
1269  * \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.
1270  *                               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.
1271  *                               And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
1272  * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
1273  * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
1274  * \param [out] addCoo - nodes to be append at the end
1275  * \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.
1276  */
1277 void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
1278                                          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)
1279 {
1280   static const int SPACEDIM=2;
1281   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1282   const int *c1(m1Desc->getNodalConnectivity()->begin()),*ci1(m1Desc->getNodalConnectivityIndex()->begin());
1283   // Build BB tree of all edges in the tool mesh (second mesh)
1284   MCAuto<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree(eps)),bbox2Arr(m2Desc->getBoundingBoxForBBTree(eps));
1285   const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1286   int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
1287   intersectEdge1.resize(nDescCell1);
1288   colinear2.resize(nDescCell2);
1289   subDiv2.resize(nDescCell2);
1290   BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
1291   BBTreePts<SPACEDIM,int> treeNodes2(m2Desc->getCoords()->begin(),0,0,m2Desc->getCoords()->getNumberOfTuples(),eps);
1292
1293   std::vector<int> candidates1(1);
1294   int offset1(m1Desc->getNumberOfNodes());
1295   int offset2(offset1+m2Desc->getNumberOfNodes());
1296   for(int i=0;i<nDescCell1;i++)  // for all edges in the first mesh
1297     {
1298       std::vector<int> candidates2; // edges of mesh2 candidate for intersection
1299       myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1300       if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
1301         {
1302           std::map<INTERP_KERNEL::Node *,int> map1,map2;
1303           std::map<int, INTERP_KERNEL::Node *> revMap2;
1304           // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
1305           INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
1306           // Build revMap2
1307           for (auto& kv : map2)
1308             revMap2[kv.second] = kv.first;
1309           candidates1[0]=i;
1310           // In the construction of pol1 we might reuse nodes from pol2, that we have identified as to be merged.
1311           INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMeshWithTree(m1Desc,candidates1,map1,treeNodes2, revMap2);
1312           // 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
1313           // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
1314           std::set<INTERP_KERNEL::Node *> nodes;
1315           pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
1316           std::size_t szz(nodes.size());
1317           std::vector< MCAuto<INTERP_KERNEL::Node> > nodesSafe(szz);
1318           std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
1319           for(std::size_t iii=0;iii<szz;iii++,itt++)
1320             { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
1321           // end of protection
1322           // Performs edge cutting:
1323           pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo,mergedNodes);
1324           delete pol2;
1325           delete pol1;
1326         }
1327       else
1328         // Copy the edge (take only the two first points, ie discard quadratic point at this stage)
1329         intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i]+3);
1330     }
1331 }
1332
1333
1334 /*!
1335  * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
1336  * It builds the descending connectivity of the two meshes, and then using a binary tree
1337  * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
1338  * Documentation about parameters  colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
1339  */
1340 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
1341                                                    std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
1342                                                    MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
1343                                                    std::vector<double>& addCoo,
1344                                                    MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
1345 {
1346   // Build desc connectivity
1347   desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
1348   desc2=DataArrayInt::New();
1349   descIndx2=DataArrayInt::New();
1350   revDesc2=DataArrayInt::New();
1351   revDescIndx2=DataArrayInt::New();
1352   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1353   MCAuto<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
1354   m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
1355   m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
1356   MCAuto<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
1357   std::map<int,int> notUsedMap;
1358   Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo,notUsedMap);
1359   m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
1360   m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
1361 }
1362
1363 /**
1364  * Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
1365  * (newly created) nodes corresponding to the edge intersections.
1366  * Output params:
1367  * @param[out] cr, crI connectivity of the resulting mesh
1368  * @param[out] cNb1, cNb2 correspondence arrays giving for the merged mesh the initial cells IDs in m1 / m2
1369  * TODO: describe input parameters
1370  */
1371 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
1372                                                          const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
1373                                                          const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
1374                                                          const std::vector<double>& addCoords,
1375                                                          std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
1376 {
1377   static const int SPACEDIM=2;
1378   const double *coo1(m1->getCoords()->begin());
1379   const int *conn1(m1->getNodalConnectivity()->begin()),*connI1(m1->getNodalConnectivityIndex()->begin());
1380   int offset1(m1->getNumberOfNodes());
1381   const double *coo2(m2->getCoords()->begin());
1382   const int *conn2(m2->getNodalConnectivity()->begin()),*connI2(m2->getNodalConnectivityIndex()->begin());
1383   int offset2(offset1+m2->getNumberOfNodes());
1384   int offset3(offset2+((int)addCoords.size())/2);
1385   MCAuto<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree(eps)),bbox2Arr(m2->getBoundingBoxForBBTree(eps));
1386   const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
1387   // Here a BBTree on 2D-cells, not on segments:
1388   BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
1389   int ncell1(m1->getNumberOfCells());
1390   crI.push_back(0);
1391   for(int i=0;i<ncell1;i++)
1392     {
1393       std::vector<int> candidates2;
1394       myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
1395       std::map<INTERP_KERNEL::Node *,int> mapp;
1396       std::map<int,INTERP_KERNEL::Node *> mappRev;
1397       INTERP_KERNEL::QuadraticPolygon pol1;
1398       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
1399       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1400       // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
1401       MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
1402       // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
1403       pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
1404           desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
1405       //
1406       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
1407       std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
1408       INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
1409       for(it1.first();!it1.finished();it1.next())
1410         edges1.insert(it1.current()->getPtr());
1411       //
1412       std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
1413       std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
1414       int ii=0;
1415       // Build, for each intersecting cell candidate from mesh2, the corresponding QP.
1416       // Again all the additional intersecting nodes are there.
1417       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1418         {
1419           INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
1420           const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
1421           // Complete mapping with elements coming from the current cell it2 in mesh2:
1422           MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
1423           // pol2 is the new QP in the final merged result.
1424           pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
1425               pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
1426         }
1427       // The cleaning below must be done after the full construction of all pol2s to correctly deal with shared edges:
1428       for (auto &p: pol2s)
1429         p.cleanDegeneratedConsecutiveEdges();
1430       edgesIn2ForShare.clear();  // removing temptation to use it further since it might now contain invalid edges.
1431       ///
1432       ii=0;
1433       // Now rebuild intersected cells from all this:
1434       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
1435         {
1436           INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
1437           pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
1438           //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
1439           pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1440         }
1441       // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
1442       // by m2 but that we still want to keep in the final result.
1443       if(!edges1.empty())
1444         {
1445           try
1446           {
1447               INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
1448           }
1449           catch(INTERP_KERNEL::Exception& e)
1450           {
1451               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();
1452               throw INTERP_KERNEL::Exception(oss.str());
1453           }
1454         }
1455       for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
1456         (*it).second->decrRef();
1457     }
1458 }
1459
1460 void InsertNodeInConnIfNecessary(int nodeIdToInsert, std::vector<int>& conn, const double *coords, double eps)
1461 {
1462   std::vector<int>::iterator it(std::find(conn.begin(),conn.end(),nodeIdToInsert));
1463   if(it!=conn.end())
1464     return ;
1465   std::size_t sz(conn.size());
1466   std::size_t found(std::numeric_limits<std::size_t>::max());
1467   for(std::size_t i=0;i<sz;i++)
1468     {
1469       int pt0(conn[i]),pt1(conn[(i+1)%sz]);
1470       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]};
1471       double normm(sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2]));
1472       std::transform(v1,v1+3,v1,std::bind2nd(std::multiplies<double>(),1./normm));
1473       std::transform(v2,v2+3,v2,std::bind2nd(std::multiplies<double>(),1./normm));
1474       double v3[3];
1475       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];
1476       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]);
1477       if(normm2<eps)
1478         if(dotTest>eps && dotTest<1.-eps)
1479           {
1480             found=i;
1481             break;
1482           }
1483     }
1484   if(found==std::numeric_limits<std::size_t>::max())
1485     throw INTERP_KERNEL::Exception("InsertNodeInConnIfNecessary : not found point !");
1486   conn.insert(conn.begin()+(found+1)%sz,nodeIdToInsert);
1487 }
1488
1489 void SplitIntoToPart(const std::vector<int>& conn, int pt0, int pt1, std::vector<int>& part0, std::vector<int>& part1)
1490 {
1491   std::size_t sz(conn.size());
1492   std::vector<int> *curPart(&part0);
1493   for(std::size_t i=0;i<sz;i++)
1494     {
1495       int nextt(conn[(i+1)%sz]);
1496       (*curPart).push_back(nextt);
1497       if(nextt==pt0 || nextt==pt1)
1498         {
1499           if(curPart==&part0)
1500             curPart=&part1;
1501           else
1502             curPart=&part0;
1503           (*curPart).push_back(nextt);
1504         }
1505     }
1506 }
1507
1508 /*!
1509  * 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.
1510  */
1511 void MEDCouplingUMesh::buildSubCellsFromCut(const std::vector< std::pair<int,int> >& cut3DSurf,
1512                                             const int *desc, const int *descIndx, const double *coords, double eps,
1513                                             std::vector<std::vector<int> >& res) const
1514 {
1515   checkFullyDefined();
1516   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1517     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1518   const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1519   int nbOfCells(getNumberOfCells());
1520   if(nbOfCells!=1)
1521     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSubCellsFromCut works only with single cell presently !");
1522   for(int i=0;i<nbOfCells;i++)
1523     {
1524       int offset(descIndx[i]),nbOfFaces(descIndx[i+1]-offset);
1525       for(int j=0;j<nbOfFaces;j++)
1526         {
1527           const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1528           const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]));
1529           int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1530           INTERP_KERNEL::AutoPtr<int> tmp(new int[sz]);
1531           INTERP_KERNEL::NormalizedCellType cmsId;
1532           unsigned nbOfNodesSon(cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId));
1533           std::vector<int> elt((int *)tmp,(int *)tmp+nbOfNodesSon);
1534           if(p.first!=-1 && p.second!=-1)
1535             {
1536               if(p.first!=-2)
1537                 {
1538                   InsertNodeInConnIfNecessary(p.first,elt,coords,eps);
1539                   InsertNodeInConnIfNecessary(p.second,elt,coords,eps);
1540                   std::vector<int> elt1,elt2;
1541                   SplitIntoToPart(elt,p.first,p.second,elt1,elt2);
1542                   res.push_back(elt1);
1543                   res.push_back(elt2);
1544                 }
1545               else
1546                 res.push_back(elt);
1547             }
1548           else
1549             res.push_back(elt);
1550         }
1551     }
1552 }
1553
1554 /*!
1555  * 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).
1556  *
1557  * \sa MEDCouplingUMesh::split2DCells
1558  */
1559 void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
1560 {
1561   checkConnectivityFullyDefined();
1562   int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+subNodesInSeg->getNumberOfTuples());
1563   MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1564   const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1565   int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1566   int prevPosOfCi(ciPtr[0]);
1567   for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1568     {
1569       int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
1570       *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
1571       for(int j=0;j<sz;j++)
1572         {
1573           int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
1574           for(int k=0;k<sz2;k++)
1575             *cPtr++=subPtr[offset2+k];
1576           if(j!=sz-1)
1577             *cPtr++=oldConn[prevPosOfCi+j+2];
1578           deltaSz+=sz2;
1579         }
1580       prevPosOfCi=ciPtr[1];
1581       ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1582     }
1583   if(c->end()!=cPtr)
1584     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
1585   _nodal_connec->decrRef();
1586   _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
1587 }
1588
1589
1590 /*!
1591  * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additional nodes can be added at the end of coordinates array object.
1592  *
1593  * \return  int - the number of new nodes created.
1594  * \sa MEDCouplingUMesh::split2DCells
1595  */
1596 int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
1597 {
1598   checkConsistencyLight();
1599   int ncells(getNumberOfCells()),lgthToReach(getNodalConnectivityArrayLen()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
1600   MCAuto<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
1601   MCAuto<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
1602   const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
1603   const int *midPtr(mid->begin()),*midIPtr(midI->begin());
1604   const double *oldCoordsPtr(getCoords()->begin());
1605   int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
1606   int prevPosOfCi(ciPtr[0]);
1607   for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
1608     {
1609       int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
1610       for(int j=0;j<sz;j++)
1611         { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
1612       *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
1613       for(int j=0;j<sz;j++)//loop over subedges of oldConn
1614         {
1615           int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
1616           if(sz2==0)
1617             {
1618               if(j<sz-1)
1619                 cPtr[1]=oldConn[prevPosOfCi+2+j];
1620               cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
1621               continue;
1622             }
1623           std::vector<INTERP_KERNEL::Node *> ns(3);
1624           ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
1625           ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
1626           ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
1627           MCAuto<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
1628           for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
1629             {
1630               cPtr[1]=subPtr[offset2+k];
1631               cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
1632             }
1633           int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
1634           if(j!=sz-1)
1635             { cPtr[1]=tmpEnd; }
1636           cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
1637         }
1638       prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
1639       ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
1640     }
1641   if(c->end()!=cPtr)
1642     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
1643   _nodal_connec->decrRef();
1644   _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
1645   addCoo->rearrange(2);
1646   MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
1647   setCoords(coo);
1648   return addCoo->getNumberOfTuples();
1649 }
1650
1651
1652 /// @endcond
1653
1654 /*!
1655  * Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
1656  * returns a result mesh constituted by polygons.
1657  * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
1658  * all nodes from m2.
1659  * The meshes should be in 2D space. In
1660  * addition, returns two arrays mapping cells of the result mesh to cells of the input
1661  * meshes.
1662  *  \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
1663  *                      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)
1664  *  \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
1665  *                      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)
1666  *  \param [in] eps - precision used to detect coincident mesh entities.
1667  *  \param [out] cellNb1 - a new instance of DataArrayInt holding for each result
1668  *         cell an id of the cell of \a m1 it comes from. The caller is to delete
1669  *         this array using decrRef() as it is no more needed.
1670  *  \param [out] cellNb2 - a new instance of DataArrayInt holding for each result
1671  *         cell an id of the cell of \a m2 it comes from. -1 value means that a
1672  *         result cell comes from a cell (or part of cell) of \a m1 not overlapped by
1673  *         any cell of \a m2. The caller is to delete this array using decrRef() as
1674  *         it is no more needed.
1675  *  \return MEDCouplingUMesh * - the result 2D mesh which is a new instance of
1676  *         MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
1677  *         is no more needed.
1678  *  \throw If the coordinates array is not set in any of the meshes.
1679  *  \throw If the nodal connectivity of cells is not defined in any of the meshes.
1680  *  \throw If any of the meshes is not a 2D mesh in 2D space.
1681  *
1682  *  \sa conformize2D, mergeNodes
1683  */
1684 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
1685                                                       double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
1686 {
1687   if(!m1 || !m2)
1688     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
1689   m1->checkFullyDefined();
1690   m2->checkFullyDefined();
1691   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1692   if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
1693     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2  with meshdim equal to 2 and spaceDim equal to 2 too!");
1694
1695   // Step 1: compute all edge intersections (new nodes)
1696   std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1697   MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
1698   DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
1699   std::vector<double> addCoo,addCoordsQuadratic;  // coordinates of newly created nodes
1700   IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
1701                               m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
1702                               addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
1703   revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
1704   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
1705   MCAuto<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
1706
1707   // Step 2: re-order newly created nodes according to the ordering found in m2
1708   std::vector< std::vector<int> > intersectEdge2;
1709   BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
1710   subDiv2.clear(); dd5=0; dd6=0;
1711
1712   // Step 3:
1713   std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
1714   std::vector<int> cNb1,cNb2; //no DataArrayInt because interface with Geometric2D
1715   BuildIntersecting2DCellsFromEdges(eps,m1,desc1->begin(),descIndx1->begin(),intersectEdge1,colinear2,m2,desc2->begin(),descIndx2->begin(),intersectEdge2,addCoo,
1716                                     /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
1717
1718   // Step 4: Prepare final result:
1719   MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1720   addCooDa->alloc((int)(addCoo.size())/2,2);
1721   std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
1722   MCAuto<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
1723   addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
1724   std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
1725   std::vector<const DataArrayDouble *> coordss(4);
1726   coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
1727   MCAuto<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
1728   MCAuto<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
1729   MCAuto<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
1730   MCAuto<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
1731   MCAuto<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
1732   MCAuto<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
1733   ret->setConnectivity(conn,connI,true);
1734   ret->setCoords(coo);
1735   cellNb1=c1.retn(); cellNb2=c2.retn();
1736   return ret.retn();
1737 }
1738
1739 /*!
1740  * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
1741  * 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
1742  * and finally, in case of quadratic polygon the centers of edges new nodes.
1743  * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
1744  *
1745  * \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
1746  *                      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)
1747  * \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
1748  *                      you can invoke orderConsecutiveCells1D on \a mesh1D.
1749  * \param [in] eps - precision used to perform intersections and localization operations.
1750  * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
1751  * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
1752  * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
1753  *                               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.
1754  * \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
1755  *                               and the cell in \a splitMesh2D on the right for the 2nt component. -1 means no cell.
1756  *                               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.
1757  *
1758  * \sa Intersect2DMeshes, orderConsecutiveCells1D, conformize2D, mergeNodes
1759  */
1760 void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
1761 {
1762   if(!mesh2D || !mesh1D)
1763     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
1764   mesh2D->checkFullyDefined();
1765   mesh1D->checkFullyDefined();
1766   const std::vector<std::string>& compNames(mesh2D->getCoords()->getInfoOnComponents());
1767   if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
1768     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
1769   // Step 1: compute all edge intersections (new nodes)
1770   std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
1771   std::vector<double> addCoo,addCoordsQuadratic;  // coordinates of newly created nodes
1772   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1773   //
1774   // Build desc connectivity
1775   DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
1776   MCAuto<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
1777   MCAuto<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
1778   std::map<int,int> mergedNodes;
1779   Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo,mergedNodes);
1780   // use mergeNodes to fix intersectEdge1
1781   for(std::vector< std::vector<int> >::iterator it0=intersectEdge1.begin();it0!=intersectEdge1.end();it0++)
1782     {
1783       std::size_t n((*it0).size()/2);
1784       int eltStart((*it0)[0]),eltEnd((*it0)[2*n-1]);
1785       std::map<int,int>::const_iterator it1;
1786       it1=mergedNodes.find(eltStart);
1787       if(it1!=mergedNodes.end())
1788         (*it0)[0]=(*it1).second;
1789       it1=mergedNodes.find(eltEnd);
1790       if(it1!=mergedNodes.end())
1791         (*it0)[2*n-1]=(*it1).second;
1792     }
1793   //
1794   MCAuto<DataArrayDouble> addCooDa(DataArrayDouble::New());
1795   addCooDa->useArray(&addCoo[0],false,DeallocType::C_DEALLOC,(int)addCoo.size()/2,2);
1796   // Step 2: re-order newly created nodes according to the ordering found in m2
1797   std::vector< std::vector<int> > intersectEdge2;
1798   BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
1799   subDiv2.clear();
1800   // Step 3: compute splitMesh1D
1801   MCAuto<DataArrayInt> idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear;
1802   MCAuto<DataArrayInt> ret2(DataArrayInt::New()); ret2->alloc(0,1);
1803   MCAuto<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo,mergedNodes,colinear2,intersectEdge1,
1804       idsInRet1Colinear,idsInDescMesh2DForIdsInRetColinear));
1805   MCAuto<DataArrayInt> ret3(DataArrayInt::New()); ret3->alloc(ret1->getNumberOfCells()*2,1); ret3->fillWithValue(std::numeric_limits<int>::max()); ret3->rearrange(2);
1806   MCAuto<DataArrayInt> idsInRet1NotColinear(idsInRet1Colinear->buildComplement(ret1->getNumberOfCells()));
1807   // deal with cells in mesh2D that are not cut but only some of their edges are
1808   MCAuto<DataArrayInt> idsInDesc2DToBeRefined(idsInDescMesh2DForIdsInRetColinear->deepCopy());
1809   idsInDesc2DToBeRefined->abs(); idsInDesc2DToBeRefined->applyLin(1,-1);
1810   idsInDesc2DToBeRefined=idsInDesc2DToBeRefined->buildUnique();
1811   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
1812   if(!idsInDesc2DToBeRefined->empty())
1813     {
1814       DataArrayInt *out0(0),*outi0(0);
1815       DataArrayInt::ExtractFromIndexedArrays(idsInDesc2DToBeRefined->begin(),idsInDesc2DToBeRefined->end(),dd3,dd4,out0,outi0);
1816       MCAuto<DataArrayInt> outi0s(outi0);
1817       out0s=out0;
1818       out0s=out0s->buildUnique();
1819       out0s->sort(true);
1820     }
1821   //
1822   MCAuto<MEDCouplingUMesh> ret1NonCol(static_cast<MEDCouplingUMesh *>(ret1->buildPartOfMySelf(idsInRet1NotColinear->begin(),idsInRet1NotColinear->end())));
1823   MCAuto<DataArrayDouble> baryRet1(ret1NonCol->computeCellCenterOfMass());
1824   MCAuto<DataArrayInt> elts,eltsIndex;
1825   mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
1826   MCAuto<DataArrayInt> eltsIndex2(DataArrayInt::New()); eltsIndex2->alloc(0,1);
1827   if (eltsIndex->getNumberOfTuples() > 1)
1828     eltsIndex2 = eltsIndex->deltaShiftIndex();
1829   MCAuto<DataArrayInt> eltsIndex3(eltsIndex2->findIdsEqual(1));
1830   if(eltsIndex2->count(0)+eltsIndex3->getNumberOfTuples()!=ret1NonCol->getNumberOfCells())
1831     throw INTERP_KERNEL::Exception("Intersect2DMeshWith1DLine : internal error 1 !");
1832   MCAuto<DataArrayInt> cellsToBeModified(elts->buildUnique());
1833   MCAuto<DataArrayInt> untouchedCells(cellsToBeModified->buildComplement(mesh2D->getNumberOfCells()));
1834   if((DataArrayInt *)out0s)
1835     untouchedCells=untouchedCells->buildSubstraction(out0s);//if some edges in ret1 are colinear to descending mesh of mesh2D remove cells from untouched one
1836   std::vector< MCAuto<MEDCouplingUMesh> > outMesh2DSplit;
1837   // OK all is ready to insert in ret2 mesh
1838   if(!untouchedCells->empty())
1839     {// the most easy part, cells in mesh2D not impacted at all
1840       outMesh2DSplit.push_back(static_cast<MEDCouplingUMesh *>(mesh2D->buildPartOfMySelf(untouchedCells->begin(),untouchedCells->end())));
1841       outMesh2DSplit.back()->setCoords(ret1->getCoords());
1842       ret2->pushBackValsSilent(untouchedCells->begin(),untouchedCells->end());
1843     }
1844   if((DataArrayInt *)out0s)
1845     {// here dealing with cells in out0s but not in cellsToBeModified
1846       MCAuto<DataArrayInt> fewModifiedCells(out0s->buildSubstraction(cellsToBeModified));
1847       const int *rdptr(dd3->begin()),*rdiptr(dd4->begin()),*dptr(dd1->begin()),*diptr(dd2->begin());
1848       for(const int *it=fewModifiedCells->begin();it!=fewModifiedCells->end();it++)
1849         {
1850           outMesh2DSplit.push_back(BuildRefined2DCell(ret1->getCoords(),mesh2D,*it,dptr+diptr[*it],dptr+diptr[*it+1],intersectEdge1));
1851           ret1->setCoords(outMesh2DSplit.back()->getCoords());
1852         }
1853       int offset(ret2->getNumberOfTuples());
1854       ret2->pushBackValsSilent(fewModifiedCells->begin(),fewModifiedCells->end());
1855       MCAuto<DataArrayInt> partOfRet3(DataArrayInt::New()); partOfRet3->alloc(2*idsInRet1Colinear->getNumberOfTuples(),1);
1856       partOfRet3->fillWithValue(std::numeric_limits<int>::max()); partOfRet3->rearrange(2);
1857       int kk(0),*ret3ptr(partOfRet3->getPointer());
1858       for(const int *it=idsInDescMesh2DForIdsInRetColinear->begin();it!=idsInDescMesh2DForIdsInRetColinear->end();it++,kk++)
1859         {
1860           int faceId(std::abs(*it)-1);
1861           for(const int *it2=rdptr+rdiptr[faceId];it2!=rdptr+rdiptr[faceId+1];it2++)
1862             {
1863               int tmp(fewModifiedCells->findIdFirstEqual(*it2));
1864               if(tmp!=-1)
1865                 {
1866                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1867                     ret3ptr[2*kk]=tmp+offset;
1868                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1869                     ret3ptr[2*kk+1]=tmp+offset;
1870                 }
1871               else
1872                 {//the current edge is shared by a 2D cell that will be split just after
1873                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],-(*it))!=dptr+diptr[*it2+1])
1874                     ret3ptr[2*kk]=-(*it2+1);
1875                   if(std::find(dptr+diptr[*it2],dptr+diptr[*it2+1],(*it))!=dptr+diptr[*it2+1])
1876                     ret3ptr[2*kk+1]=-(*it2+1);
1877                 }
1878             }
1879         }
1880       m1Desc->setCoords(ret1->getCoords());
1881       ret1NonCol->setCoords(ret1->getCoords());
1882       ret3->setPartOfValues3(partOfRet3,idsInRet1Colinear->begin(),idsInRet1Colinear->end(),0,2,1,true);
1883       if(!outMesh2DSplit.empty())
1884         {
1885           DataArrayDouble *da(outMesh2DSplit.back()->getCoords());
1886           for(std::vector< MCAuto<MEDCouplingUMesh> >::iterator itt=outMesh2DSplit.begin();itt!=outMesh2DSplit.end();itt++)
1887             (*itt)->setCoords(da);
1888         }
1889     }
1890   cellsToBeModified=cellsToBeModified->buildUniqueNotSorted();
1891   for(const int *it=cellsToBeModified->begin();it!=cellsToBeModified->end();it++)
1892     {
1893       MCAuto<DataArrayInt> idsNonColPerCell(elts->findIdsEqual(*it));
1894       idsNonColPerCell->transformWithIndArr(eltsIndex3->begin(),eltsIndex3->end());
1895       MCAuto<DataArrayInt> idsNonColPerCell2(idsInRet1NotColinear->selectByTupleId(idsNonColPerCell->begin(),idsNonColPerCell->end()));
1896       MCAuto<MEDCouplingUMesh> partOfMesh1CuttingCur2DCell(static_cast<MEDCouplingUMesh *>(ret1NonCol->buildPartOfMySelf(idsNonColPerCell->begin(),idsNonColPerCell->end())));
1897       MCAuto<DataArrayInt> partOfRet3;
1898       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));
1899       ret3->setPartOfValues3(partOfRet3,idsNonColPerCell2->begin(),idsNonColPerCell2->end(),0,2,1,true);
1900       outMesh2DSplit.push_back(splitOfOneCell);
1901       for(std::size_t i=0;i<splitOfOneCell->getNumberOfCells();i++)
1902         ret2->pushBackSilent(*it);
1903     }
1904   //
1905   std::size_t nbOfMeshes(outMesh2DSplit.size());
1906   std::vector<const MEDCouplingUMesh *> tmp(nbOfMeshes);
1907   for(std::size_t i=0;i<nbOfMeshes;i++)
1908     tmp[i]=outMesh2DSplit[i];
1909   //
1910   ret1->getCoords()->setInfoOnComponents(compNames);
1911   MCAuto<MEDCouplingUMesh> ret2D(MEDCouplingUMesh::MergeUMeshesOnSameCoords(tmp));
1912   // To finish - filter ret3 - std::numeric_limits<int>::max() -> -1 - negate values must be resolved.
1913   ret3->rearrange(1);
1914   MCAuto<DataArrayInt> edgesToDealWith(ret3->findIdsStrictlyNegative());
1915   for(const int *it=edgesToDealWith->begin();it!=edgesToDealWith->end();it++)
1916     {
1917       int old2DCellId(-ret3->getIJ(*it,0)-1);
1918       MCAuto<DataArrayInt> candidates(ret2->findIdsEqual(old2DCellId));
1919       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
1920     }
1921   ret3->changeValue(std::numeric_limits<int>::max(),-1);
1922   ret3->rearrange(2);
1923   //
1924   splitMesh1D=ret1.retn();
1925   splitMesh2D=ret2D.retn();
1926   cellIdInMesh2D=ret2.retn();
1927   cellIdInMesh1D=ret3.retn();
1928 }
1929
1930 /*!
1931  * \b WARNING this method is \b potentially \b non \b const (if returned array is empty).
1932  * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
1933  * 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
1934  * 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).
1935  * 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.
1936  *
1937  * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
1938  * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
1939  *
1940  * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 too.
1941  * This method expects that all nodes in \a this are not closer than \a eps.
1942  * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
1943  *
1944  * \param [in] eps the relative error to detect merged edges.
1945  * \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
1946  *                           that the user is expected to deal with.
1947  *
1948  * \throw If \a this is not coherent.
1949  * \throw If \a this has not spaceDim equal to 2.
1950  * \throw If \a this has not meshDim equal to 2.
1951  * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
1952  */
1953 DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
1954 {
1955   static const int SPACEDIM=2;
1956   checkConsistencyLight();
1957   if(getSpaceDimension()!=2 || getMeshDimension()!=2)
1958     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
1959   MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
1960   MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
1961   const int *c(mDesc->getNodalConnectivity()->begin()),*ci(mDesc->getNodalConnectivityIndex()->begin()),*rd(revDesc1->begin()),*rdi(revDescIndx1->begin());
1962   MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
1963   const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
1964   int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
1965   std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
1966   std::vector<double> addCoo;
1967   BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
1968   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
1969   for(int i=0;i<nDescCell;i++)
1970     {
1971       std::vector<int> candidates;
1972       myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
1973       for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
1974         if(*it>i)  // we're dealing with pair of edges, no need to treat the same pair twice
1975           {
1976             std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
1977             INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
1978                 *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
1979             INTERP_KERNEL::MergePoints merge;
1980             INTERP_KERNEL::QuadraticPolygon c1,c2;
1981             e1->intersectWith(e2,merge,c1,c2);
1982             e1->decrRef(); e2->decrRef();
1983             if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
1984               overlapEdge[i].push_back(*it);
1985             if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
1986               overlapEdge[*it].push_back(i);
1987           }
1988     }
1989   // splitting done. sort intersect point in intersectEdge.
1990   std::vector< std::vector<int> > middle(nDescCell);
1991   int nbOf2DCellsToBeSplit(0);
1992   bool middleNeedsToBeUsed(false);
1993   std::vector<bool> cells2DToTreat(nDescCell,false);
1994   for(int i=0;i<nDescCell;i++)
1995     {
1996       std::vector<int>& isect(intersectEdge[i]);
1997       int sz((int)isect.size());
1998       if(sz>1)
1999         {
2000           std::map<MCAuto<INTERP_KERNEL::Node>,int> m;
2001           INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
2002           e->sortSubNodesAbs(coords,isect);
2003           e->decrRef();
2004         }
2005       if(sz!=0)
2006         {
2007           int idx0(rdi[i]),idx1(rdi[i+1]);
2008           if(idx1-idx0!=1)
2009             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
2010           if(!cells2DToTreat[rd[idx0]])
2011             {
2012               cells2DToTreat[rd[idx0]]=true;
2013               nbOf2DCellsToBeSplit++;
2014             }
2015           // try to reuse at most eventual 'middle' of SEG3
2016           std::vector<int>& mid(middle[i]);
2017           mid.resize(sz+1,-1);
2018           if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
2019             {
2020               middleNeedsToBeUsed=true;
2021               const std::vector<int>& candidates(overlapEdge[i]);
2022               std::vector<int> trueCandidates;
2023               for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
2024                 if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
2025                   trueCandidates.push_back(*itc);
2026               int stNode(c[ci[i]+1]),endNode(isect[0]);
2027               for(int j=0;j<sz+1;j++)
2028                 {
2029                   for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
2030                     {
2031                       int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
2032                       if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
2033                         { mid[j]=*itc; break; }
2034                     }
2035                   stNode=endNode;
2036                   endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
2037                 }
2038             }
2039         }
2040     }
2041   MCAuto<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
2042   if(nbOf2DCellsToBeSplit==0)
2043     return ret.retn();
2044   //
2045   int *retPtr(ret->getPointer());
2046   for(int i=0;i<nCell;i++)
2047     if(cells2DToTreat[i])
2048       *retPtr++=i;
2049   //
2050   MCAuto<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
2051   DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
2052   DataArrayInt::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
2053   DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
2054   if(middleNeedsToBeUsed)
2055     { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
2056   MCAuto<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
2057   int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
2058   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.
2059   setPartOfMySelf(ret->begin(),ret->end(),*modif);
2060   {
2061     bool areNodesMerged; int newNbOfNodes;
2062     if(nbOfNodesCreated!=0)
2063       MCAuto<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
2064   }
2065   return ret.retn();
2066 }
2067
2068 /*!
2069  * 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.
2070  * 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).
2071  * 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
2072  * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
2073  * 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
2074  * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
2075  *
2076  * 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
2077  * using new instance, idem for coordinates.
2078  *
2079  * 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.
2080  *
2081  * \return DataArrayInt  * - The list of cellIds in \a this that have at least one edge colinearized.
2082  *
2083  * \throw If \a this is not coherent.
2084  * \throw If \a this has not spaceDim equal to 2.
2085  * \throw If \a this has not meshDim equal to 2.
2086  *
2087  * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
2088  */
2089 DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
2090 {
2091   return internalColinearize2D(eps, false);
2092 }
2093
2094 /*!
2095  * Performs exactly the same job as colinearize2D, except that this function does not create new non-conformal cells.
2096  * 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
2097  * merged, contrary to colinearize2D().
2098  *
2099  * \sa MEDCouplingUMesh::colinearize2D
2100  */
2101 DataArrayInt *MEDCouplingUMesh::colinearizeKeepingConform2D(double eps)
2102 {
2103   return internalColinearize2D(eps, true);
2104 }
2105
2106
2107 /*!
2108  * \param stayConform is set to True, will not fuse two edges sharing a node that has (strictly) more than 2 egdes connected to it
2109  */
2110 DataArrayInt *MEDCouplingUMesh::internalColinearize2D(double eps, bool stayConform)
2111 {
2112   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2113   checkConsistencyLight();
2114   if(getSpaceDimension()!=2 || getMeshDimension()!=2)
2115     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
2116   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
2117   int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
2118   const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
2119   MCAuto<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
2120   std::map<int, bool> forbiddenPoints;  // list of points that can not be removed (or it will break conformity)
2121   if(stayConform) //
2122     {
2123       // A point that is used by more than 2 edges can not be removed without breaking conformity:
2124       MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descI1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescI1(DataArrayInt::New());
2125       MCAuto<MEDCouplingUMesh> mDesc1D(buildDescendingConnectivity(desc1,descI1,revDesc1,revDescI1));
2126       MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2127       MCAuto<MEDCouplingUMesh> mDesc0D(mDesc1D->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2));
2128       MCAuto<DataArrayInt> dsi(revDescI2->deltaShiftIndex());
2129       MCAuto<DataArrayInt> ids(dsi->findIdsGreaterThan(2));
2130       const int * cPtr(mDesc0D->getNodalConnectivity()->begin());
2131       for(auto it = ids->begin(); it != ids->end(); it++)
2132          forbiddenPoints[cPtr[2*(*it)+1]] = true;  // we know that a 0D mesh has a connectivity of the form [NORM_POINT1, i1, NORM_POINT1, i2, ...]
2133     }
2134
2135   MCAuto<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
2136   const double *coords(_coords->begin());
2137   int *newciptr(newci->getPointer());
2138   for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
2139     {
2140       if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,forbiddenPoints, /*out*/ newc,appendedCoords))
2141         ret->pushBackSilent(i);
2142       newciptr[1]=newc->getNumberOfTuples();
2143     }
2144   //
2145   if(ret->empty())
2146     return ret.retn();
2147   if(!appendedCoords->empty())
2148     {
2149       appendedCoords->rearrange(2);
2150       MCAuto<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
2151       //non const part
2152       setCoords(newCoords);
2153     }
2154   //non const part
2155   setConnectivity(newc,newci,true);
2156   return ret.retn();
2157 }
2158
2159 ///@cond INTERNAL
2160 /**
2161  * c, cI describe a wire mesh in 3D space, in local numbering
2162  * startNode, endNode in global numbering
2163  *\return true if the segment is indeed split
2164  */
2165 bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
2166                                             const int * c, const int * cI, const int *idsBg, const int *endBg,
2167                                             std::vector<int> & pointIds, std::vector<int> & hitSegs)
2168 {
2169   using namespace std;
2170
2171   const int SPACEDIM=3;
2172   typedef pair<double, int> PairDI;
2173   set< PairDI > x;
2174   for (const int * it = idsBg; it != endBg; ++it)
2175     {
2176       assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
2177       int start = c[cI[*it]+1], end = c[cI[*it]+2];
2178       x.insert(make_pair(coo[start*SPACEDIM], start));  // take only X coords
2179       x.insert(make_pair(coo[end*SPACEDIM], end));
2180     }
2181
2182   vector<PairDI> xx(x.begin(), x.end());
2183   sort(xx.begin(),xx.end());
2184   pointIds.reserve(xx.size());
2185
2186   // Keep what is inside [startNode, endNode]:
2187   int go = 0;
2188   for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
2189     {
2190       const int idx = (*it).second;
2191       if (!go)
2192         {
2193           if (idx == startNode)   go = 1;
2194           if (idx == endNode)     go = 2;
2195           if (go)                 pointIds.push_back(idx);
2196           continue;
2197         }
2198       pointIds.push_back(idx);
2199       if (idx == endNode || idx == startNode)
2200         break;
2201     }
2202
2203 //  vector<int> pointIds2(pointIds.size()+2);
2204 //  copy(pointIds.begin(), pointIds.end(), pointIds2.data()+1);
2205 //  pointIds2[0] = startNode;
2206 //  pointIds2[pointIds2.size()-1] = endNode;
2207
2208   if (go == 2)
2209     reverse(pointIds.begin(), pointIds.end());
2210
2211   // Now identify smaller segments that are not sub-divided - those won't need any further treatment:
2212   for (const int * it = idsBg; it != endBg; ++it)
2213     {
2214       int start = c[cI[*it]+1], end = c[cI[*it]+2];
2215       vector<int>::const_iterator itStart = find(pointIds.begin(), pointIds.end(), start);
2216       if (itStart == pointIds.end()) continue;
2217       vector<int>::const_iterator itEnd = find(pointIds.begin(), pointIds.end(), end);
2218       if (itEnd == pointIds.end())               continue;
2219       if (abs(distance(itEnd, itStart)) != 1)    continue;
2220       hitSegs.push_back(*it);   // segment is undivided.
2221     }
2222
2223   return (pointIds.size() > 2); // something else apart start and end node
2224 }
2225
2226 void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
2227                                           const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
2228 {
2229   using namespace std;
2230   int dst = distance(sIdxConn, sIdxConnE);
2231   modifiedFace.reserve(dst + insidePoints.size()-2);
2232   modifiedFace.resize(dst);
2233   copy(sIdxConn, sIdxConnE, modifiedFace.data());
2234
2235   vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
2236   vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
2237   if (startPos == shortEnd)
2238     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2239   vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
2240   if (endPos == shortEnd)
2241     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
2242   int d = distance(startPos, endPos);
2243   if (d == 1 || d == (1-dst)) // don't use modulo, for neg numbers, result is implementation defined ...
2244     modifiedFace.insert(++startPos, ++insidePoints.begin(), --insidePoints.end());  // insidePoints also contains start and end node. Those don't need to be inserted.
2245   else
2246     modifiedFace.insert(++endPos, ++insidePoints.rbegin(), --insidePoints.rend());
2247 }
2248
2249 ///@endcond
2250
2251
2252 /*!
2253  * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
2254  * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
2255  * This method performs a conformization of \b this.
2256  *
2257  * Only polyhedron cells are supported. You can call convertAllToPoly()
2258  *
2259  * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
2260  * This method expects that all nodes in \a this are not closer than \a eps.
2261  * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
2262  *
2263  * \param [in] eps the relative error to detect merged edges.
2264  * \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
2265  *                           that the user is expected to deal with.
2266  *
2267  * \throw If \a this is not coherent.
2268  * \throw If \a this has not spaceDim equal to 3.
2269  * \throw If \a this has not meshDim equal to 3.
2270  * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
2271  */
2272 DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
2273 {
2274   using namespace std;
2275
2276   static const int SPACEDIM=3;
2277   checkConsistencyLight();
2278   if(getSpaceDimension()!=3 || getMeshDimension()!=3)
2279     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
2280   if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
2281     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
2282
2283   MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
2284   const double * coo(_coords->begin());
2285   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
2286
2287   {
2288     /*************************
2289      *  STEP 1  -- faces
2290      *************************/
2291     MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2292     MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
2293     const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2294     const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2295     MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
2296
2297     // Build BBTree
2298     MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree(eps));
2299     const double *bbox(bboxArr->begin()); getCoords()->begin();
2300     int nDescCell(mDesc->getNumberOfCells());
2301     BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
2302     // Surfaces - handle biggest first
2303     MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
2304     DataArrayDouble * surfs = surfF->getArray();
2305     // Normal field
2306     MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
2307     DataArrayDouble * normals = normalsF->getArray();
2308     const double * normalsP = normals->getConstPointer();
2309
2310     // Sort faces by decreasing surface:
2311     vector< pair<double,int> > S;
2312     for(std::size_t i=0;i < surfs->getNumberOfTuples();i++)
2313       {
2314         pair<double,int> p = make_pair(surfs->begin()[i], i);
2315         S.push_back(p);
2316       }
2317     sort(S.rbegin(),S.rend()); // reverse sort
2318     vector<bool> hit(nDescCell);
2319     fill(hit.begin(), hit.end(), false);
2320     vector<int> hitPoly; // the final result: which 3D cells have been modified.
2321
2322     for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2323       {
2324         int faceIdx = (*it).second;
2325         if (hit[faceIdx]) continue;
2326
2327         vector<int> candidates, cands2;
2328         myTree.getIntersectingElems(bbox+faceIdx*2*SPACEDIM,candidates);
2329         // Keep only candidates whose normal matches the normal of current face
2330         for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2331           {
2332             bool col = INTERP_KERNEL::isColinear3D(normalsP + faceIdx*SPACEDIM, normalsP + *(it2)*SPACEDIM, eps);
2333             if (*it2 != faceIdx && col)
2334               cands2.push_back(*it2);
2335           }
2336         if (!cands2.size())
2337           continue;
2338
2339         // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
2340         INTERP_KERNEL::TranslationRotationMatrix rotation;
2341         INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+1]),
2342                                                                    coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+2]),
2343                                                                    coo+SPACEDIM*(cDesc[cIDesc[faceIdx]+3]), rotation);
2344
2345         MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(faceIdx, faceIdx+1,1,false));  // false=zipCoords is called
2346         MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
2347         double * cooPartRef(mPartRef->_coords->getPointer());
2348         double * cooPartCand(mPartCand->_coords->getPointer());
2349         for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2350           rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2351         for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2352           rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2353
2354         // Localize faces in 2D thanks to barycenters
2355         MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2356         vector<int> compo; compo.push_back(2);
2357         MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2358         MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(-eps, +eps);
2359         if (!idsGoodPlane->getNumberOfTuples())
2360           continue;
2361
2362         baryPart = baryPart->selectByTupleId(*idsGoodPlane);
2363
2364         compo[0] = 0; compo.push_back(1);
2365         MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
2366         mPartRef->changeSpaceDimension(2,0.0);
2367         MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
2368         mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
2369
2370         if (!cc->getNumberOfTuples())
2371           continue;
2372         MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
2373
2374         {
2375           MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
2376           if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
2377             {
2378               ostringstream oss;
2379               oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2380               throw INTERP_KERNEL::Exception(oss.str());
2381             }
2382         }
2383
2384         MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
2385         // Boundary face:
2386         if (!ids->getNumberOfTuples())
2387           continue;
2388
2389         double checkSurf=0.0;
2390         const int * idsGoodPlaneP(idsGoodPlane->begin());
2391         for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2392           {
2393             int faceIdx2 = cands2[idsGoodPlaneP[*ii]];
2394             hit[faceIdx2] = true;
2395             checkSurf += surfs->begin()[faceIdx2];
2396           }
2397         if (fabs(checkSurf - surfs->begin()[faceIdx]) > eps)
2398           {
2399             ostringstream oss;
2400             oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << faceIdx << " violates this condition!";
2401             throw INTERP_KERNEL::Exception(oss.str());
2402           }
2403
2404         // For all polyhedrons using this face, replace connectivity:
2405         vector<int> polyIndices, packsIds, facePack;
2406         for (int ii=revDescIP[faceIdx]; ii < revDescIP[faceIdx+1]; ii++)
2407           polyIndices.push_back(revDescP[ii]);
2408         ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
2409
2410         // Current face connectivity
2411         const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2412         const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2413         connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
2414         // Deletion of old faces
2415         int jj=0;
2416         for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2, ++jj)
2417           {
2418             if (packsIds[jj] == -1)
2419               // The below should never happen - if a face is used several times, with a different layout of the nodes
2420               // it means that it is already conform, so it is *not* hit by the algorithm. The algorithm only hits
2421               // faces which are actually used only once, by a single cell. This is different for edges below.
2422               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
2423             else
2424               connSla->deletePack(*it2, packsIds[jj]);
2425           }
2426         // Insertion of new faces:
2427         for (const int * ii = ids->begin(); ii != ids->end(); ii++)
2428           {
2429             // Build pack from the face to insert:
2430             int faceIdx2 = cands2[idsGoodPlane->getIJ(*ii,0)];
2431             int facePack2Sz;
2432             const int * facePack2 = connSlaDesc->getSimplePackSafePtr(faceIdx2, facePack2Sz); // contains the type!
2433             // Insert it in all hit polyhedrons:
2434             for (vector<int>::const_iterator it2=polyIndices.begin(); it2!=polyIndices.end(); ++it2)
2435               connSla->pushBackPack(*it2, facePack2+1, facePack2+facePack2Sz);  // without the type
2436           }
2437       }
2438   }  // end step1
2439
2440   // Set back modified connectivity
2441   MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
2442   MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
2443   connSla->convertToPolyhedronConn(cAuto, cIAuto);
2444
2445   {
2446     /************************
2447      *  STEP 2 -- edges
2448      ************************/
2449     // Now we have a face-conform mesh.
2450
2451     // Recompute descending
2452     MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
2453     // Rebuild desc connectivity with orientation this time!!
2454     MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc,descI,revDesc,revDescI));
2455     const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
2456     const int *descIP(descI->getConstPointer()), *descP(desc->getConstPointer());
2457     const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
2458     MCAuto<DataArrayInt> ciDeepC(mDesc->getNodalConnectivityIndex()->deepCopy());
2459     MCAuto<DataArrayInt> cDeepC(mDesc->getNodalConnectivity()->deepCopy());
2460     MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(ciDeepC, cDeepC));
2461     MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
2462     MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
2463 //    std::cout << "writing!\n";
2464 //    mDesc->writeVTK("/tmp/toto_desc_confInter.vtu");
2465 //    mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu");
2466     const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
2467     const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
2468     MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree(eps));
2469     const double *bbox2(bboxArr->begin());
2470     int nDesc2Cell=mDesc2->getNumberOfCells();
2471     BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
2472
2473     // Edges - handle longest first
2474     MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
2475     DataArrayDouble * lens = lenF->getArray();
2476
2477     // Sort edges by decreasing length:
2478     vector<pair<double,int> > S;
2479     for(std::size_t i=0;i < lens->getNumberOfTuples();i++)
2480       {
2481         pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
2482         S.push_back(p);
2483       }
2484     sort(S.rbegin(),S.rend()); // reverse sort
2485
2486     vector<bool> hit(nDesc2Cell);
2487     fill(hit.begin(), hit.end(), false);
2488
2489     for( vector<pair<double,int> >::const_iterator it = S.begin(); it != S.end(); it++)
2490       {
2491         int eIdx = (*it).second;
2492         if (hit[eIdx])
2493           continue;
2494
2495         vector<int> candidates, cands2;
2496         myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
2497         // Keep only candidates colinear with current edge
2498         double vCurr[3];
2499         unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
2500         for (int i3=0; i3 < 3; i3++)  // TODO: use fillSonCellNodalConnectivity2 or similar?
2501           vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
2502         for(vector<int>::const_iterator it2=candidates.begin();it2!=candidates.end();it2++)
2503           {
2504             double vOther[3];
2505             unsigned start2 = cDesc2[cIDesc2[*it2]+1], end2 = cDesc2[cIDesc2[*it2]+2];
2506             for (int i3=0; i3 < 3; i3++)
2507               vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
2508             bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
2509             // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
2510             // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
2511             if (col)
2512               cands2.push_back(*it2);
2513           }
2514         if (cands2.size() == 1 && cands2[0] == eIdx)  // see warning above
2515           continue;
2516
2517         // Now rotate edges to bring them on Ox
2518         int startNode = cDesc2[cIDesc2[eIdx]+1];
2519         int endNode = cDesc2[cIDesc2[eIdx]+2];
2520         INTERP_KERNEL::TranslationRotationMatrix rotation;
2521         INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
2522         MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false));  // false=zipCoords is called
2523         MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
2524         MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
2525         int nbElemsNotM1;
2526         {
2527           MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
2528           nbElemsNotM1 = tmp->getNbOfElems();
2529         }
2530         MCAuto<DataArrayInt>  nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
2531         double * cooPartRef(mPartRef->_coords->getPointer());
2532         double * cooPartCand(mPartCand->_coords->getPointer());
2533         for (std::size_t ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
2534           rotation.transform_vector(cooPartRef+SPACEDIM*ii);
2535         for (std::size_t ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
2536           rotation.transform_vector(cooPartCand+SPACEDIM*ii);
2537
2538
2539         // Eliminate all edges for which y or z is not null
2540         MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
2541         vector<int> compo; compo.push_back(1);
2542         MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
2543         compo[0] = 2;
2544         MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
2545         MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
2546         MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
2547         MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
2548         if (!idsGoodLine->getNumberOfTuples())
2549           continue;
2550
2551         // Now the ordering along the Ox axis:
2552         std::vector<int> insidePoints, hitSegs;
2553         bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
2554             mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
2555             idsGoodLine->begin(), idsGoodLine->end(),
2556             /*out*/insidePoints, hitSegs);
2557         // Optim: smaller segments completely included in eIdx and not split won't need any further treatment:
2558         for (vector<int>::const_iterator its=hitSegs.begin(); its != hitSegs.end(); ++its)
2559           hit[cands2[*its]] = true;
2560
2561         if (!isSplit)  // current segment remains in one piece
2562           continue;
2563
2564         // Get original node IDs in global coords array
2565         for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
2566           *iit = nodeMapInv->begin()[*iit];
2567
2568         vector<int> polyIndices, packsIds, facePack;
2569         // For each face implying this edge
2570         for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
2571           {
2572             int faceIdx = revDescP2[ii];
2573             // each cell where this face is involved connectivity will be modified:
2574             ret->pushBackValsSilent(revDescP + revDescIP[faceIdx], revDescP + revDescIP[faceIdx+1]);
2575
2576             // Current face connectivity
2577             const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
2578             const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
2579
2580             vector<int> modifiedFace;
2581             ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
2582             modifiedFace.insert(modifiedFace.begin(), INTERP_KERNEL::NORM_POLYGON);
2583             connSlaDesc->replaceSimplePack(faceIdx, modifiedFace.data(), modifiedFace.data()+modifiedFace.size());
2584           }
2585       }
2586
2587     // Rebuild 3D connectivity from descending:
2588     MCAuto<MEDCouplingSkyLineArray> newConn(MEDCouplingSkyLineArray::New());
2589     MCAuto<DataArrayInt> superIdx(DataArrayInt::New());  superIdx->alloc(getNumberOfCells()+1); superIdx->fillWithValue(0);
2590     MCAuto<DataArrayInt> idx(DataArrayInt::New());       idx->alloc(1);                         idx->fillWithValue(0);
2591     MCAuto<DataArrayInt> vals(DataArrayInt::New());      vals->alloc(0);
2592     newConn->set3(superIdx, idx, vals);
2593     for(std::size_t ii = 0; ii < getNumberOfCells(); ii++)
2594       for (int jj=descIP[ii]; jj < descIP[ii+1]; jj++)
2595         {
2596           int sz, faceIdx = abs(descP[jj])-1;
2597           bool orient = descP[jj]>0;
2598           const int * p = connSlaDesc->getSimplePackSafePtr(faceIdx, sz);
2599           if (orient)
2600             newConn->pushBackPack(ii, p+1, p+sz);  // +1 to skip type
2601           else
2602             {
2603               vector<int> rev(sz-1);
2604               for (int kk=0; kk<sz-1; kk++) rev[kk]=*(p+sz-kk-1);
2605               newConn->pushBackPack(ii, rev.data(), rev.data()+sz-1);
2606             }
2607         }
2608     // And finally:
2609     newConn->convertToPolyhedronConn(cAuto, cIAuto);
2610   } // end step2
2611
2612   ret = ret->buildUniqueNotSorted();
2613   return ret.retn();
2614 }
2615
2616