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