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