src/MEDCoupling/MEDCouplingUMesh_intersection.cxx

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