]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/MEDCoupling/MEDCouplingUMesh_internal.cxx
Salome HOME
10c0858400c273f6f5907d74884e262e7bf4bbbd
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh_internal.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 "MEDCouplingCMesh.hxx"
23 #include "MEDCoupling1GTUMesh.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
25 #include "MEDCouplingSkyLineArray.hxx"
26 #include "CellModel.hxx"
27 #include "VolSurfUser.txx"
28 #include "InterpolationUtils.hxx"
29 #include "PointLocatorAlgos.txx"
30 #include "BBTree.txx"
31 #include "BBTreeDst.txx"
32 #include "SplitterTetra.hxx"
33 #include "DiameterCalculator.hxx"
34 #include "DirectedBoundingBox.hxx"
35 #include "InterpKernelMatrixTools.hxx"
36 #include "InterpKernelMeshQuality.hxx"
37 #include "InterpKernelCellSimplify.hxx"
38 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
39 #include "InterpKernelAutoPtr.hxx"
40 #include "InterpKernelGeo2DNode.hxx"
41 #include "InterpKernelGeo2DEdgeLin.hxx"
42 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
43 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
44 #include "MEDCouplingUMesh_internal.hxx"
45
46 #include <sstream>
47 #include <fstream>
48 #include <numeric>
49 #include <cstring>
50 #include <limits>
51 #include <list>
52
53 using namespace MEDCoupling;
54
55 /*!
56  * This method checks that all arrays are set. If yes nothing done if no an exception is thrown.
57  */
58 void MEDCouplingUMesh::checkFullyDefined() const
59 {
60   if(!_nodal_connec_index || !_nodal_connec || !_coords)
61     throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity and coordinates set in unstructured mesh.");
62 }
63
64 /*!
65  * This method checks that all connectivity arrays are set. If yes nothing done if no an exception is thrown.
66  */
67 void MEDCouplingUMesh::checkConnectivityFullyDefined() const
68 {
69   if(!_nodal_connec_index || !_nodal_connec)
70     throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity set in unstructured mesh.");
71 }
72
73 void MEDCouplingUMesh::reprConnectivityOfThisLL(std::ostringstream& stream) const
74 {
75   if(_nodal_connec!=0 && _nodal_connec_index!=0)
76     {
77       int nbOfCells=getNumberOfCells();
78       const int *c=_nodal_connec->getConstPointer();
79       const int *ci=_nodal_connec_index->getConstPointer();
80       for(int i=0;i<nbOfCells;i++)
81         {
82           const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[i]]);
83           stream << "Cell #" << i << " " << cm.getRepr() << " : ";
84           std::copy(c+ci[i]+1,c+ci[i+1],std::ostream_iterator<int>(stream," "));
85           stream << "\n";
86         }
87     }
88   else
89     stream << "Connectivity not defined !\n";
90 }
91
92
93 /*!
94  * This method implements policy 0 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
95  */
96 DataArrayInt *MEDCouplingUMesh::simplexizePol0()
97 {
98   checkConnectivityFullyDefined();
99   if(getMeshDimension()!=2)
100     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePol0 : this policy is only available for mesh with meshdim == 2 !");
101   int nbOfCells=getNumberOfCells();
102   MCAuto<DataArrayInt> ret=DataArrayInt::New();
103   int nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
104   ret->alloc(nbOfCells+nbOfCutCells,1);
105   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
106   int *retPt=ret->getPointer();
107   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
108   MCAuto<DataArrayInt> newConnI=DataArrayInt::New();
109   newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
110   newConn->alloc(getNodalConnectivityArrayLen()+3*nbOfCutCells,1);
111   int *pt=newConn->getPointer();
112   int *ptI=newConnI->getPointer();
113   ptI[0]=0;
114   const int *oldc=_nodal_connec->begin();
115   const int *ci=_nodal_connec_index->begin();
116   for(int i=0;i<nbOfCells;i++,ci++)
117     {
118       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
119         {
120           const int tmp[8]={(int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+3],
121             (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+3],oldc[ci[0]+4]};
122           pt=std::copy(tmp,tmp+8,pt);
123           ptI[1]=ptI[0]+4;
124           ptI[2]=ptI[0]+8;
125           *retPt++=i;
126           *retPt++=i;
127           ptI+=2;
128         }
129       else
130         {
131           pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
132           ptI[1]=ptI[0]+ci[1]-ci[0];
133           ptI++;
134           *retPt++=i;
135         }
136     }
137   _nodal_connec->decrRef();
138   _nodal_connec=newConn.retn();
139   _nodal_connec_index->decrRef();
140   _nodal_connec_index=newConnI.retn();
141   computeTypes();
142   updateTime();
143   return ret.retn();
144 }
145
146 /*!
147  * This method implements policy 1 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
148  */
149 DataArrayInt *MEDCouplingUMesh::simplexizePol1()
150 {
151   checkConnectivityFullyDefined();
152   if(getMeshDimension()!=2)
153     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePol0 : this policy is only available for mesh with meshdim == 2 !");
154   int nbOfCells=getNumberOfCells();
155   MCAuto<DataArrayInt> ret=DataArrayInt::New();
156   int nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
157   ret->alloc(nbOfCells+nbOfCutCells,1);
158   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
159   int *retPt=ret->getPointer();
160   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
161   MCAuto<DataArrayInt> newConnI=DataArrayInt::New();
162   newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
163   newConn->alloc(getNodalConnectivityArrayLen()+3*nbOfCutCells,1);
164   int *pt=newConn->getPointer();
165   int *ptI=newConnI->getPointer();
166   ptI[0]=0;
167   const int *oldc=_nodal_connec->begin();
168   const int *ci=_nodal_connec_index->begin();
169   for(int i=0;i<nbOfCells;i++,ci++)
170     {
171       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
172         {
173           const int tmp[8]={(int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+4],
174             (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+2],oldc[ci[0]+3],oldc[ci[0]+4]};
175           pt=std::copy(tmp,tmp+8,pt);
176           ptI[1]=ptI[0]+4;
177           ptI[2]=ptI[0]+8;
178           *retPt++=i;
179           *retPt++=i;
180           ptI+=2;
181         }
182       else
183         {
184           pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
185           ptI[1]=ptI[0]+ci[1]-ci[0];
186           ptI++;
187           *retPt++=i;
188         }
189     }
190   _nodal_connec->decrRef();
191   _nodal_connec=newConn.retn();
192   _nodal_connec_index->decrRef();
193   _nodal_connec_index=newConnI.retn();
194   computeTypes();
195   updateTime();
196   return ret.retn();
197 }
198
199 /*!
200  * This method implements policy INTERP_KERNEL::PLANAR_FACE_5 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
201  */
202 DataArrayInt *MEDCouplingUMesh::simplexizePlanarFace5()
203 {
204   checkConnectivityFullyDefined();
205   if(getMeshDimension()!=3)
206     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePlanarFace5 : this policy is only available for mesh with meshdim == 3 !");
207   int nbOfCells=getNumberOfCells();
208   MCAuto<DataArrayInt> ret=DataArrayInt::New();
209   int nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_HEXA8);
210   ret->alloc(nbOfCells+4*nbOfCutCells,1);
211   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
212   int *retPt=ret->getPointer();
213   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
214   MCAuto<DataArrayInt> newConnI=DataArrayInt::New();
215   newConnI->alloc(nbOfCells+4*nbOfCutCells+1,1);
216   newConn->alloc(getNodalConnectivityArrayLen()+16*nbOfCutCells,1);//21
217   int *pt=newConn->getPointer();
218   int *ptI=newConnI->getPointer();
219   ptI[0]=0;
220   const int *oldc=_nodal_connec->begin();
221   const int *ci=_nodal_connec_index->begin();
222   for(int i=0;i<nbOfCells;i++,ci++)
223     {
224       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_HEXA8)
225         {
226           for(int j=0;j<5;j++,pt+=5,ptI++)
227             {
228               pt[0]=(int)INTERP_KERNEL::NORM_TETRA4;
229               pt[1]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+0]+1]; pt[2]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+1]+1]; pt[3]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+2]+1]; pt[4]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_5_WO[4*j+3]+1];
230               *retPt++=i;
231               ptI[1]=ptI[0]+5;
232             }
233         }
234       else
235         {
236           pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
237           ptI[1]=ptI[0]+ci[1]-ci[0];
238           ptI++;
239           *retPt++=i;
240         }
241     }
242   _nodal_connec->decrRef();
243   _nodal_connec=newConn.retn();
244   _nodal_connec_index->decrRef();
245   _nodal_connec_index=newConnI.retn();
246   computeTypes();
247   updateTime();
248   return ret.retn();
249 }
250
251 /*!
252  * This method implements policy INTERP_KERNEL::PLANAR_FACE_6 of virtual method MEDCoupling::MEDCouplingUMesh::simplexize.
253  */
254 DataArrayInt *MEDCouplingUMesh::simplexizePlanarFace6()
255 {
256   checkConnectivityFullyDefined();
257   if(getMeshDimension()!=3)
258     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePlanarFace6 : this policy is only available for mesh with meshdim == 3 !");
259   int nbOfCells=getNumberOfCells();
260   MCAuto<DataArrayInt> ret=DataArrayInt::New();
261   int nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_HEXA8);
262   ret->alloc(nbOfCells+5*nbOfCutCells,1);
263   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
264   int *retPt=ret->getPointer();
265   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
266   MCAuto<DataArrayInt> newConnI=DataArrayInt::New();
267   newConnI->alloc(nbOfCells+5*nbOfCutCells+1,1);
268   newConn->alloc(getNodalConnectivityArrayLen()+21*nbOfCutCells,1);
269   int *pt=newConn->getPointer();
270   int *ptI=newConnI->getPointer();
271   ptI[0]=0;
272   const int *oldc=_nodal_connec->begin();
273   const int *ci=_nodal_connec_index->begin();
274   for(int i=0;i<nbOfCells;i++,ci++)
275     {
276       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_HEXA8)
277         {
278           for(int j=0;j<6;j++,pt+=5,ptI++)
279             {
280               pt[0]=(int)INTERP_KERNEL::NORM_TETRA4;
281               pt[1]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+0]+1]; pt[2]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+1]+1]; pt[3]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+2]+1]; pt[4]=oldc[ci[0]+INTERP_KERNEL::SPLIT_NODES_6_WO[4*j+3]+1];
282               *retPt++=i;
283               ptI[1]=ptI[0]+5;
284             }
285         }
286       else
287         {
288           pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
289           ptI[1]=ptI[0]+ci[1]-ci[0];
290           ptI++;
291           *retPt++=i;
292         }
293     }
294   _nodal_connec->decrRef();
295   _nodal_connec=newConn.retn();
296   _nodal_connec_index->decrRef();
297   _nodal_connec_index=newConnI.retn();
298   computeTypes();
299   updateTime();
300   return ret.retn();
301 }
302
303 /*!
304  * Tessellates \a this 2D mesh by dividing not straight edges of quadratic faces,
305  * so that the number of cells remains the same. Quadratic faces are converted to
306  * polygons. This method works only for 2D meshes in
307  * 2D space. If no cells are quadratic (INTERP_KERNEL::NORM_QUAD8,
308  * INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_QPOLYG ), \a this mesh remains unchanged.
309  * \warning This method can lead to a huge amount of nodes if \a eps is very low.
310  *  \param [in] eps - specifies the maximal angle (in radians) between 2 sub-edges of
311  *         a polylinized edge constituting the input polygon.
312  *  \throw If the coordinates array is not set.
313  *  \throw If the nodal connectivity of cells is not defined.
314  *  \throw If \a this->getMeshDimension() != 2.
315  *  \throw If \a this->getSpaceDimension() != 2.
316  */
317 void MEDCouplingUMesh::tessellate2DInternal(double eps)
318 {
319   checkFullyDefined();
320   if(getMeshDimension()!=2 || getSpaceDimension()!=2)
321     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DInternal works on umeshes with meshdim equal to 2 and spaceDim equal to 2 too!");
322   double epsa=fabs(eps);
323   if(epsa<std::numeric_limits<double>::min())
324     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DInternal : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !");
325   MCAuto<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
326   MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
327   revDesc1=0; revDescIndx1=0;
328   mDesc->tessellate2D(eps);
329   subDivide2DMesh(mDesc->_nodal_connec->begin(),mDesc->_nodal_connec_index->begin(),desc1->begin(),descIndx1->begin());
330   setCoords(mDesc->getCoords());
331 }
332
333 /*!
334  * Tessellates \a this 1D mesh in 2D space by dividing not straight quadratic edges.
335  * \warning This method can lead to a huge amount of nodes if \a eps is very low.
336  *  \param [in] eps - specifies the maximal angle (in radian) between 2 sub-edges of
337  *         a sub-divided edge.
338  *  \throw If the coordinates array is not set.
339  *  \throw If the nodal connectivity of cells is not defined.
340  *  \throw If \a this->getMeshDimension() != 1.
341  *  \throw If \a this->getSpaceDimension() != 2.
342  */
343 void MEDCouplingUMesh::tessellate2DCurveInternal(double eps)
344 {
345   checkFullyDefined();
346   if(getMeshDimension()!=1 || getSpaceDimension()!=2)
347     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurveInternal works on umeshes with meshdim equal to 1 and spaceDim equal to 2 too!");
348   double epsa=fabs(eps);
349   if(epsa<std::numeric_limits<double>::min())
350     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurveInternal : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !");
351   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=1.e-10;
352   int nbCells=getNumberOfCells();
353   int nbNodes=getNumberOfNodes();
354   const int *conn=_nodal_connec->begin();
355   const int *connI=_nodal_connec_index->begin();
356   const double *coords=_coords->begin();
357   std::vector<double> addCoo;
358   std::vector<int> newConn;//no direct DataArrayInt because interface with Geometric2D
359   MCAuto<DataArrayInt> newConnI(DataArrayInt::New());
360   newConnI->alloc(nbCells+1,1);
361   int *newConnIPtr=newConnI->getPointer();
362   *newConnIPtr=0;
363   int tmp1[3];
364   INTERP_KERNEL::Node *tmp2[3];
365   std::set<INTERP_KERNEL::NormalizedCellType> types;
366   for(int i=0;i<nbCells;i++,newConnIPtr++)
367     {
368       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
369       if(cm.isQuadratic())
370         {//assert(connI[i+1]-connI[i]-1==3)
371           tmp1[0]=conn[connI[i]+1+0]; tmp1[1]=conn[connI[i]+1+1]; tmp1[2]=conn[connI[i]+1+2];
372           tmp2[0]=new INTERP_KERNEL::Node(coords[2*tmp1[0]],coords[2*tmp1[0]+1]);
373           tmp2[1]=new INTERP_KERNEL::Node(coords[2*tmp1[1]],coords[2*tmp1[1]+1]);
374           tmp2[2]=new INTERP_KERNEL::Node(coords[2*tmp1[2]],coords[2*tmp1[2]+1]);
375           INTERP_KERNEL::EdgeArcCircle *eac=INTERP_KERNEL::EdgeArcCircle::BuildFromNodes(tmp2[0],tmp2[2],tmp2[1]);
376           if(eac)
377             {
378               eac->tesselate(tmp1,nbNodes,epsa,newConn,addCoo);
379               types.insert((INTERP_KERNEL::NormalizedCellType)newConn[newConnIPtr[0]]);
380               delete eac;
381               newConnIPtr[1]=(int)newConn.size();
382             }
383           else
384             {
385               types.insert(INTERP_KERNEL::NORM_SEG2);
386               newConn.push_back(INTERP_KERNEL::NORM_SEG2);
387               newConn.insert(newConn.end(),conn+connI[i]+1,conn+connI[i]+3);
388               newConnIPtr[1]=newConnIPtr[0]+3;
389             }
390         }
391       else
392         {
393           types.insert((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
394           newConn.insert(newConn.end(),conn+connI[i],conn+connI[i+1]);
395           newConnIPtr[1]=newConnIPtr[0]+3;
396         }
397     }
398   if(addCoo.empty() && newConn.size()==_nodal_connec->getNumberOfTuples())//nothing happens during tessellation : no update needed
399     return ;
400   _types=types;
401   DataArrayInt::SetArrayIn(newConnI,_nodal_connec_index);
402   MCAuto<DataArrayInt> newConnArr=DataArrayInt::New();
403   newConnArr->alloc((int)newConn.size(),1);
404   std::copy(newConn.begin(),newConn.end(),newConnArr->getPointer());
405   DataArrayInt::SetArrayIn(newConnArr,_nodal_connec);
406   MCAuto<DataArrayDouble> newCoords=DataArrayDouble::New();
407   newCoords->alloc(nbNodes+((int)addCoo.size())/2,2);
408   double *work=std::copy(_coords->begin(),_coords->end(),newCoords->getPointer());
409   std::copy(addCoo.begin(),addCoo.end(),work);
410   DataArrayDouble::SetArrayIn(newCoords,_coords);
411   updateTime();
412 }
413
414
415 /*!
416  * This private method is used to subdivide edges of a mesh with meshdim==2. If \a this has no a meshdim equal to 2 an exception will be thrown.
417  * This method completly ignore coordinates.
418  * \param nodeSubdived is the nodal connectivity of subdivision of edges
419  * \param nodeIndxSubdived is the nodal connectivity index of subdivision of edges
420  * \param desc is descending connectivity in format specified in MEDCouplingUMesh::buildDescendingConnectivity2
421  * \param descIndex is descending connectivity index in format specified in MEDCouplingUMesh::buildDescendingConnectivity2
422  */
423 void MEDCouplingUMesh::subDivide2DMesh(const int *nodeSubdived, const int *nodeIndxSubdived, const int *desc, const int *descIndex)
424 {
425   checkFullyDefined();
426   if(getMeshDimension()!=2)
427     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : works only on umesh with meshdim==2 !");
428   int nbOfCells=getNumberOfCells();
429   int *connI=_nodal_connec_index->getPointer();
430   int newConnLgth=0;
431   for(int i=0;i<nbOfCells;i++,connI++)
432     {
433       int offset=descIndex[i];
434       int nbOfEdges=descIndex[i+1]-offset;
435       //
436       bool ddirect=desc[offset+nbOfEdges-1]>0;
437       int eedgeId=std::abs(desc[offset+nbOfEdges-1])-1;
438       int ref=ddirect?nodeSubdived[nodeIndxSubdived[eedgeId+1]-1]:nodeSubdived[nodeIndxSubdived[eedgeId]+1];
439       for(int j=0;j<nbOfEdges;j++)
440         {
441           bool direct=desc[offset+j]>0;
442           int edgeId=std::abs(desc[offset+j])-1;
443           if(!INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodeSubdived[nodeIndxSubdived[edgeId]]).isQuadratic())
444             {
445               int id1=nodeSubdived[nodeIndxSubdived[edgeId]+1];
446               int id2=nodeSubdived[nodeIndxSubdived[edgeId+1]-1];
447               int ref2=direct?id1:id2;
448               if(ref==ref2)
449                 {
450                   int nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1;
451                   newConnLgth+=nbOfSubNodes-1;
452                   ref=direct?id2:id1;
453                 }
454               else
455                 {
456                   std::ostringstream oss; oss << "MEDCouplingUMesh::subDivide2DMesh : On polygon #" << i << " edgeid #" << j << " subedges mismatch : end subedge k!=start subedge k+1 !";
457                   throw INTERP_KERNEL::Exception(oss.str());
458                 }
459             }
460           else
461             {
462               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : this method only subdivides into linear edges !");
463             }
464         }
465       newConnLgth++;//+1 is for cell type
466       connI[1]=newConnLgth;
467     }
468   //
469   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
470   newConn->alloc(newConnLgth,1);
471   int *work=newConn->getPointer();
472   for(int i=0;i<nbOfCells;i++)
473     {
474       *work++=INTERP_KERNEL::NORM_POLYGON;
475       int offset=descIndex[i];
476       int nbOfEdges=descIndex[i+1]-offset;
477       for(int j=0;j<nbOfEdges;j++)
478         {
479           bool direct=desc[offset+j]>0;
480           int edgeId=std::abs(desc[offset+j])-1;
481           if(direct)
482             work=std::copy(nodeSubdived+nodeIndxSubdived[edgeId]+1,nodeSubdived+nodeIndxSubdived[edgeId+1]-1,work);
483           else
484             {
485               int nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1;
486               std::reverse_iterator<const int *> it(nodeSubdived+nodeIndxSubdived[edgeId+1]);
487               work=std::copy(it,it+nbOfSubNodes-1,work);
488             }
489         }
490     }
491   DataArrayInt::SetArrayIn(newConn,_nodal_connec);
492   _types.clear();
493   if(nbOfCells>0)
494     _types.insert(INTERP_KERNEL::NORM_POLYGON);
495 }
496
497 /*!
498  * Keeps from \a this only cells which constituing point id are in the ids specified by [ \a begin,\a end ).
499  * The resulting cell ids are stored at the end of the 'cellIdsKept' parameter.
500  * Parameter \a fullyIn specifies if a cell that has part of its nodes in ids array is kept or not.
501  * If \a fullyIn is true only cells whose ids are \b fully contained in [ \a begin,\a end ) tab will be kept.
502  *
503  * \param [in] begin input start of array of node ids.
504  * \param [in] end input end of array of node ids.
505  * \param [in] fullyIn input that specifies if all node ids must be in [ \a begin,\a end ) array to consider cell to be in.
506  * \param [in,out] cellIdsKeptArr array where all candidate cell ids are put at the end.
507  */
508 void MEDCouplingUMesh::fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const
509 {
510   MCAuto<DataArrayInt> cellIdsKept=DataArrayInt::New(); cellIdsKept->alloc(0,1);
511   checkConnectivityFullyDefined();
512   int tmp=-1;
513   int sz=getNodalConnectivity()->getMaxValue(tmp); sz=std::max(sz,0)+1;
514   std::vector<bool> fastFinder(sz,false);
515   for(const int *work=begin;work!=end;work++)
516     if(*work>=0 && *work<sz)
517       fastFinder[*work]=true;
518   int nbOfCells=getNumberOfCells();
519   const int *conn=getNodalConnectivity()->getConstPointer();
520   const int *connIndex=getNodalConnectivityIndex()->getConstPointer();
521   for(int i=0;i<nbOfCells;i++)
522     {
523       int ref=0,nbOfHit=0;
524       for(const int *work2=conn+connIndex[i]+1;work2!=conn+connIndex[i+1];work2++)
525         if(*work2>=0)
526           {
527             ref++;
528             if(fastFinder[*work2])
529               nbOfHit++;
530           }
531       if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
532         cellIdsKept->pushBackSilent(i);
533     }
534   cellIdsKeptArr=cellIdsKept.retn();
535 }
536
537 /*!
538  * This method works on a 3D curve linear mesh that is to say (meshDim==1 and spaceDim==3).
539  * If it is not the case an exception will be thrown.
540  * This method is non const because the coordinate of \a this can be appended with some new points issued from
541  * intersection of plane defined by ('origin','vec').
542  * This method has one in/out parameter : 'cut3DCurve'.
543  * Param 'cut3DCurve' is expected to be of size 'this->getNumberOfCells()'. For each i in [0,'this->getNumberOfCells()')
544  * if cut3DCurve[i]==-2, it means that for cell #i in \a this nothing has been detected previously.
545  * if cut3DCurve[i]==-1, it means that cell#i has been already detected to be fully part of plane defined by ('origin','vec').
546  * This method will throw an exception if \a this contains a non linear segment.
547  */
548 void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double *vec, double eps, std::vector<int>& cut3DCurve)
549 {
550   checkFullyDefined();
551   if(getMeshDimension()!=1 || getSpaceDimension()!=3)
552     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane works on umeshes with meshdim equal to 1 and spaceDim equal to 3 !");
553   int ncells=getNumberOfCells();
554   int nnodes=getNumberOfNodes();
555   double vec2[3],vec3[3],vec4[3];
556   double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
557   if(normm<1e-6)
558     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : parameter 'vec' should have a norm2 greater than 1e-6 !");
559   vec2[0]=vec[0]/normm; vec2[1]=vec[1]/normm; vec2[2]=vec[2]/normm;
560   const int *conn=_nodal_connec->getConstPointer();
561   const int *connI=_nodal_connec_index->getConstPointer();
562   const double *coo=_coords->getConstPointer();
563   std::vector<double> addCoo;
564   for(int i=0;i<ncells;i++)
565     {
566       if(conn[connI[i]]==(int)INTERP_KERNEL::NORM_SEG2)
567         {
568           if(cut3DCurve[i]==-2)
569             {
570               int st=conn[connI[i]+1],endd=conn[connI[i]+2];
571               vec3[0]=coo[3*endd]-coo[3*st]; vec3[1]=coo[3*endd+1]-coo[3*st+1]; vec3[2]=coo[3*endd+2]-coo[3*st+2];
572               double normm2=sqrt(vec3[0]*vec3[0]+vec3[1]*vec3[1]+vec3[2]*vec3[2]);
573               double colin=std::abs((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2])/normm2);
574               if(colin>eps)//if colin<=eps -> current SEG2 is colinear to the input plane
575                 {
576                   const double *st2=coo+3*st;
577                   vec4[0]=st2[0]-origin[0]; vec4[1]=st2[1]-origin[1]; vec4[2]=st2[2]-origin[2];
578                   double pos=-(vec4[0]*vec2[0]+vec4[1]*vec2[1]+vec4[2]*vec2[2])/((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2]));
579                   if(pos>eps && pos<1-eps)
580                     {
581                       int nNode=((int)addCoo.size())/3;
582                       vec4[0]=st2[0]+pos*vec3[0]; vec4[1]=st2[1]+pos*vec3[1]; vec4[2]=st2[2]+pos*vec3[2];
583                       addCoo.insert(addCoo.end(),vec4,vec4+3);
584                       cut3DCurve[i]=nnodes+nNode;
585                     }
586                 }
587             }
588         }
589       else
590         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : this method is only available for linear cell (NORM_SEG2) !");
591     }
592   if(!addCoo.empty())
593     {
594       int newNbOfNodes=nnodes+((int)addCoo.size())/3;
595       MCAuto<DataArrayDouble> coo2=DataArrayDouble::New();
596       coo2->alloc(newNbOfNodes,3);
597       double *tmp=coo2->getPointer();
598       tmp=std::copy(_coords->begin(),_coords->end(),tmp);
599       std::copy(addCoo.begin(),addCoo.end(),tmp);
600       DataArrayDouble::SetArrayIn(coo2,_coords);
601     }
602 }
603
604 /*!
605  * This method incarnates the policy 0 for MEDCouplingUMesh::buildExtrudedMesh method.
606  * \param mesh1D is the input 1D mesh used for translation computation.
607  * \return newCoords new coords filled by this method.
608  */
609 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
610 {
611   int oldNbOfNodes=getNumberOfNodes();
612   int nbOf1DCells=mesh1D->getNumberOfCells();
613   int spaceDim=getSpaceDimension();
614   DataArrayDouble *ret=DataArrayDouble::New();
615   std::vector<bool> isQuads;
616   int nbOfLevsInVec=isQuad?2*nbOf1DCells+1:nbOf1DCells+1;
617   ret->alloc(oldNbOfNodes*nbOfLevsInVec,spaceDim);
618   double *retPtr=ret->getPointer();
619   const double *coords=getCoords()->getConstPointer();
620   double *work=std::copy(coords,coords+spaceDim*oldNbOfNodes,retPtr);
621   std::vector<int> v;
622   std::vector<double> c;
623   double vec[3];
624   v.reserve(3);
625   c.reserve(6);
626   for(int i=0;i<nbOf1DCells;i++)
627     {
628       v.resize(0);
629       mesh1D->getNodeIdsOfCell(i,v);
630       c.resize(0);
631       mesh1D->getCoordinatesOfNode(v[isQuad?2:1],c);
632       mesh1D->getCoordinatesOfNode(v[0],c);
633       std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
634       for(int j=0;j<oldNbOfNodes;j++)
635         work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
636       if(isQuad)
637         {
638           c.resize(0);
639           mesh1D->getCoordinatesOfNode(v[1],c);
640           mesh1D->getCoordinatesOfNode(v[0],c);
641           std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
642           for(int j=0;j<oldNbOfNodes;j++)
643             work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
644         }
645     }
646   ret->copyStringInfoFrom(*getCoords());
647   return ret;
648 }
649
650 /*!
651  * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
652  * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
653  * \return newCoords new coords filled by this method.
654  */
655 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
656 {
657   if(mesh1D->getSpaceDimension()==2)
658     return fillExtCoordsUsingTranslAndAutoRotation2D(mesh1D,isQuad);
659   if(mesh1D->getSpaceDimension()==3)
660     return fillExtCoordsUsingTranslAndAutoRotation3D(mesh1D,isQuad);
661   throw INTERP_KERNEL::Exception("Not implemented rotation and translation alg. for spacedim other than 2 and 3 !");
662 }
663
664 /*!
665  * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
666  * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
667  * \return newCoords new coords filled by this method.
668  */
669 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
670 {
671   if(isQuad)
672     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : not implemented for quadratic cells !");
673   int oldNbOfNodes=getNumberOfNodes();
674   int nbOf1DCells=mesh1D->getNumberOfCells();
675   if(nbOf1DCells<2)
676     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
677   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
678   int nbOfLevsInVec=nbOf1DCells+1;
679   ret->alloc(oldNbOfNodes*nbOfLevsInVec,2);
680   double *retPtr=ret->getPointer();
681   retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
682   MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
683   MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
684   tmp->setCoords(tmp2);
685   const double *coo1D=mesh1D->getCoords()->getConstPointer();
686   const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
687   const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
688   for(int i=1;i<nbOfLevsInVec;i++)
689     {
690       const double *begin=coo1D+2*conn1D[connI1D[i-1]+1];
691       const double *end=coo1D+2*conn1D[connI1D[i-1]+2];
692       const double *third=i+1<nbOfLevsInVec?coo1D+2*conn1D[connI1D[i]+2]:coo1D+2*conn1D[connI1D[i-2]+1];
693       const double vec[2]={end[0]-begin[0],end[1]-begin[1]};
694       tmp->translate(vec);
695       double tmp3[2],radius,alpha,alpha0;
696       const double *p0=i+1<nbOfLevsInVec?begin:third;
697       const double *p1=i+1<nbOfLevsInVec?end:begin;
698       const double *p2=i+1<nbOfLevsInVec?third:end;
699       INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
700       double cosangle=i+1<nbOfLevsInVec?(p0[0]-tmp3[0])*(p1[0]-tmp3[0])+(p0[1]-tmp3[1])*(p1[1]-tmp3[1]):(p2[0]-tmp3[0])*(p1[0]-tmp3[0])+(p2[1]-tmp3[1])*(p1[1]-tmp3[1]);
701       double angle=acos(cosangle/(radius*radius));
702       tmp->rotate(end,0,angle);
703       retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
704     }
705   return ret.retn();
706 }
707
708 /*!
709  * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
710  * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
711  * \return newCoords new coords filled by this method.
712  */
713 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
714 {
715   if(isQuad)
716     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : not implemented for quadratic cells !");
717   int oldNbOfNodes=getNumberOfNodes();
718   int nbOf1DCells=mesh1D->getNumberOfCells();
719   if(nbOf1DCells<2)
720     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
721   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
722   int nbOfLevsInVec=nbOf1DCells+1;
723   ret->alloc(oldNbOfNodes*nbOfLevsInVec,3);
724   double *retPtr=ret->getPointer();
725   retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
726   MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
727   MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
728   tmp->setCoords(tmp2);
729   const double *coo1D=mesh1D->getCoords()->getConstPointer();
730   const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
731   const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
732   for(int i=1;i<nbOfLevsInVec;i++)
733     {
734       const double *begin=coo1D+3*conn1D[connI1D[i-1]+1];
735       const double *end=coo1D+3*conn1D[connI1D[i-1]+2];
736       const double *third=i+1<nbOfLevsInVec?coo1D+3*conn1D[connI1D[i]+2]:coo1D+3*conn1D[connI1D[i-2]+1];
737       const double vec[3]={end[0]-begin[0],end[1]-begin[1],end[2]-begin[2]};
738       tmp->translate(vec);
739       double tmp3[2],radius,alpha,alpha0;
740       const double *p0=i+1<nbOfLevsInVec?begin:third;
741       const double *p1=i+1<nbOfLevsInVec?end:begin;
742       const double *p2=i+1<nbOfLevsInVec?third:end;
743       double vecPlane[3]={
744         (p1[1]-p0[1])*(p2[2]-p1[2])-(p1[2]-p0[2])*(p2[1]-p1[1]),
745         (p1[2]-p0[2])*(p2[0]-p1[0])-(p1[0]-p0[0])*(p2[2]-p1[2]),
746         (p1[0]-p0[0])*(p2[1]-p1[1])-(p1[1]-p0[1])*(p2[0]-p1[0]),
747       };
748       double norm=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]+vecPlane[2]*vecPlane[2]);
749       if(norm>1.e-7)
750         {
751           vecPlane[0]/=norm; vecPlane[1]/=norm; vecPlane[2]/=norm;
752           double norm2=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]);
753           double vec2[2]={vecPlane[1]/norm2,-vecPlane[0]/norm2};
754           double s2=norm2;
755           double c2=cos(asin(s2));
756           double m[3][3]={
757             {vec2[0]*vec2[0]*(1-c2)+c2, vec2[0]*vec2[1]*(1-c2), vec2[1]*s2},
758             {vec2[0]*vec2[1]*(1-c2), vec2[1]*vec2[1]*(1-c2)+c2, -vec2[0]*s2},
759             {-vec2[1]*s2, vec2[0]*s2, c2}
760           };
761           double p0r[3]={m[0][0]*p0[0]+m[0][1]*p0[1]+m[0][2]*p0[2], m[1][0]*p0[0]+m[1][1]*p0[1]+m[1][2]*p0[2], m[2][0]*p0[0]+m[2][1]*p0[1]+m[2][2]*p0[2]};
762           double p1r[3]={m[0][0]*p1[0]+m[0][1]*p1[1]+m[0][2]*p1[2], m[1][0]*p1[0]+m[1][1]*p1[1]+m[1][2]*p1[2], m[2][0]*p1[0]+m[2][1]*p1[1]+m[2][2]*p1[2]};
763           double p2r[3]={m[0][0]*p2[0]+m[0][1]*p2[1]+m[0][2]*p2[2], m[1][0]*p2[0]+m[1][1]*p2[1]+m[1][2]*p2[2], m[2][0]*p2[0]+m[2][1]*p2[1]+m[2][2]*p2[2]};
764           INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0);
765           double cosangle=i+1<nbOfLevsInVec?(p0r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p0r[1]-tmp3[1])*(p1r[1]-tmp3[1]):(p2r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p2r[1]-tmp3[1])*(p1r[1]-tmp3[1]);
766           double angle=acos(cosangle/(radius*radius));
767           tmp->rotate(end,vecPlane,angle);
768         }
769       retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
770     }
771   return ret.retn();
772 }
773
774 /*!
775  * This method is private because not easy to use for end user. This method is const contrary to
776  * MEDCouplingUMesh::buildExtrudedMesh method because this->_coords are expected to contain
777  * the coords sorted slice by slice.
778  * \param isQuad specifies presence of quadratic cells.
779  */
780 MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(int nbOfNodesOf1Lev, bool isQuad) const
781 {
782   int nbOf1DCells(getNumberOfNodes()/nbOfNodesOf1Lev-1);
783   int nbOf2DCells(getNumberOfCells());
784   int nbOf3DCells(nbOf2DCells*nbOf1DCells);
785   MEDCouplingUMesh *ret(MEDCouplingUMesh::New("Extruded",getMeshDimension()+1));
786   const int *conn(_nodal_connec->begin()),*connI(_nodal_connec_index->begin());
787   MCAuto<DataArrayInt> newConn(DataArrayInt::New()),newConnI(DataArrayInt::New());
788   newConnI->alloc(nbOf3DCells+1,1);
789   int *newConnIPtr(newConnI->getPointer());
790   *newConnIPtr++=0;
791   std::vector<int> newc;
792   for(int j=0;j<nbOf2DCells;j++)
793     {
794       AppendExtrudedCell(conn+connI[j],conn+connI[j+1],nbOfNodesOf1Lev,isQuad,newc);
795       *newConnIPtr++=(int)newc.size();
796     }
797   newConn->alloc((int)(newc.size())*nbOf1DCells,1);
798   int *newConnPtr(newConn->getPointer());
799   int deltaPerLev(isQuad?2*nbOfNodesOf1Lev:nbOfNodesOf1Lev);
800   newConnIPtr=newConnI->getPointer();
801   for(int iz=0;iz<nbOf1DCells;iz++)
802     {
803       if(iz!=0)
804         std::transform(newConnIPtr+1,newConnIPtr+1+nbOf2DCells,newConnIPtr+1+iz*nbOf2DCells,std::bind2nd(std::plus<int>(),newConnIPtr[iz*nbOf2DCells]));
805       const int *posOfTypeOfCell(newConnIPtr);
806       for(std::vector<int>::const_iterator iter=newc.begin();iter!=newc.end();iter++,newConnPtr++)
807         {
808           int icell((int)(iter-newc.begin()));//std::distance unfortunately cannot been called here in C++98
809           if(icell!=*posOfTypeOfCell)
810             {
811               if(*iter!=-1)
812                 *newConnPtr=(*iter)+iz*deltaPerLev;
813               else
814                 *newConnPtr=-1;
815             }
816           else
817             {
818               *newConnPtr=*iter;
819               posOfTypeOfCell++;
820             }
821         }
822     }
823   ret->setConnectivity(newConn,newConnI,true);
824   ret->setCoords(getCoords());
825   return ret;
826 }
827
828
829 /*!
830  * This method find in candidate pool defined by 'candidates' the cells equal following the polycy 'compType'.
831  * If any true is returned and the results will be put at the end of 'result' output parameter. If not false is returned
832  * and result remains unchanged.
833  * The semantic of 'compType' is specified in MEDCouplingPointSet::zipConnectivityTraducer method.
834  * If in 'candidates' pool -1 value is considered as an empty value.
835  * WARNING this method returns only ONE set of result !
836  */
837 bool MEDCouplingUMesh::AreCellsEqualInPool(const std::vector<int>& candidates, int compType, const int *conn, const int *connI, DataArrayInt *result)
838 {
839   if(candidates.size()<1)
840     return false;
841   bool ret=false;
842   std::vector<int>::const_iterator iter=candidates.begin();
843   int start=(*iter++);
844   for(;iter!=candidates.end();iter++)
845     {
846       int status=AreCellsEqual(conn,connI,start,*iter,compType);
847       if(status!=0)
848         {
849           if(!ret)
850             {
851               result->pushBackSilent(start);
852               ret=true;
853             }
854           if(status==1)
855             result->pushBackSilent(*iter);
856           else
857             result->pushBackSilent(status==2?(*iter+1):-(*iter+1));
858         }
859     }
860   return ret;
861 }
862
863 /*!
864  * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelf.
865  * Keeps from \a this only cells which constituing point id are in the ids specified by [ \a begin,\a end ).
866  * The return newly allocated mesh will share the same coordinates as \a this.
867  */
868 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *begin, const int *end) const
869 {
870   checkConnectivityFullyDefined();
871   int ncell=getNumberOfCells();
872   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
873   ret->_mesh_dim=_mesh_dim;
874   ret->setCoords(_coords);
875   std::size_t nbOfElemsRet=std::distance(begin,end);
876   int *connIndexRet=(int *)malloc((nbOfElemsRet+1)*sizeof(int));
877   connIndexRet[0]=0;
878   const int *conn=_nodal_connec->getConstPointer();
879   const int *connIndex=_nodal_connec_index->getConstPointer();
880   int newNbring=0;
881   for(const int *work=begin;work!=end;work++,newNbring++)
882     {
883       if(*work>=0 && *work<ncell)
884         connIndexRet[newNbring+1]=connIndexRet[newNbring]+connIndex[*work+1]-connIndex[*work];
885       else
886         {
887           free(connIndexRet);
888           std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoords : On pos #" << std::distance(begin,work) << " input cell id =" << *work << " should be in [0," << ncell << ") !";
889           throw INTERP_KERNEL::Exception(oss.str());
890         }
891     }
892   int *connRet=(int *)malloc(connIndexRet[nbOfElemsRet]*sizeof(int));
893   int *connRetWork=connRet;
894   std::set<INTERP_KERNEL::NormalizedCellType> types;
895   for(const int *work=begin;work!=end;work++)
896     {
897       types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*work]]);
898       connRetWork=std::copy(conn+connIndex[*work],conn+connIndex[*work+1],connRetWork);
899     }
900   MCAuto<DataArrayInt> connRetArr=DataArrayInt::New();
901   connRetArr->useArray(connRet,true,C_DEALLOC,connIndexRet[nbOfElemsRet],1);
902   MCAuto<DataArrayInt> connIndexRetArr=DataArrayInt::New();
903   connIndexRetArr->useArray(connIndexRet,true,C_DEALLOC,(int)nbOfElemsRet+1,1);
904   ret->setConnectivity(connRetArr,connIndexRetArr,false);
905   ret->_types=types;
906   ret->copyTinyInfoFrom(this);
907   return ret.retn();
908 }
909
910 /*!
911  * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelfSlice.
912  * CellIds are given using range specified by a start an end and step.
913  */
914 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice(int start, int end, int step) const
915 {
916   checkFullyDefined();
917   int ncell=getNumberOfCells();
918   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
919   ret->_mesh_dim=_mesh_dim;
920   ret->setCoords(_coords);
921   int newNbOfCells=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : ");
922   MCAuto<DataArrayInt> newConnI=DataArrayInt::New(); newConnI->alloc(newNbOfCells+1,1);
923   int *newConnIPtr=newConnI->getPointer(); *newConnIPtr=0;
924   int work=start;
925   const int *conn=_nodal_connec->getConstPointer();
926   const int *connIndex=_nodal_connec_index->getConstPointer();
927   for(int i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
928     {
929       if(work>=0 && work<ncell)
930         {
931           newConnIPtr[1]=newConnIPtr[0]+connIndex[work+1]-connIndex[work];
932         }
933       else
934         {
935           std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : On pos #" << i << " input cell id =" << work << " should be in [0," << ncell << ") !";
936           throw INTERP_KERNEL::Exception(oss.str());
937         }
938     }
939   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(newConnIPtr[0],1);
940   int *newConnPtr=newConn->getPointer();
941   std::set<INTERP_KERNEL::NormalizedCellType> types;
942   work=start;
943   for(int i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
944     {
945       types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[work]]);
946       newConnPtr=std::copy(conn+connIndex[work],conn+connIndex[work+1],newConnPtr);
947     }
948   ret->setConnectivity(newConn,newConnI,false);
949   ret->_types=types;
950   ret->copyTinyInfoFrom(this);
951   return ret.retn();
952 }
953
954
955 int MEDCouplingFastNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2)
956 {
957   return id;
958 }
959
960 int MEDCouplingOrientationSensitiveNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2)
961 {
962   if(!compute)
963     return id+1;
964   else
965     {
966       if(cm.getOrientationStatus(nb,conn1,conn2))
967         return id+1;
968       else
969         return -(id+1);
970     }
971 }
972
973
974 /*!
975  * Implementes \a conversionType 0 for meshes with meshDim = 1, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
976  * \return a newly created DataArrayInt instance that the caller should deal with containing cell ids of converted cells.
977  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
978  */
979 DataArrayInt *MEDCouplingUMesh::convertLinearCellsToQuadratic1D0(DataArrayInt *&conn, DataArrayInt *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
980 {
981   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
982   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(0,1);
983   MCAuto<DataArrayInt> newConnI=DataArrayInt::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
984   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(0,1);
985   int nbOfCells=getNumberOfCells();
986   int nbOfNodes=getNumberOfNodes();
987   const int *cPtr=_nodal_connec->begin();
988   const int *icPtr=_nodal_connec_index->begin();
989   int lastVal=0,offset=nbOfNodes;
990   for(int i=0;i<nbOfCells;i++,icPtr++)
991     {
992       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
993       if(type==INTERP_KERNEL::NORM_SEG2)
994         {
995           types.insert(INTERP_KERNEL::NORM_SEG3);
996           newConn->pushBackSilent((int)INTERP_KERNEL::NORM_SEG3);
997           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[0]+3);
998           newConn->pushBackSilent(offset++);
999           lastVal+=4;
1000           newConnI->pushBackSilent(lastVal);
1001           ret->pushBackSilent(i);
1002         }
1003       else
1004         {
1005           types.insert(type);
1006           lastVal+=(icPtr[1]-icPtr[0]);
1007           newConnI->pushBackSilent(lastVal);
1008           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1009         }
1010     }
1011   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1012   coords=DataArrayDouble::Aggregate(getCoords(),tmp); conn=newConn.retn(); connI=newConnI.retn();
1013   return ret.retn();
1014 }
1015
1016 DataArrayInt *MEDCouplingUMesh::convertLinearCellsToQuadratic2DAnd3D0(const MEDCouplingUMesh *m1D, const DataArrayInt *desc, const DataArrayInt *descI, DataArrayInt *&conn, DataArrayInt *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1017 {
1018   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(0,1);
1019   MCAuto<DataArrayInt> newConnI=DataArrayInt::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1020   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(0,1);
1021   //
1022   const int *descPtr(desc->begin()),*descIPtr(descI->begin());
1023   DataArrayInt *conn1D=0,*conn1DI=0;
1024   std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1025   DataArrayDouble *coordsTmp=0;
1026   MCAuto<DataArrayInt> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1027   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1028   MCAuto<DataArrayInt> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1029   const int *c1DPtr=conn1D->begin();
1030   const int *c1DIPtr=conn1DI->begin();
1031   int nbOfCells=getNumberOfCells();
1032   const int *cPtr=_nodal_connec->begin();
1033   const int *icPtr=_nodal_connec_index->begin();
1034   int lastVal=0;
1035   for(int i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1036     {
1037       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1038       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1039       if(!cm.isQuadratic())
1040         {
1041           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType();
1042           types.insert(typ2); newConn->pushBackSilent(typ2);
1043           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1044           for(const int *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1045             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1046           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0]);
1047           newConnI->pushBackSilent(lastVal);
1048           ret->pushBackSilent(i);
1049         }
1050       else
1051         {
1052           types.insert(typ);
1053           lastVal+=(icPtr[1]-icPtr[0]);
1054           newConnI->pushBackSilent(lastVal);
1055           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1056         }
1057     }
1058   conn=newConn.retn(); connI=newConnI.retn(); coords=coordsTmpSafe.retn();
1059   return ret.retn();
1060 }
1061
1062 /*!
1063  * Implementes \a conversionType 0 for meshes with meshDim = 2, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1064  * \return a newly created DataArrayInt instance that the caller should deal with containing cell ids of converted cells.
1065  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1066  */
1067 DataArrayInt *MEDCouplingUMesh::convertLinearCellsToQuadratic2D0(DataArrayInt *&conn, DataArrayInt *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1068 {
1069   MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),tmp2(DataArrayInt::New()),tmp3(DataArrayInt::New());
1070   MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1071   return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1072 }
1073
1074 DataArrayInt *MEDCouplingUMesh::convertLinearCellsToQuadratic2D1(DataArrayInt *&conn, DataArrayInt *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1075 {
1076   MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),tmp2(DataArrayInt::New()),tmp3(DataArrayInt::New());
1077   MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1078   //
1079   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(0,1);
1080   MCAuto<DataArrayInt> newConnI=DataArrayInt::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1081   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(0,1);
1082   //
1083   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1084   const int *descPtr(desc->begin()),*descIPtr(descI->begin());
1085   DataArrayInt *conn1D=0,*conn1DI=0;
1086   std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1087   DataArrayDouble *coordsTmp=0;
1088   MCAuto<DataArrayInt> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1089   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1090   MCAuto<DataArrayInt> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1091   const int *c1DPtr=conn1D->begin();
1092   const int *c1DIPtr=conn1DI->begin();
1093   int nbOfCells=getNumberOfCells();
1094   const int *cPtr=_nodal_connec->begin();
1095   const int *icPtr=_nodal_connec_index->begin();
1096   int lastVal=0,offset=coordsTmpSafe->getNumberOfTuples();
1097   for(int i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1098     {
1099       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1100       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1101       if(!cm.isQuadratic())
1102         {
1103           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1104           types.insert(typ2); newConn->pushBackSilent(typ2);
1105           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1106           for(const int *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1107             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1108           newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1109           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+1;
1110           newConnI->pushBackSilent(lastVal);
1111           ret->pushBackSilent(i);
1112         }
1113       else
1114         {
1115           types.insert(typ);
1116           lastVal+=(icPtr[1]-icPtr[0]);
1117           newConnI->pushBackSilent(lastVal);
1118           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1119         }
1120     }
1121   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1122   coords=DataArrayDouble::Aggregate(coordsTmpSafe,tmp); conn=newConn.retn(); connI=newConnI.retn();
1123   return ret.retn();
1124 }
1125
1126 /*!
1127  * Implementes \a conversionType 0 for meshes with meshDim = 3, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1128  * \return a newly created DataArrayInt instance that the caller should deal with containing cell ids of converted cells.
1129  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1130  */
1131 DataArrayInt *MEDCouplingUMesh::convertLinearCellsToQuadratic3D0(DataArrayInt *&conn, DataArrayInt *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1132 {
1133   MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),tmp2(DataArrayInt::New()),tmp3(DataArrayInt::New());
1134   MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1135   return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1136 }
1137
1138 DataArrayInt *MEDCouplingUMesh::convertLinearCellsToQuadratic3D1(DataArrayInt *&conn, DataArrayInt *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1139 {
1140   MCAuto<DataArrayInt> desc2(DataArrayInt::New()),desc2I(DataArrayInt::New()),tmp2(DataArrayInt::New()),tmp3(DataArrayInt::New());
1141   MCAuto<MEDCouplingUMesh> m2D=buildDescendingConnectivityGen<MinusOneSonsGeneratorBiQuadratic>(desc2,desc2I,tmp2,tmp3,MEDCouplingFastNbrer); tmp2=0; tmp3=0;
1142   MCAuto<DataArrayInt> desc1(DataArrayInt::New()),desc1I(DataArrayInt::New()),tmp4(DataArrayInt::New()),tmp5(DataArrayInt::New());
1143   MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc1,desc1I,tmp4,tmp5); tmp4=0; tmp5=0;
1144   //
1145   MCAuto<DataArrayInt> newConn=DataArrayInt::New(); newConn->alloc(0,1);
1146   MCAuto<DataArrayInt> newConnI=DataArrayInt::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1147   MCAuto<DataArrayInt> ret=DataArrayInt::New(),ret2=DataArrayInt::New(); ret->alloc(0,1); ret2->alloc(0,1);
1148   //
1149   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1150   const int *descPtr(desc1->begin()),*descIPtr(desc1I->begin()),*desc2Ptr(desc2->begin()),*desc2IPtr(desc2I->begin());
1151   DataArrayInt *conn1D=0,*conn1DI=0,*conn2D=0,*conn2DI=0;
1152   std::set<INTERP_KERNEL::NormalizedCellType> types1D,types2D;
1153   DataArrayDouble *coordsTmp=0,*coordsTmp2=0;
1154   MCAuto<DataArrayInt> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=DataArrayInt::New(); ret1D->alloc(0,1);
1155   MCAuto<DataArrayInt> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1156   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1157   MCAuto<DataArrayInt> ret2D=m2D->convertLinearCellsToQuadratic2D1(conn2D,conn2DI,coordsTmp2,types2D); ret2D=DataArrayInt::New(); ret2D->alloc(0,1);
1158   MCAuto<DataArrayDouble> coordsTmp2Safe(coordsTmp2);
1159   MCAuto<DataArrayInt> conn2DSafe(conn2D),conn2DISafe(conn2DI);
1160   const int *c1DPtr=conn1D->begin(),*c1DIPtr=conn1DI->begin(),*c2DPtr=conn2D->begin(),*c2DIPtr=conn2DI->begin();
1161   int nbOfCells=getNumberOfCells();
1162   const int *cPtr=_nodal_connec->begin();
1163   const int *icPtr=_nodal_connec_index->begin();
1164   int lastVal=0,offset=coordsTmpSafe->getNumberOfTuples();
1165   for(int i=0;i<nbOfCells;i++,icPtr++,descIPtr++,desc2IPtr++)
1166     {
1167       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1168       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1169       if(!cm.isQuadratic())
1170         {
1171           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1172           if(typ2==INTERP_KERNEL::NORM_ERROR)
1173             {
1174               std::ostringstream oss; oss << "MEDCouplingUMesh::convertLinearCellsToQuadratic3D1 : On cell #" << i << " the linear cell type does not support advanced quadratization !";
1175               throw INTERP_KERNEL::Exception(oss.str());
1176             }
1177           types.insert(typ2); newConn->pushBackSilent(typ2);
1178           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1179           for(const int *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1180             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1181           for(const int *d=desc2Ptr+desc2IPtr[0];d!=desc2Ptr+desc2IPtr[1];d++)
1182             {
1183               int nodeId2=c2DPtr[c2DIPtr[(*d)+1]-1];
1184               int tmpPos=newConn->getNumberOfTuples();
1185               newConn->pushBackSilent(nodeId2);
1186               ret2D->pushBackSilent(nodeId2); ret1D->pushBackSilent(tmpPos);
1187             }
1188           newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1189           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+(desc2IPtr[1]-desc2IPtr[0])+1;
1190           newConnI->pushBackSilent(lastVal);
1191           ret->pushBackSilent(i);
1192         }
1193       else
1194         {
1195           types.insert(typ);
1196           lastVal+=(icPtr[1]-icPtr[0]);
1197           newConnI->pushBackSilent(lastVal);
1198           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1199         }
1200     }
1201   MCAuto<DataArrayInt> diffRet2D=ret2D->getDifferentValues();
1202   MCAuto<DataArrayInt> o2nRet2D=diffRet2D->invertArrayN2O2O2N(coordsTmp2Safe->getNumberOfTuples());
1203   coordsTmp2Safe=coordsTmp2Safe->selectByTupleId(diffRet2D->begin(),diffRet2D->end());
1204   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1205   std::vector<const DataArrayDouble *> v(3); v[0]=coordsTmpSafe; v[1]=coordsTmp2Safe; v[2]=tmp;
1206   int *c=newConn->getPointer();
1207   const int *cI(newConnI->begin());
1208   for(const int *elt=ret1D->begin();elt!=ret1D->end();elt++)
1209     c[*elt]=o2nRet2D->getIJ(c[*elt],0)+offset;
1210   offset=coordsTmp2Safe->getNumberOfTuples();
1211   for(const int *elt=ret->begin();elt!=ret->end();elt++)
1212     c[cI[(*elt)+1]-1]+=offset;
1213   coords=DataArrayDouble::Aggregate(v); conn=newConn.retn(); connI=newConnI.retn();
1214   return ret.retn();
1215 }
1216
1217 DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMeshLinear(const MEDCouplingUMesh *skin, const DataArrayInt *n2o) const
1218 {
1219   int nbOfNodesExpected(skin->getNumberOfNodes());
1220   const int *n2oPtr(n2o->begin());
1221   MCAuto<DataArrayInt> revNodal(DataArrayInt::New()),revNodalI(DataArrayInt::New());
1222   skin->getReverseNodalConnectivity(revNodal,revNodalI);
1223   const int *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1224   const int *nodalPtr(skin->getNodalConnectivity()->begin());
1225   const int *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1226   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(nbOfNodesExpected+1,1);
1227   int *work(ret->getPointer());  *work++=INTERP_KERNEL::NORM_POLYGON;
1228   if(nbOfNodesExpected<1)
1229     return ret.retn();
1230   int prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1231   *work++=n2oPtr[prevNode];
1232   for(int i=1;i<nbOfNodesExpected;i++)
1233     {
1234       if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==3)
1235         {
1236           std::set<int> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1237           conn.erase(prevNode);
1238           if(conn.size()==1)
1239             {
1240               int curNode(*(conn.begin()));
1241               *work++=n2oPtr[curNode];
1242               std::set<int> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1243               shar.erase(prevCell);
1244               if(shar.size()==1)
1245                 {
1246                   prevCell=*(shar.begin());
1247                   prevNode=curNode;
1248                 }
1249               else
1250                 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 2 !");
1251             }
1252           else
1253             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 1 !");
1254         }
1255       else
1256         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected cell !");
1257     }
1258   return ret.retn();
1259 }
1260
1261 DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayInt *n2o) const
1262 {
1263   int nbOfNodesExpected(skin->getNumberOfNodes());
1264   int nbOfTurn(nbOfNodesExpected/2);
1265   const int *n2oPtr(n2o->begin());
1266   MCAuto<DataArrayInt> revNodal(DataArrayInt::New()),revNodalI(DataArrayInt::New());
1267   skin->getReverseNodalConnectivity(revNodal,revNodalI);
1268   const int *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1269   const int *nodalPtr(skin->getNodalConnectivity()->begin());
1270   const int *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1271   MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(nbOfNodesExpected+1,1);
1272   int *work(ret->getPointer());  *work++=INTERP_KERNEL::NORM_QPOLYG;
1273   if(nbOfNodesExpected<1)
1274     return ret.retn();
1275   int prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1276   *work=n2oPtr[prevNode]; work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[0]+3]]; work++;
1277   for(int i=1;i<nbOfTurn;i++)
1278     {
1279       if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==4)
1280         {
1281           std::set<int> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1282           conn.erase(prevNode);
1283           if(conn.size()==1)
1284             {
1285               int curNode(*(conn.begin()));
1286               *work=n2oPtr[curNode];
1287               std::set<int> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1288               shar.erase(prevCell);
1289               if(shar.size()==1)
1290                 {
1291                   int curCell(*(shar.begin()));
1292                   work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[curCell]+3]];
1293                   prevCell=curCell;
1294                   prevNode=curNode;
1295                   work++;
1296                 }
1297               else
1298                 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 2 !");
1299             }
1300           else
1301             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 1 !");
1302         }
1303       else
1304         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected cell !");
1305     }
1306   return ret.retn();
1307 }
1308
1309 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector<const MEDCouplingUMesh *>& a)
1310 {
1311   if(a.empty())
1312     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : input array must be NON EMPTY !");
1313   std::vector<const MEDCouplingUMesh *>::const_iterator it=a.begin();
1314   int meshDim=(*it)->getMeshDimension();
1315   int nbOfCells=(*it)->getNumberOfCells();
1316   int meshLgth=(*it++)->getNodalConnectivityArrayLen();
1317   for(;it!=a.end();it++)
1318     {
1319       if(meshDim!=(*it)->getMeshDimension())
1320         throw INTERP_KERNEL::Exception("Mesh dimensions mismatches, MergeUMeshes impossible !");
1321       nbOfCells+=(*it)->getNumberOfCells();
1322       meshLgth+=(*it)->getNodalConnectivityArrayLen();
1323     }
1324   std::vector<const MEDCouplingPointSet *> aps(a.size());
1325   std::copy(a.begin(),a.end(),aps.begin());
1326   MCAuto<DataArrayDouble> pts=MergeNodesArray(aps);
1327   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("merge",meshDim);
1328   ret->setCoords(pts);
1329   MCAuto<DataArrayInt> c=DataArrayInt::New();
1330   c->alloc(meshLgth,1);
1331   int *cPtr=c->getPointer();
1332   MCAuto<DataArrayInt> cI=DataArrayInt::New();
1333   cI->alloc(nbOfCells+1,1);
1334   int *cIPtr=cI->getPointer();
1335   *cIPtr++=0;
1336   int offset=0;
1337   int offset2=0;
1338   for(it=a.begin();it!=a.end();it++)
1339     {
1340       int curNbOfCell=(*it)->getNumberOfCells();
1341       const int *curCI=(*it)->_nodal_connec_index->begin();
1342       const int *curC=(*it)->_nodal_connec->begin();
1343       cIPtr=std::transform(curCI+1,curCI+curNbOfCell+1,cIPtr,std::bind2nd(std::plus<int>(),offset));
1344       for(int j=0;j<curNbOfCell;j++)
1345         {
1346           const int *src=curC+curCI[j];
1347           *cPtr++=*src++;
1348           for(;src!=curC+curCI[j+1];src++,cPtr++)
1349             {
1350               if(*src!=-1)
1351                 *cPtr=*src+offset2;
1352               else
1353                 *cPtr=-1;
1354             }
1355         }
1356       offset+=curCI[curNbOfCell];
1357       offset2+=(*it)->getNumberOfNodes();
1358     }
1359   //
1360   ret->setConnectivity(c,cI,true);
1361   return ret.retn();
1362 }
1363
1364
1365 /*!
1366  * \param [in] pt the start pointer (included) of the coordinates of the point
1367  * \param [in] cellIdsBg the start pointer (included) of cellIds
1368  * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1369  * \param [in] nc nodal connectivity
1370  * \param [in] ncI nodal connectivity index
1371  * \param [in,out] ret0 the min distance between \a this and the external input point
1372  * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1373  * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1374  */
1375 void MEDCouplingUMesh::DistanceToPoint3DSurfAlg(const double *pt, const int *cellIdsBg, const int *cellIdsEnd, const double *coords, const int *nc, const int *ncI, double& ret0, int& cellId)
1376 {
1377   cellId=-1;
1378   ret0=std::numeric_limits<double>::max();
1379   for(const int *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1380     {
1381       switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1382       {
1383         case INTERP_KERNEL::NORM_TRI3:
1384           {
1385             double tmp=INTERP_KERNEL::DistanceFromPtToTriInSpaceDim3(pt,coords+3*nc[ncI[*zeCell]+1],coords+3*nc[ncI[*zeCell]+2],coords+3*nc[ncI[*zeCell]+3]);
1386             if(tmp<ret0)
1387               { ret0=tmp; cellId=*zeCell; }
1388             break;
1389           }
1390         case INTERP_KERNEL::NORM_QUAD4:
1391         case INTERP_KERNEL::NORM_POLYGON:
1392           {
1393             double tmp=INTERP_KERNEL::DistanceFromPtToPolygonInSpaceDim3(pt,nc+ncI[*zeCell]+1,nc+ncI[*zeCell+1],coords);
1394             if(tmp<ret0)
1395               { ret0=tmp; cellId=*zeCell; }
1396             break;
1397           }
1398         default:
1399           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint3DSurfAlg : not managed cell type ! Supporting TRI3, QUAD4 and POLYGON !");
1400       }
1401     }
1402 }
1403
1404 /*!
1405  * \param [in] pt the start pointer (included) of the coordinates of the point
1406  * \param [in] cellIdsBg the start pointer (included) of cellIds
1407  * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1408  * \param [in] nc nodal connectivity
1409  * \param [in] ncI nodal connectivity index
1410  * \param [in,out] ret0 the min distance between \a this and the external input point
1411  * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1412  * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1413  */
1414 void MEDCouplingUMesh::DistanceToPoint2DCurveAlg(const double *pt, const int *cellIdsBg, const int *cellIdsEnd, const double *coords, const int *nc, const int *ncI, double& ret0, int& cellId)
1415 {
1416   cellId=-1;
1417   ret0=std::numeric_limits<double>::max();
1418   for(const int *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1419     {
1420       switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1421       {
1422         case INTERP_KERNEL::NORM_SEG2:
1423           {
1424             std::size_t uselessEntry=0;
1425             double tmp=INTERP_KERNEL::SquareDistanceFromPtToSegInSpaceDim2(pt,coords+2*nc[ncI[*zeCell]+1],coords+2*nc[ncI[*zeCell]+2],uselessEntry);
1426             tmp=sqrt(tmp);
1427             if(tmp<ret0)
1428               { ret0=tmp; cellId=*zeCell; }
1429             break;
1430           }
1431         default:
1432           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint2DCurveAlg : not managed cell type ! Supporting SEG2 !");
1433       }
1434     }
1435 }
1436 DataArrayInt *MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg(std::vector<bool>& fetched, const int *seedBg, const int *seedEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn, int nbOfDepthPeeling, int& nbOfDepthPeelingPerformed)
1437 {
1438   nbOfDepthPeelingPerformed=0;
1439   if(!seedBg || !seedEnd || !arrIn || !arrIndxIn)
1440     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg : some input pointer is NULL !");
1441   int nbOfTuples=arrIndxIn->getNumberOfTuples()-1;
1442   std::vector<bool> fetched2(nbOfTuples,false);
1443   int i=0;
1444   for(const int *seedElt=seedBg;seedElt!=seedEnd;seedElt++,i++)
1445     {
1446       if(*seedElt>=0 && *seedElt<nbOfTuples)
1447         { fetched[*seedElt]=true; fetched2[*seedElt]=true; }
1448       else
1449         { std::ostringstream oss; oss << "MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg : At pos #" << i << " of seeds value is " << *seedElt << "! Should be in [0," << nbOfTuples << ") !"; throw INTERP_KERNEL::Exception(oss.str()); }
1450     }
1451   const int *arrInPtr=arrIn->begin();
1452   const int *arrIndxPtr=arrIndxIn->begin();
1453   int targetNbOfDepthPeeling=nbOfDepthPeeling!=-1?nbOfDepthPeeling:std::numeric_limits<int>::max();
1454   std::vector<int> idsToFetch1(seedBg,seedEnd);
1455   std::vector<int> idsToFetch2;
1456   std::vector<int> *idsToFetch=&idsToFetch1;
1457   std::vector<int> *idsToFetchOther=&idsToFetch2;
1458   while(!idsToFetch->empty() && nbOfDepthPeelingPerformed<targetNbOfDepthPeeling)
1459     {
1460       for(std::vector<int>::const_iterator it=idsToFetch->begin();it!=idsToFetch->end();it++)
1461         for(const int *it2=arrInPtr+arrIndxPtr[*it];it2!=arrInPtr+arrIndxPtr[*it+1];it2++)
1462           if(!fetched[*it2])
1463             { fetched[*it2]=true; fetched2[*it2]=true; idsToFetchOther->push_back(*it2); }
1464       std::swap(idsToFetch,idsToFetchOther);
1465       idsToFetchOther->clear();
1466       nbOfDepthPeelingPerformed++;
1467     }
1468   int lgth=(int)std::count(fetched2.begin(),fetched2.end(),true);
1469   i=0;
1470   MCAuto<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(lgth,1);
1471   int *retPtr=ret->getPointer();
1472   for(std::vector<bool>::const_iterator it=fetched2.begin();it!=fetched2.end();it++,i++)
1473     if(*it)
1474       *retPtr++=i;
1475   return ret.retn();
1476 }
1477
1478 /*!
1479  * This method put in zip format into parameter 'zipFrmt' in full interlace mode.
1480  * This format is often asked by INTERP_KERNEL algorithms to avoid many indirections into coordinates array.
1481  */
1482 void MEDCouplingUMesh::FillInCompact3DMode(int spaceDim, int nbOfNodesInCell, const int *conn, const double *coo, double *zipFrmt)
1483 {
1484   double *w=zipFrmt;
1485   if(spaceDim==3)
1486     for(int i=0;i<nbOfNodesInCell;i++)
1487       w=std::copy(coo+3*conn[i],coo+3*conn[i]+3,w);
1488   else if(spaceDim==2)
1489     {
1490       for(int i=0;i<nbOfNodesInCell;i++)
1491         {
1492           w=std::copy(coo+2*conn[i],coo+2*conn[i]+2,w);
1493           *w++=0.;
1494         }
1495     }
1496   else
1497     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::FillInCompact3DMode : Invalid spaceDim specified : must be 2 or 3 !");
1498 }
1499
1500 /*!
1501  * This method takes in input a cell defined by its MEDcouplingUMesh connectivity [ \a connBg , \a connEnd ) and returns its extruded cell by inserting the result at the end of ret.
1502  * \param nbOfNodesPerLev in parameter that specifies the number of nodes of one slice of global dataset
1503  * \param isQuad specifies the policy of connectivity.
1504  * @ret in/out parameter in which the result will be append
1505  */
1506 void MEDCouplingUMesh::AppendExtrudedCell(const int *connBg, const int *connEnd, int nbOfNodesPerLev, bool isQuad, std::vector<int>& ret)
1507 {
1508   INTERP_KERNEL::NormalizedCellType flatType=(INTERP_KERNEL::NormalizedCellType)connBg[0];
1509   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(flatType);
1510   ret.push_back(cm.getExtrudedType());
1511   int deltaz=isQuad?2*nbOfNodesPerLev:nbOfNodesPerLev;
1512   switch(flatType)
1513   {
1514     case INTERP_KERNEL::NORM_POINT1:
1515       {
1516         ret.push_back(connBg[1]);
1517         ret.push_back(connBg[1]+nbOfNodesPerLev);
1518         break;
1519       }
1520     case INTERP_KERNEL::NORM_SEG2:
1521       {
1522         int conn[4]={connBg[1],connBg[2],connBg[2]+deltaz,connBg[1]+deltaz};
1523         ret.insert(ret.end(),conn,conn+4);
1524         break;
1525       }
1526     case INTERP_KERNEL::NORM_SEG3:
1527       {
1528         int conn[8]={connBg[1],connBg[3],connBg[3]+deltaz,connBg[1]+deltaz,connBg[2],connBg[3]+nbOfNodesPerLev,connBg[2]+deltaz,connBg[1]+nbOfNodesPerLev};
1529         ret.insert(ret.end(),conn,conn+8);
1530         break;
1531       }
1532     case INTERP_KERNEL::NORM_QUAD4:
1533       {
1534         int conn[8]={connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz};
1535         ret.insert(ret.end(),conn,conn+8);
1536         break;
1537       }
1538     case INTERP_KERNEL::NORM_TRI3:
1539       {
1540         int conn[6]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz};
1541         ret.insert(ret.end(),conn,conn+6);
1542         break;
1543       }
1544     case INTERP_KERNEL::NORM_TRI6:
1545       {
1546         int conn[15]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4],connBg[5],connBg[6],connBg[4]+deltaz,connBg[5]+deltaz,connBg[6]+deltaz,
1547           connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
1548         ret.insert(ret.end(),conn,conn+15);
1549         break;
1550       }
1551     case INTERP_KERNEL::NORM_QUAD8:
1552       {
1553         int conn[20]={
1554           connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz,
1555           connBg[5],connBg[6],connBg[7],connBg[8],connBg[5]+deltaz,connBg[6]+deltaz,connBg[7]+deltaz,connBg[8]+deltaz,
1556           connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev,connBg[4]+nbOfNodesPerLev
1557         };
1558         ret.insert(ret.end(),conn,conn+20);
1559         break;
1560       }
1561     case INTERP_KERNEL::NORM_POLYGON:
1562       {
1563         std::back_insert_iterator< std::vector<int> > ii(ret);
1564         std::copy(connBg+1,connEnd,ii);
1565         *ii++=-1;
1566         std::reverse_iterator<const int *> rConnBg(connEnd);
1567         std::reverse_iterator<const int *> rConnEnd(connBg+1);
1568         std::transform(rConnBg,rConnEnd,ii,std::bind2nd(std::plus<int>(),deltaz));
1569         std::size_t nbOfRadFaces=std::distance(connBg+1,connEnd);
1570         for(std::size_t i=0;i<nbOfRadFaces;i++)
1571           {
1572             *ii++=-1;
1573             int conn[4]={connBg[(i+1)%nbOfRadFaces+1],connBg[i+1],connBg[i+1]+deltaz,connBg[(i+1)%nbOfRadFaces+1]+deltaz};
1574             std::copy(conn,conn+4,ii);
1575           }
1576         break;
1577       }
1578     default:
1579       throw INTERP_KERNEL::Exception("A flat type has been detected that has not its extruded representation !");
1580   }
1581 }
1582
1583
1584 /*!
1585  * This method is part of the Slice3D algorithm. It is the first step of assembly process, ones coordinates have been computed (by MEDCouplingUMesh::split3DCurveWithPlane method).
1586  * This method allows to compute given the status of 3D curve cells and the descending connectivity 3DSurf->3DCurve to deduce the intersection of each 3D surf cells
1587  * with a plane. The result will be put in 'cut3DSuf' out parameter.
1588  * \param [in] cut3DCurve  input paramter that gives for each 3DCurve cell if it owns fully to the plane or partially.
1589  * \param [out] nodesOnPlane, returns all the nodes that are on the plane.
1590  * \param [in] nodal3DSurf is the nodal connectivity of 3D surf mesh.
1591  * \param [in] nodalIndx3DSurf is the nodal connectivity index of 3D surf mesh.
1592  * \param [in] nodal3DCurve is the nodal connectivity of 3D curve mesh.
1593  * \param [in] nodal3DIndxCurve is the nodal connectivity index of 3D curve mesh.
1594  * \param [in] desc is the descending connectivity 3DSurf->3DCurve
1595  * \param [in] descIndx is the descending connectivity index 3DSurf->3DCurve
1596  * \param [out] cut3DSuf input/output param.
1597  */
1598 void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<int>& cut3DCurve, std::vector<int>& nodesOnPlane, const int *nodal3DSurf, const int *nodalIndx3DSurf,
1599                                                    const int *nodal3DCurve, const int *nodalIndx3DCurve,
1600                                                    const int *desc, const int *descIndx,
1601                                                    std::vector< std::pair<int,int> >& cut3DSurf)
1602 {
1603   std::set<int> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
1604   int nbOf3DSurfCell=(int)cut3DSurf.size();
1605   for(int i=0;i<nbOf3DSurfCell;i++)
1606     {
1607       std::vector<int> res;
1608       int offset=descIndx[i];
1609       int nbOfSeg=descIndx[i+1]-offset;
1610       for(int j=0;j<nbOfSeg;j++)
1611         {
1612           int edgeId=desc[offset+j];
1613           int status=cut3DCurve[edgeId];
1614           if(status!=-2)
1615             {
1616               if(status>-1)
1617                 res.push_back(status);
1618               else
1619                 {
1620                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
1621                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
1622                 }
1623             }
1624         }
1625       switch(res.size())
1626       {
1627         case 2:
1628           {
1629             cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1630             break;
1631           }
1632         case 1:
1633         case 0:
1634           {
1635             std::set<int> s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]);
1636             std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector<int> >(res));
1637             if(res.size()==2)
1638               {
1639                 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1640               }
1641             else
1642               {
1643                 cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
1644               }
1645             break;
1646           }
1647         default:
1648           {// case when plane is on a multi colinear edge of a polyhedron
1649             if((int)res.size()==2*nbOfSeg)
1650               {
1651                 cut3DSurf[i].first=-2; cut3DSurf[i].second=i;
1652               }
1653             else
1654               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
1655           }
1656       }
1657     }
1658 }
1659
1660
1661 /*!
1662  * \a this is expected to be a mesh with spaceDim==3 and meshDim==3. If not an exception will be thrown.
1663  * This method is part of the Slice3D algorithm. It is the second step of assembly process, ones coordinates have been computed (by MEDCouplingUMesh::split3DCurveWithPlane method).
1664  * This method allows to compute given the result of 3D surf cells with plane and the descending connectivity 3D->3DSurf to deduce the intersection of each 3D cells
1665  * with a plane. The result will be put in 'nodalRes' 'nodalResIndx' and 'cellIds' out parameters.
1666  * \param cut3DSurf  input paramter that gives for each 3DSurf its intersection with plane (result of MEDCouplingUMesh::AssemblyForSplitFrom3DCurve).
1667  * \param desc is the descending connectivity 3D->3DSurf
1668  * \param descIndx is the descending connectivity index 3D->3DSurf
1669  */
1670 void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<int,int> >& cut3DSurf,
1671                                                   const int *desc, const int *descIndx,
1672                                                   DataArrayInt *nodalRes, DataArrayInt *nodalResIndx, DataArrayInt *cellIds) const
1673 {
1674   checkFullyDefined();
1675   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1676     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1677   const int *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1678   int nbOfCells(getNumberOfCells());
1679   for(int i=0;i<nbOfCells;i++)
1680     {
1681       std::map<int, std::set<int> > m;
1682       int offset=descIndx[i];
1683       int nbOfFaces=descIndx[i+1]-offset;
1684       int start=-1;
1685       int end=-1;
1686       for(int j=0;j<nbOfFaces;j++)
1687         {
1688           const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
1689           if(p.first!=-1 && p.second!=-1)
1690             {
1691               if(p.first!=-2)
1692                 {
1693                   start=p.first; end=p.second;
1694                   m[p.first].insert(p.second);
1695                   m[p.second].insert(p.first);
1696                 }
1697               else
1698                 {
1699                   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]);
1700                   int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1701                   INTERP_KERNEL::AutoPtr<int> tmp=new int[sz];
1702                   INTERP_KERNEL::NormalizedCellType cmsId;
1703                   unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId);
1704                   start=tmp[0]; end=tmp[nbOfNodesSon-1];
1705                   for(unsigned k=0;k<nbOfNodesSon;k++)
1706                     {
1707                       m[tmp[k]].insert(tmp[(k+1)%nbOfNodesSon]);
1708                       m[tmp[(k+1)%nbOfNodesSon]].insert(tmp[k]);
1709                     }
1710                 }
1711             }
1712         }
1713       if(m.empty())
1714         continue;
1715       std::vector<int> conn(1,(int)INTERP_KERNEL::NORM_POLYGON);
1716       int prev=end;
1717       while(end!=start)
1718         {
1719           std::map<int, std::set<int> >::const_iterator it=m.find(start);
1720           const std::set<int>& s=(*it).second;
1721           std::set<int> s2; s2.insert(prev);
1722           std::set<int> s3;
1723           std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
1724           if(s3.size()==1)
1725             {
1726               int val=*s3.begin();
1727               conn.push_back(start);
1728               prev=start;
1729               start=val;
1730             }
1731           else
1732             start=end;
1733         }
1734       conn.push_back(end);
1735       if(conn.size()>3)
1736         {
1737           nodalRes->insertAtTheEnd(conn.begin(),conn.end());
1738           nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples());
1739           cellIds->pushBackSilent(i);
1740         }
1741     }
1742 }
1743
1744
1745 void MEDCouplingUMesh::ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayInt *nodalConnec, const DataArrayInt *nodalConnecIndex)
1746 {
1747   if(nodalConnec && nodalConnecIndex)
1748     {
1749       types.clear();
1750       const int *conn(nodalConnec->begin()),*connIndex(nodalConnecIndex->begin());
1751       int nbOfElem(nodalConnecIndex->getNbOfElems()-1);
1752       if(nbOfElem>0)
1753         for(const int *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
1754           types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
1755     }
1756 }