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