src/MEDCoupling/MEDCouplingUMesh_intersection.cxx

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