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::QuadraticPlanarArcDetectionPrecision arcPrec(1.e-10);  // RAII
 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 }