src/MEDCoupling/MEDCouplingUMesh_internal.cxx

   1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingUMesh.hxx"
  22 #include "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       mcIdType nbOfCells=getNumberOfCells();
  78       const mcIdType *c=_nodal_connec->getConstPointer();
  79       const mcIdType *ci=_nodal_connec_index->getConstPointer();
  80       for(mcIdType 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 DataArrayIdType *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   mcIdType nbOfCells=getNumberOfCells();
 102   MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
 103   mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
 104   ret->alloc(nbOfCells+nbOfCutCells,1);
 105   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
 106   mcIdType *retPt=ret->getPointer();
 107   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
 108   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
 109   newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
 110   newConn->alloc(getNodalConnectivityArrayLen()+3*nbOfCutCells,1);
 111   mcIdType *pt=newConn->getPointer();
 112   mcIdType *ptI=newConnI->getPointer();
 113   ptI[0]=0;
 114   const mcIdType *oldc=_nodal_connec->begin();
 115   const mcIdType *ci=_nodal_connec_index->begin();
 116   for(mcIdType i=0;i<nbOfCells;i++,ci++)
 117     {
 118       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
 119         {
 120           const mcIdType tmp[8]={ToIdType(INTERP_KERNEL::NORM_TRI3),oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+3],
 121                                  ToIdType(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 DataArrayIdType *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   mcIdType nbOfCells=getNumberOfCells();
 155   MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
 156   mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
 157   ret->alloc(nbOfCells+nbOfCutCells,1);
 158   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
 159   mcIdType *retPt=ret->getPointer();
 160   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
 161   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
 162   newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
 163   newConn->alloc(getNodalConnectivityArrayLen()+3*nbOfCutCells,1);
 164   mcIdType *pt=newConn->getPointer();
 165   mcIdType *ptI=newConnI->getPointer();
 166   ptI[0]=0;
 167   const mcIdType *oldc=_nodal_connec->begin();
 168   const mcIdType *ci=_nodal_connec_index->begin();
 169   for(mcIdType i=0;i<nbOfCells;i++,ci++)
 170     {
 171       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
 172         {
 173           const mcIdType tmp[8]={ToIdType(INTERP_KERNEL::NORM_TRI3),oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+4],
 174                                  ToIdType(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 DataArrayIdType *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   mcIdType nbOfCells=getNumberOfCells();
 208   MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
 209   mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_HEXA8);
 210   ret->alloc(nbOfCells+4*nbOfCutCells,1);
 211   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
 212   mcIdType *retPt=ret->getPointer();
 213   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
 214   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
 215   newConnI->alloc(nbOfCells+4*nbOfCutCells+1,1);
 216   newConn->alloc(getNodalConnectivityArrayLen()+16*nbOfCutCells,1);//21
 217   mcIdType *pt=newConn->getPointer();
 218   mcIdType *ptI=newConnI->getPointer();
 219   ptI[0]=0;
 220   const mcIdType *oldc=_nodal_connec->begin();
 221   const mcIdType *ci=_nodal_connec_index->begin();
 222   for(mcIdType 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]=ToIdType(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 DataArrayIdType *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   mcIdType nbOfCells=getNumberOfCells();
 260   MCAuto<DataArrayIdType> ret=DataArrayIdType::New();
 261   mcIdType nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_HEXA8);
 262   ret->alloc(nbOfCells+5*nbOfCutCells,1);
 263   if(nbOfCutCells==0) { ret->iota(0); return ret.retn(); }
 264   mcIdType *retPt=ret->getPointer();
 265   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New();
 266   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New();
 267   newConnI->alloc(nbOfCells+5*nbOfCutCells+1,1);
 268   newConn->alloc(getNodalConnectivityArrayLen()+21*nbOfCutCells,1);
 269   mcIdType *pt=newConn->getPointer();
 270   mcIdType *ptI=newConnI->getPointer();
 271   ptI[0]=0;
 272   const mcIdType *oldc=_nodal_connec->begin();
 273   const mcIdType *ci=_nodal_connec_index->begin();
 274   for(mcIdType 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]=ToIdType(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<DataArrayIdType> desc1(DataArrayIdType::New()),descIndx1(DataArrayIdType::New()),revDesc1(DataArrayIdType::New()),revDescIndx1(DataArrayIdType::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::QuadraticPlanarPrecision arcPrec(1.e-10);  // RAII
 352   mcIdType nbCells=getNumberOfCells();
 353   mcIdType nbNodes=getNumberOfNodes();
 354   const mcIdType *conn=_nodal_connec->begin();
 355   const mcIdType *connI=_nodal_connec_index->begin();
 356   const double *coords=_coords->begin();
 357   std::vector<double> addCoo;
 358   std::vector<mcIdType> newConn;//no direct DataArrayIdType because interface with Geometric2D
 359   MCAuto<DataArrayIdType> newConnI(DataArrayIdType::New());
 360   newConnI->alloc(nbCells+1,1);
 361   mcIdType *newConnIPtr=newConnI->getPointer();
 362   *newConnIPtr=0;
 363   mcIdType tmp1[3];
 364   INTERP_KERNEL::Node *tmp2[3];
 365   std::set<INTERP_KERNEL::NormalizedCellType> types;
 366   for(mcIdType 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]=ToIdType(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() && ToIdType(newConn.size())==_nodal_connec->getNumberOfTuples())//nothing happens during tessellation : no update needed
 399     return ;
 400   _types=types;
 401   DataArrayIdType::SetArrayIn(newConnI,_nodal_connec_index);
 402   MCAuto<DataArrayIdType> newConnArr=DataArrayIdType::New();
 403   newConnArr->alloc(newConn.size(),1);
 404   std::copy(newConn.begin(),newConn.end(),newConnArr->getPointer());
 405   DataArrayIdType::SetArrayIn(newConnArr,_nodal_connec);
 406   MCAuto<DataArrayDouble> newCoords=DataArrayDouble::New();
 407   newCoords->alloc(nbNodes+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 completely ignores 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 mcIdType *nodeSubdived, const mcIdType *nodeIndxSubdived, const mcIdType *desc, const mcIdType *descIndex)
 424 {
 425   checkFullyDefined();
 426   if(getMeshDimension()!=2)
 427     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : works only on umesh with meshdim==2 !");
 428   mcIdType nbOfCells=getNumberOfCells();
 429   mcIdType *connI=_nodal_connec_index->getPointer();
 430   mcIdType newConnLgth=0;
 431   for(mcIdType i=0;i<nbOfCells;i++,connI++)
 432     {
 433       mcIdType offset=descIndex[i];
 434       mcIdType nbOfEdges=descIndex[i+1]-offset;
 435       //
 436       bool ddirect=desc[offset+nbOfEdges-1]>0;
 437       mcIdType eedgeId=std::abs(desc[offset+nbOfEdges-1])-1;
 438       mcIdType ref=ddirect?nodeSubdived[nodeIndxSubdived[eedgeId+1]-1]:nodeSubdived[nodeIndxSubdived[eedgeId]+1];
 439       for(mcIdType j=0;j<nbOfEdges;j++)
 440         {
 441           bool direct=desc[offset+j]>0;
 442           mcIdType edgeId=std::abs(desc[offset+j])-1;
 443           if(!INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodeSubdived[nodeIndxSubdived[edgeId]]).isQuadratic())
 444             {
 445               mcIdType id1=nodeSubdived[nodeIndxSubdived[edgeId]+1];
 446               mcIdType id2=nodeSubdived[nodeIndxSubdived[edgeId+1]-1];
 447               mcIdType ref2=direct?id1:id2;
 448               if(ref==ref2)
 449                 {
 450                   mcIdType 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<DataArrayIdType> newConn=DataArrayIdType::New();
 470   newConn->alloc(newConnLgth,1);
 471   mcIdType *work=newConn->getPointer();
 472   for(mcIdType i=0;i<nbOfCells;i++)
 473     {
 474       *work++=INTERP_KERNEL::NORM_POLYGON;
 475       mcIdType offset=descIndex[i];
 476       mcIdType nbOfEdges=descIndex[i+1]-offset;
 477       for(mcIdType j=0;j<nbOfEdges;j++)
 478         {
 479           bool direct=desc[offset+j]>0;
 480           mcIdType 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               mcIdType nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1;
 486               std::reverse_iterator<const mcIdType *> it(nodeSubdived+nodeIndxSubdived[edgeId+1]);
 487               work=std::copy(it,it+nbOfSubNodes-1,work);
 488             }
 489         }
 490     }
 491   DataArrayIdType::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 mcIdType *begin, const mcIdType *end, bool fullyIn, DataArrayIdType *&cellIdsKeptArr) const
 509 {
 510   MCAuto<DataArrayIdType> cellIdsKept=DataArrayIdType::New(); cellIdsKept->alloc(0,1);
 511   checkConnectivityFullyDefined();
 512   mcIdType tmp=-1;
 513   mcIdType sz=getNodalConnectivity()->getMaxValue(tmp); sz=std::max(sz,ToIdType(0))+1;
 514   std::vector<bool> fastFinder(sz,false);
 515   for(const mcIdType *work=begin;work!=end;work++)
 516     if(*work>=0 && *work<sz)
 517       fastFinder[*work]=true;
 518   mcIdType nbOfCells=getNumberOfCells();
 519   const mcIdType *conn=getNodalConnectivity()->getConstPointer();
 520   const mcIdType *connIndex=getNodalConnectivityIndex()->getConstPointer();
 521   for(mcIdType i=0;i<nbOfCells;i++)
 522     {
 523       mcIdType ref=0,nbOfHit=0;
 524       for(const mcIdType *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<mcIdType>& 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   mcIdType ncells=getNumberOfCells();
 554   mcIdType 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 mcIdType *conn=_nodal_connec->getConstPointer();
 561   const mcIdType *connI=_nodal_connec_index->getConstPointer();
 562   const double *coo=_coords->getConstPointer();
 563   std::vector<double> addCoo;
 564   for(mcIdType i=0;i<ncells;i++)
 565     {
 566       if(conn[connI[i]]==ToIdType(INTERP_KERNEL::NORM_SEG2))
 567         {
 568           if(cut3DCurve[i]==-2)
 569             {
 570               mcIdType 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                       mcIdType nNode=ToIdType(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       mcIdType newNbOfNodes=nnodes+ToIdType(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   mcIdType oldNbOfNodes=getNumberOfNodes();
 612   mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
 613   std::size_t spaceDim=getSpaceDimension();
 614   DataArrayDouble *ret=DataArrayDouble::New();
 615   std::vector<bool> isQuads;
 616   mcIdType 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<mcIdType> v;
 622   std::vector<double> c;
 623   double vec[3];
 624   v.reserve(3);
 625   c.reserve(6);
 626   for(mcIdType 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(mcIdType 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(mcIdType 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   mcIdType oldNbOfNodes=getNumberOfNodes();
 674   mcIdType nbOf1DCells=ToIdType(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   mcIdType 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 mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
 687   const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
 688   for(mcIdType 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   mcIdType oldNbOfNodes=getNumberOfNodes();
 718   mcIdType nbOf1DCells=ToIdType(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   mcIdType 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 mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
 731   const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
 732   for(mcIdType 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(mcIdType nbOfNodesOf1Lev, bool isQuad) const
 781 {
 782   mcIdType nbOf1DCells(getNumberOfNodes()/nbOfNodesOf1Lev-1);
 783   mcIdType nbOf2DCells=getNumberOfCells();
 784   mcIdType nbOf3DCells(nbOf2DCells*nbOf1DCells);
 785   MEDCouplingUMesh *ret(MEDCouplingUMesh::New("Extruded",getMeshDimension()+1));
 786   const mcIdType *conn(_nodal_connec->begin()),*connI(_nodal_connec_index->begin());
 787   MCAuto<DataArrayIdType> newConn(DataArrayIdType::New()),newConnI(DataArrayIdType::New());
 788   newConnI->alloc(nbOf3DCells+1,1);
 789   mcIdType *newConnIPtr(newConnI->getPointer());
 790   *newConnIPtr++=0;
 791   std::vector<mcIdType> newc;
 792   for(mcIdType j=0;j<nbOf2DCells;j++)
 793     {
 794       AppendExtrudedCell(conn+connI[j],conn+connI[j+1],nbOfNodesOf1Lev,isQuad,newc);
 795       *newConnIPtr++=ToIdType(newc.size());
 796     }
 797   newConn->alloc(newc.size()*nbOf1DCells,1);
 798   mcIdType *newConnPtr(newConn->getPointer());
 799   mcIdType deltaPerLev(isQuad?2*nbOfNodesOf1Lev:nbOfNodesOf1Lev);
 800   newConnIPtr=newConnI->getPointer();
 801   for(mcIdType 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<mcIdType>(),newConnIPtr[iz*nbOf2DCells]));
 805       const mcIdType *posOfTypeOfCell(newConnIPtr);
 806       for(std::vector<mcIdType>::const_iterator iter=newc.begin();iter!=newc.end();iter++,newConnPtr++)
 807         {
 808           mcIdType icell(ToIdType(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<mcIdType>& candidates, int compType, const mcIdType *conn, const mcIdType *connI, DataArrayIdType *result)
 838 {
 839   if(candidates.size()<1)
 840     return false;
 841   bool ret=false;
 842   std::vector<mcIdType>::const_iterator iter=candidates.begin();
 843   mcIdType 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 mcIdType *begin, const mcIdType *end) const
 869 {
 870   checkConnectivityFullyDefined();
 871   mcIdType 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   mcIdType *connIndexRet=(mcIdType *)malloc((nbOfElemsRet+1)*sizeof(mcIdType));
 877   connIndexRet[0]=0;
 878   const mcIdType *conn=_nodal_connec->getConstPointer();
 879   const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
 880   mcIdType newNbring=0;
 881   for(const mcIdType *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   mcIdType *connRet=(mcIdType *)malloc(connIndexRet[nbOfElemsRet]*sizeof(mcIdType));
 893   mcIdType *connRetWork=connRet;
 894   std::set<INTERP_KERNEL::NormalizedCellType> types;
 895   for(const mcIdType *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<DataArrayIdType> connRetArr=DataArrayIdType::New();
 901   connRetArr->useArray(connRet,true,DeallocType::C_DEALLOC,connIndexRet[nbOfElemsRet],1);
 902   MCAuto<DataArrayIdType> connIndexRetArr=DataArrayIdType::New();
 903   connIndexRetArr->useArray(connIndexRet,true,DeallocType::C_DEALLOC,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(mcIdType start, mcIdType end, mcIdType step) const
 915 {
 916   checkFullyDefined();
 917   mcIdType ncell=getNumberOfCells();
 918   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
 919   ret->_mesh_dim=_mesh_dim;
 920   ret->setCoords(_coords);
 921   mcIdType newNbOfCells=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : ");
 922   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(newNbOfCells+1,1);
 923   mcIdType *newConnIPtr=newConnI->getPointer(); *newConnIPtr=0;
 924   mcIdType work=start;
 925   const mcIdType *conn=_nodal_connec->getConstPointer();
 926   const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
 927   for(mcIdType 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<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(newConnIPtr[0],1);
 940   mcIdType *newConnPtr=newConn->getPointer();
 941   std::set<INTERP_KERNEL::NormalizedCellType> types;
 942   work=start;
 943   for(mcIdType 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 mcIdType MEDCouplingFastNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2)
 956 {
 957   return id;
 958 }
 959
 960 mcIdType MEDCouplingOrientationSensitiveNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *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  * Implements \a conversionType 0 for meshes with meshDim = 1, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
 976  * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
 977  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
 978  */
 979 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic1D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
 980 {
 981   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
 982   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
 983   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
 984   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
 985   mcIdType nbOfCells=getNumberOfCells();
 986   mcIdType nbOfNodes=getNumberOfNodes();
 987   const mcIdType *cPtr=_nodal_connec->begin();
 988   const mcIdType *icPtr=_nodal_connec_index->begin();
 989   mcIdType lastVal=0,offset=nbOfNodes;
 990   for(mcIdType 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(ToIdType(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 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2DAnd3D0(const MEDCouplingUMesh *m1D, const DataArrayIdType *desc, const DataArrayIdType *descI, DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1017 {
1018   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1019   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1020   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1021   //
1022   const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1023   DataArrayIdType *conn1D=0,*conn1DI=0;
1024   std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1025   DataArrayDouble *coordsTmp=0;
1026   MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1027   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1028   MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1029   const mcIdType *c1DPtr=conn1D->begin();
1030   const mcIdType *c1DIPtr=conn1DI->begin();
1031   mcIdType nbOfCells=getNumberOfCells();
1032   const mcIdType *cPtr=_nodal_connec->begin();
1033   const mcIdType *icPtr=_nodal_connec_index->begin();
1034   mcIdType lastVal=0;
1035   for(mcIdType 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 mcIdType *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  * Implements \a conversionType 0 for meshes with meshDim = 2, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1064  * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1065  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1066  */
1067 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1068 {
1069   MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::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 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1075 {
1076   MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1077   MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1078   //
1079   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1080   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1081   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1082   //
1083   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1084   const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1085   DataArrayIdType *conn1D=0,*conn1DI=0;
1086   std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1087   DataArrayDouble *coordsTmp=0;
1088   MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1089   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1090   MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1091   const mcIdType *c1DPtr=conn1D->begin();
1092   const mcIdType *c1DIPtr=conn1DI->begin();
1093   mcIdType nbOfCells=getNumberOfCells();
1094   const mcIdType *cPtr=_nodal_connec->begin();
1095   const mcIdType *icPtr=_nodal_connec_index->begin();
1096   mcIdType lastVal=0;
1097   mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1098   for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1099     {
1100       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1101       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1102       if(!cm.isQuadratic())
1103         {
1104           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1105           types.insert(typ2); newConn->pushBackSilent(typ2);
1106           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1107           for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1108             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1109           newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1110           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+1;
1111           newConnI->pushBackSilent(lastVal);
1112           ret->pushBackSilent(i);
1113         }
1114       else
1115         {
1116           types.insert(typ);
1117           lastVal+=(icPtr[1]-icPtr[0]);
1118           newConnI->pushBackSilent(lastVal);
1119           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1120         }
1121     }
1122   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1123   coords=DataArrayDouble::Aggregate(coordsTmpSafe,tmp); conn=newConn.retn(); connI=newConnI.retn();
1124   return ret.retn();
1125 }
1126
1127 /*!
1128  * Implements \a conversionType 0 for meshes with meshDim = 3, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1129  * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1130  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1131  */
1132 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1133 {
1134   MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1135   MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1136   return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1137 }
1138
1139 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1140 {
1141   MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),desc2I(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1142   MCAuto<MEDCouplingUMesh> m2D=buildDescendingConnectivityGen<MinusOneSonsGeneratorBiQuadratic>(desc2,desc2I,tmp2,tmp3,MEDCouplingFastNbrer); tmp2=0; tmp3=0;
1143   MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),desc1I(DataArrayIdType::New()),tmp4(DataArrayIdType::New()),tmp5(DataArrayIdType::New());
1144   MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc1,desc1I,tmp4,tmp5); tmp4=0; tmp5=0;
1145   //
1146   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1147   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1148   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(),ret2=DataArrayIdType::New(); ret->alloc(0,1); ret2->alloc(0,1);
1149   //
1150   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1151   const mcIdType *descPtr(desc1->begin()),*descIPtr(desc1I->begin()),*desc2Ptr(desc2->begin()),*desc2IPtr(desc2I->begin());
1152   DataArrayIdType *conn1D=0,*conn1DI=0,*conn2D=0,*conn2DI=0;
1153   std::set<INTERP_KERNEL::NormalizedCellType> types1D,types2D;
1154   DataArrayDouble *coordsTmp=0,*coordsTmp2=0;
1155   MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=DataArrayIdType::New(); ret1D->alloc(0,1);
1156   MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1157   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1158   MCAuto<DataArrayIdType> ret2D=m2D->convertLinearCellsToQuadratic2D1(conn2D,conn2DI,coordsTmp2,types2D); ret2D=DataArrayIdType::New(); ret2D->alloc(0,1);
1159   MCAuto<DataArrayDouble> coordsTmp2Safe(coordsTmp2);
1160   MCAuto<DataArrayIdType> conn2DSafe(conn2D),conn2DISafe(conn2DI);
1161   const mcIdType *c1DPtr=conn1D->begin(),*c1DIPtr=conn1DI->begin(),*c2DPtr=conn2D->begin(),*c2DIPtr=conn2DI->begin();
1162   mcIdType nbOfCells=getNumberOfCells();
1163   const mcIdType *cPtr=_nodal_connec->begin();
1164   const mcIdType *icPtr=_nodal_connec_index->begin();
1165   mcIdType lastVal=0;
1166   mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1167   for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++,desc2IPtr++)
1168     {
1169       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1170       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1171       if(!cm.isQuadratic())
1172         {
1173           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1174           if(typ2==INTERP_KERNEL::NORM_ERROR)
1175             {
1176               std::ostringstream oss; oss << "MEDCouplingUMesh::convertLinearCellsToQuadratic3D1 : On cell #" << i << " the linear cell type does not support advanced quadratization !";
1177               throw INTERP_KERNEL::Exception(oss.str());
1178             }
1179           types.insert(typ2); newConn->pushBackSilent(typ2);
1180           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1181           for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1182             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1183           for(const mcIdType *d=desc2Ptr+desc2IPtr[0];d!=desc2Ptr+desc2IPtr[1];d++)
1184             {
1185               mcIdType nodeId2=c2DPtr[c2DIPtr[(*d)+1]-1];
1186               mcIdType tmpPos=newConn->getNumberOfTuples();
1187               newConn->pushBackSilent(nodeId2);
1188               ret2D->pushBackSilent(nodeId2); ret1D->pushBackSilent(tmpPos);
1189             }
1190           newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1191           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+(desc2IPtr[1]-desc2IPtr[0])+1;
1192           newConnI->pushBackSilent(lastVal);
1193           ret->pushBackSilent(i);
1194         }
1195       else
1196         {
1197           types.insert(typ);
1198           lastVal+=(icPtr[1]-icPtr[0]);
1199           newConnI->pushBackSilent(lastVal);
1200           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1201         }
1202     }
1203   MCAuto<DataArrayIdType> diffRet2D=ret2D->getDifferentValues();
1204   MCAuto<DataArrayIdType> o2nRet2D=diffRet2D->invertArrayN2O2O2N(coordsTmp2Safe->getNumberOfTuples());
1205   coordsTmp2Safe=coordsTmp2Safe->selectByTupleId(diffRet2D->begin(),diffRet2D->end());
1206   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1207   std::vector<const DataArrayDouble *> v(3); v[0]=coordsTmpSafe; v[1]=coordsTmp2Safe; v[2]=tmp;
1208   mcIdType *c=newConn->getPointer();
1209   const mcIdType *cI(newConnI->begin());
1210   for(const mcIdType *elt=ret1D->begin();elt!=ret1D->end();elt++)
1211     c[*elt]=o2nRet2D->getIJ(c[*elt],0)+offset;
1212   offset=coordsTmp2Safe->getNumberOfTuples();
1213   for(const mcIdType *elt=ret->begin();elt!=ret->end();elt++)
1214     c[cI[(*elt)+1]-1]+=offset;
1215   coords=DataArrayDouble::Aggregate(v); conn=newConn.retn(); connI=newConnI.retn();
1216   return ret.retn();
1217 }
1218
1219 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshLinear(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1220 {
1221   mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1222   const mcIdType *n2oPtr(n2o->begin());
1223   MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1224   skin->getReverseNodalConnectivity(revNodal,revNodalI);
1225   const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1226   const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1227   const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1228   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1229   mcIdType *work(ret->getPointer());  *work++=INTERP_KERNEL::NORM_POLYGON;
1230   if(nbOfNodesExpected<1)
1231     return ret.retn();
1232   mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1233   *work++=n2oPtr[prevNode];
1234   for(mcIdType i=1;i<nbOfNodesExpected;i++)
1235     {
1236       if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==3)
1237         {
1238           std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1239           conn.erase(prevNode);
1240           if(conn.size()==1)
1241             {
1242               mcIdType curNode(*(conn.begin()));
1243               *work++=n2oPtr[curNode];
1244               std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1245               shar.erase(prevCell);
1246               if(shar.size()==1)
1247                 {
1248                   prevCell=*(shar.begin());
1249                   prevNode=curNode;
1250                 }
1251               else
1252                 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 2 !");
1253             }
1254           else
1255             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 1 !");
1256         }
1257       else
1258         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected cell !");
1259     }
1260   return ret.retn();
1261 }
1262
1263 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1264 {
1265   mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1266   mcIdType nbOfTurn(nbOfNodesExpected/2);
1267   const mcIdType *n2oPtr(n2o->begin());
1268   MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1269   skin->getReverseNodalConnectivity(revNodal,revNodalI);
1270   const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1271   const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1272   const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1273   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1274   mcIdType *work(ret->getPointer());  *work++=INTERP_KERNEL::NORM_QPOLYG;
1275   if(nbOfNodesExpected<1)
1276     return ret.retn();
1277   mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1278   *work=n2oPtr[prevNode]; work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[0]+3]]; work++;
1279   for(mcIdType i=1;i<nbOfTurn;i++)
1280     {
1281       if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==4)
1282         {
1283           std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1284           conn.erase(prevNode);
1285           if(conn.size()==1)
1286             {
1287               mcIdType curNode(*(conn.begin()));
1288               *work=n2oPtr[curNode];
1289               std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1290               shar.erase(prevCell);
1291               if(shar.size()==1)
1292                 {
1293                   mcIdType curCell(*(shar.begin()));
1294                   work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[curCell]+3]];
1295                   prevCell=curCell;
1296                   prevNode=curNode;
1297                   work++;
1298                 }
1299               else
1300                 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 2 !");
1301             }
1302           else
1303             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 1 !");
1304         }
1305       else
1306         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected cell !");
1307     }
1308   return ret.retn();
1309 }
1310
1311 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector<const MEDCouplingUMesh *>& a)
1312 {
1313   if(a.empty())
1314     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : input array must be NON EMPTY !");
1315   std::vector<const MEDCouplingUMesh *>::const_iterator it=a.begin();
1316   int meshDim=(*it)->getMeshDimension();
1317   mcIdType nbOfCells=ToIdType((*it)->getNumberOfCells());
1318   mcIdType meshLgth=(*it++)->getNodalConnectivityArrayLen();
1319   for(;it!=a.end();it++)
1320     {
1321       if(meshDim!=(*it)->getMeshDimension())
1322         throw INTERP_KERNEL::Exception("Mesh dimensions mismatches, MergeUMeshes impossible !");
1323       nbOfCells+=ToIdType((*it)->getNumberOfCells());
1324       meshLgth+=(*it)->getNodalConnectivityArrayLen();
1325     }
1326   std::vector<const MEDCouplingPointSet *> aps(a.size());
1327   std::copy(a.begin(),a.end(),aps.begin());
1328   MCAuto<DataArrayDouble> pts=MergeNodesArray(aps);
1329   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("merge",meshDim);
1330   ret->setCoords(pts);
1331   MCAuto<DataArrayIdType> c=DataArrayIdType::New();
1332   c->alloc(meshLgth,1);
1333   mcIdType *cPtr=c->getPointer();
1334   MCAuto<DataArrayIdType> cI=DataArrayIdType::New();
1335   cI->alloc(nbOfCells+1,1);
1336   mcIdType *cIPtr=cI->getPointer();
1337   *cIPtr++=0;
1338   mcIdType offset=0;
1339   mcIdType offset2=0;
1340   for(it=a.begin();it!=a.end();it++)
1341     {
1342       mcIdType curNbOfCell=ToIdType((*it)->getNumberOfCells());
1343       const mcIdType *curCI=(*it)->_nodal_connec_index->begin();
1344       const mcIdType *curC=(*it)->_nodal_connec->begin();
1345       cIPtr=std::transform(curCI+1,curCI+curNbOfCell+1,cIPtr,std::bind2nd(std::plus<mcIdType>(),offset));
1346       for(mcIdType j=0;j<curNbOfCell;j++)
1347         {
1348           const mcIdType *src=curC+curCI[j];
1349           *cPtr++=*src++;
1350           for(;src!=curC+curCI[j+1];src++,cPtr++)
1351             {
1352               if(*src!=-1)
1353                 *cPtr=*src+offset2;
1354               else
1355                 *cPtr=-1;
1356             }
1357         }
1358       offset+=curCI[curNbOfCell];
1359       offset2+=(*it)->getNumberOfNodes();
1360     }
1361   //
1362   ret->setConnectivity(c,cI,true);
1363   return ret.retn();
1364 }
1365
1366
1367 /*!
1368  * \param [in] pt the start pointer (included) of the coordinates of the point
1369  * \param [in] cellIdsBg the start pointer (included) of cellIds
1370  * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1371  * \param [in] nc nodal connectivity
1372  * \param [in] ncI nodal connectivity index
1373  * \param [in,out] ret0 the min distance between \a this and the external input point
1374  * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1375  * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1376  */
1377 void MEDCouplingUMesh::DistanceToPoint3DSurfAlg(const double *pt, const mcIdType *cellIdsBg, const mcIdType *cellIdsEnd, const double *coords, const mcIdType *nc, const mcIdType *ncI, double& ret0, mcIdType& cellId)
1378 {
1379   cellId=-1;
1380   ret0=std::numeric_limits<double>::max();
1381   for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1382     {
1383       switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1384       {
1385         case INTERP_KERNEL::NORM_TRI3:
1386           {
1387             double tmp=INTERP_KERNEL::DistanceFromPtToTriInSpaceDim3(pt,coords+3*nc[ncI[*zeCell]+1],coords+3*nc[ncI[*zeCell]+2],coords+3*nc[ncI[*zeCell]+3]);
1388             if(tmp<ret0)
1389               { ret0=tmp; cellId=*zeCell; }
1390             break;
1391           }
1392         case INTERP_KERNEL::NORM_QUAD4:
1393         case INTERP_KERNEL::NORM_POLYGON:
1394           {
1395             double tmp=INTERP_KERNEL::DistanceFromPtToPolygonInSpaceDim3(pt,nc+ncI[*zeCell]+1,nc+ncI[*zeCell+1],coords);
1396             if(tmp<ret0)
1397               { ret0=tmp; cellId=*zeCell; }
1398             break;
1399           }
1400         default:
1401           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint3DSurfAlg : not managed cell type ! Supporting TRI3, QUAD4 and POLYGON !");
1402       }
1403     }
1404 }
1405
1406 /*!
1407  * \param [in] pt the start pointer (included) of the coordinates of the point
1408  * \param [in] cellIdsBg the start pointer (included) of cellIds
1409  * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1410  * \param [in] nc nodal connectivity
1411  * \param [in] ncI nodal connectivity index
1412  * \param [in,out] ret0 the min distance between \a this and the external input point
1413  * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1414  * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1415  */
1416 void MEDCouplingUMesh::DistanceToPoint2DCurveAlg(const double *pt, const mcIdType *cellIdsBg, const mcIdType *cellIdsEnd, const double *coords, const mcIdType *nc, const mcIdType *ncI, double& ret0, mcIdType& cellId)
1417 {
1418   cellId=-1;
1419   ret0=std::numeric_limits<double>::max();
1420   for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1421     {
1422       switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1423       {
1424         case INTERP_KERNEL::NORM_SEG2:
1425           {
1426             std::size_t uselessEntry=0;
1427             double tmp=INTERP_KERNEL::SquareDistanceFromPtToSegInSpaceDim2(pt,coords+2*nc[ncI[*zeCell]+1],coords+2*nc[ncI[*zeCell]+2],uselessEntry);
1428             tmp=sqrt(tmp);
1429             if(tmp<ret0)
1430               { ret0=tmp; cellId=*zeCell; }
1431             break;
1432           }
1433         default:
1434           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint2DCurveAlg : not managed cell type ! Supporting SEG2 !");
1435       }
1436     }
1437 }
1438 DataArrayIdType *MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg(std::vector<bool>& fetched, const mcIdType *seedBg, const mcIdType *seedEnd, const DataArrayIdType *arrIn, const DataArrayIdType *arrIndxIn, mcIdType nbOfDepthPeeling, mcIdType& nbOfDepthPeelingPerformed)
1439 {
1440   nbOfDepthPeelingPerformed=0;
1441   if(!seedBg || !seedEnd || !arrIn || !arrIndxIn)
1442     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg : some input pointer is NULL !");
1443   mcIdType nbOfTuples=arrIndxIn->getNumberOfTuples()-1;
1444   std::vector<bool> fetched2(nbOfTuples,false);
1445   int i=0;
1446   for(const mcIdType *seedElt=seedBg;seedElt!=seedEnd;seedElt++,i++)
1447     {
1448       if(*seedElt>=0 && *seedElt<nbOfTuples)
1449         { fetched[*seedElt]=true; fetched2[*seedElt]=true; }
1450       else
1451         { 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()); }
1452     }
1453   const mcIdType *arrInPtr=arrIn->begin();
1454   const mcIdType *arrIndxPtr=arrIndxIn->begin();
1455   mcIdType targetNbOfDepthPeeling=nbOfDepthPeeling!=-1?nbOfDepthPeeling:std::numeric_limits<mcIdType>::max();
1456   std::vector<mcIdType> idsToFetch1(seedBg,seedEnd);
1457   std::vector<mcIdType> idsToFetch2;
1458   std::vector<mcIdType> *idsToFetch=&idsToFetch1;
1459   std::vector<mcIdType> *idsToFetchOther=&idsToFetch2;
1460   while(!idsToFetch->empty() && nbOfDepthPeelingPerformed<targetNbOfDepthPeeling)
1461     {
1462       for(std::vector<mcIdType>::const_iterator it=idsToFetch->begin();it!=idsToFetch->end();it++)
1463         for(const mcIdType *it2=arrInPtr+arrIndxPtr[*it];it2!=arrInPtr+arrIndxPtr[*it+1];it2++)
1464           if(!fetched[*it2])
1465             { fetched[*it2]=true; fetched2[*it2]=true; idsToFetchOther->push_back(*it2); }
1466       std::swap(idsToFetch,idsToFetchOther);
1467       idsToFetchOther->clear();
1468       nbOfDepthPeelingPerformed++;
1469     }
1470   mcIdType lgth=ToIdType(std::count(fetched2.begin(),fetched2.end(),true));
1471   i=0;
1472   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(lgth,1);
1473   mcIdType *retPtr=ret->getPointer();
1474   for(std::vector<bool>::const_iterator it=fetched2.begin();it!=fetched2.end();it++,i++)
1475     if(*it)
1476       *retPtr++=i;
1477   return ret.retn();
1478 }
1479
1480 /*!
1481  * This method put in zip format into parameter 'zipFrmt' in full interlace mode.
1482  * This format is often asked by INTERP_KERNEL algorithms to avoid many indirections into coordinates array.
1483  */
1484 void MEDCouplingUMesh::FillInCompact3DMode(int spaceDim, mcIdType nbOfNodesInCell, const mcIdType *conn, const double *coo, double *zipFrmt)
1485 {
1486   double *w=zipFrmt;
1487   if(spaceDim==3)
1488     for(int i=0;i<nbOfNodesInCell;i++)
1489       w=std::copy(coo+3*conn[i],coo+3*conn[i]+3,w);
1490   else if(spaceDim==2)
1491     {
1492       for(int i=0;i<nbOfNodesInCell;i++)
1493         {
1494           w=std::copy(coo+2*conn[i],coo+2*conn[i]+2,w);
1495           *w++=0.;
1496         }
1497     }
1498   else
1499     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::FillInCompact3DMode : Invalid spaceDim specified : must be 2 or 3 !");
1500 }
1501
1502 /*!
1503  * 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.
1504  * \param nbOfNodesPerLev in parameter that specifies the number of nodes of one slice of global dataset
1505  * \param isQuad specifies the policy of connectivity.
1506  * @ret in/out parameter in which the result will be append
1507  */
1508 void MEDCouplingUMesh::AppendExtrudedCell(const mcIdType *connBg, const mcIdType *connEnd, mcIdType nbOfNodesPerLev, bool isQuad, std::vector<mcIdType>& ret)
1509 {
1510   INTERP_KERNEL::NormalizedCellType flatType=(INTERP_KERNEL::NormalizedCellType)connBg[0];
1511   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(flatType);
1512   ret.push_back(cm.getExtrudedType());
1513   mcIdType deltaz=isQuad?2*nbOfNodesPerLev:nbOfNodesPerLev;
1514   switch(flatType)
1515   {
1516     case INTERP_KERNEL::NORM_POINT1:
1517       {
1518         ret.push_back(connBg[1]);
1519         ret.push_back(connBg[1]+nbOfNodesPerLev);
1520         break;
1521       }
1522     case INTERP_KERNEL::NORM_SEG2:
1523       {
1524         mcIdType conn[4]={connBg[1],connBg[2],connBg[2]+deltaz,connBg[1]+deltaz};
1525         ret.insert(ret.end(),conn,conn+4);
1526         break;
1527       }
1528     case INTERP_KERNEL::NORM_SEG3:
1529       {
1530         mcIdType conn[8]={connBg[1],connBg[3],connBg[3]+deltaz,connBg[1]+deltaz,connBg[2],connBg[3]+nbOfNodesPerLev,connBg[2]+deltaz,connBg[1]+nbOfNodesPerLev};
1531         ret.insert(ret.end(),conn,conn+8);
1532         break;
1533       }
1534     case INTERP_KERNEL::NORM_QUAD4:
1535       {
1536         mcIdType conn[8]={connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz};
1537         ret.insert(ret.end(),conn,conn+8);
1538         break;
1539       }
1540     case INTERP_KERNEL::NORM_TRI3:
1541       {
1542         mcIdType conn[6]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz};
1543         ret.insert(ret.end(),conn,conn+6);
1544         break;
1545       }
1546     case INTERP_KERNEL::NORM_TRI6:
1547       {
1548         mcIdType 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,
1549           connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
1550         ret.insert(ret.end(),conn,conn+15);
1551         break;
1552       }
1553     case INTERP_KERNEL::NORM_QUAD8:
1554       {
1555         mcIdType conn[20]={
1556           connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz,
1557           connBg[5],connBg[6],connBg[7],connBg[8],connBg[5]+deltaz,connBg[6]+deltaz,connBg[7]+deltaz,connBg[8]+deltaz,
1558           connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev,connBg[4]+nbOfNodesPerLev
1559         };
1560         ret.insert(ret.end(),conn,conn+20);
1561         break;
1562       }
1563     case INTERP_KERNEL::NORM_POLYGON:
1564       {
1565         std::back_insert_iterator< std::vector<mcIdType> > ii(ret);
1566         std::copy(connBg+1,connEnd,ii);
1567         *ii++=-1;
1568         std::reverse_iterator<const mcIdType *> rConnBg(connEnd);
1569         std::reverse_iterator<const mcIdType *> rConnEnd(connBg+1);
1570         std::transform(rConnBg,rConnEnd,ii,std::bind2nd(std::plus<mcIdType>(),deltaz));
1571         std::size_t nbOfRadFaces=std::distance(connBg+1,connEnd);
1572         for(std::size_t i=0;i<nbOfRadFaces;i++)
1573           {
1574             *ii++=-1;
1575             mcIdType conn[4]={connBg[(i+1)%nbOfRadFaces+1],connBg[i+1],connBg[i+1]+deltaz,connBg[(i+1)%nbOfRadFaces+1]+deltaz};
1576             std::copy(conn,conn+4,ii);
1577           }
1578         break;
1579       }
1580     default:
1581       throw INTERP_KERNEL::Exception("A flat type has been detected that has not its extruded representation !");
1582   }
1583 }
1584
1585
1586 /*!
1587  * 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).
1588  * 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
1589  * with a plane. The result will be put in 'cut3DSuf' out parameter.
1590  * \param [in] cut3DCurve  input parameter that gives for each 3DCurve cell if it owns fully to the plane or partially.
1591  * \param [out] nodesOnPlane, returns all the nodes that are on the plane.
1592  * \param [in] nodal3DSurf is the nodal connectivity of 3D surf mesh.
1593  * \param [in] nodalIndx3DSurf is the nodal connectivity index of 3D surf mesh.
1594  * \param [in] nodal3DCurve is the nodal connectivity of 3D curve mesh.
1595  * \param [in] nodal3DIndxCurve is the nodal connectivity index of 3D curve mesh.
1596  * \param [in] desc is the descending connectivity 3DSurf->3DCurve
1597  * \param [in] descIndx is the descending connectivity index 3DSurf->3DCurve
1598  * \param [out] cut3DSuf input/output param.
1599  */
1600 void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<mcIdType>& cut3DCurve, std::vector<mcIdType>& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf,
1601                                                    const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve,
1602                                                    const mcIdType *desc, const mcIdType *descIndx,
1603                                                    std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf)
1604 {
1605   std::set<mcIdType> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
1606   mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size());
1607   for(mcIdType i=0;i<nbOf3DSurfCell;i++)
1608     {
1609       std::vector<mcIdType> res;
1610       mcIdType offset=descIndx[i];
1611       mcIdType nbOfSeg=descIndx[i+1]-offset;
1612       for(mcIdType j=0;j<nbOfSeg;j++)
1613         {
1614           mcIdType edgeId=desc[offset+j];
1615           mcIdType status=cut3DCurve[edgeId];
1616           if(status!=-2)
1617             {
1618               if(status>-1)
1619                 res.push_back(status);
1620               else
1621                 {
1622                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
1623                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
1624                 }
1625             }
1626         }
1627       switch(res.size())
1628       {
1629         case 2:
1630           {
1631             cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1632             break;
1633           }
1634         case 1:
1635         case 0:
1636           {
1637             std::set<mcIdType> s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]);
1638             std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector<mcIdType> >(res));
1639             if(res.size()==2)
1640               {
1641                 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1642               }
1643             else
1644               {
1645                 cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
1646               }
1647             break;
1648           }
1649         default:
1650           {// case when plane is on a multi colinear edge of a polyhedron
1651             if(ToIdType(res.size())==2*nbOfSeg)
1652               {
1653                 cut3DSurf[i].first=-2; cut3DSurf[i].second=i;
1654               }
1655             else
1656               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
1657           }
1658       }
1659     }
1660 }
1661
1662
1663 /*!
1664  * \a this is expected to be a mesh with spaceDim==3 and meshDim==3. If not an exception will be thrown.
1665  * 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).
1666  * 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
1667  * with a plane. The result will be put in 'nodalRes' 'nodalResIndx' and 'cellIds' out parameters.
1668  * \param cut3DSurf  input parameter that gives for each 3DSurf its intersection with plane (result of MEDCouplingUMesh::AssemblyForSplitFrom3DCurve).
1669  * \param desc is the descending connectivity 3D->3DSurf
1670  * \param descIndx is the descending connectivity index 3D->3DSurf
1671  */
1672 void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
1673                                                   const mcIdType *desc, const mcIdType *descIndx,
1674                                                   DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const
1675 {
1676   checkFullyDefined();
1677   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1678     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1679   const mcIdType *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1680   mcIdType nbOfCells=getNumberOfCells();
1681   for(mcIdType i=0;i<nbOfCells;i++)
1682     {
1683       std::map<mcIdType, std::set<mcIdType> > m;
1684       mcIdType offset=descIndx[i];
1685       mcIdType nbOfFaces=descIndx[i+1]-offset;
1686       mcIdType start=-1;
1687       mcIdType end=-1;
1688       for(int j=0;j<nbOfFaces;j++)
1689         {
1690           const std::pair<mcIdType,mcIdType>& p=cut3DSurf[desc[offset+j]];
1691           if(p.first!=-1 && p.second!=-1)
1692             {
1693               if(p.first!=-2)
1694                 {
1695                   start=p.first; end=p.second;
1696                   m[p.first].insert(p.second);
1697                   m[p.second].insert(p.first);
1698                 }
1699               else
1700                 {
1701                   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]);
1702                   mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1703                   INTERP_KERNEL::AutoPtr<mcIdType> tmp=new mcIdType[sz];
1704                   INTERP_KERNEL::NormalizedCellType cmsId;
1705                   unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId);
1706                   start=tmp[0]; end=tmp[nbOfNodesSon-1];
1707                   for(unsigned k=0;k<nbOfNodesSon;k++)
1708                     {
1709                       m[tmp[k]].insert(tmp[(k+1)%nbOfNodesSon]);
1710                       m[tmp[(k+1)%nbOfNodesSon]].insert(tmp[k]);
1711                     }
1712                 }
1713             }
1714         }
1715       if(m.empty())
1716         continue;
1717       std::vector<mcIdType> conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON));
1718       mcIdType prev=end;
1719       while(end!=start)
1720         {
1721           std::map<mcIdType, std::set<mcIdType> >::const_iterator it=m.find(start);
1722           const std::set<mcIdType>& s=(*it).second;
1723           std::set<mcIdType> s2; s2.insert(prev);
1724           std::set<mcIdType> s3;
1725           std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
1726           if(s3.size()==1)
1727             {
1728               mcIdType val=*s3.begin();
1729               conn.push_back(start);
1730               prev=start;
1731               start=val;
1732             }
1733           else
1734             start=end;
1735         }
1736       conn.push_back(end);
1737       if(conn.size()>3)
1738         {
1739           nodalRes->insertAtTheEnd(conn.begin(),conn.end());
1740           nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples());
1741           cellIds->pushBackSilent(i);
1742         }
1743     }
1744 }
1745
1746
1747 void MEDCouplingUMesh::ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayIdType *nodalConnec, const DataArrayIdType *nodalConnecIndex)
1748 {
1749   if(nodalConnec && nodalConnecIndex)
1750     {
1751       types.clear();
1752       const mcIdType *conn(nodalConnec->begin()),*connIndex(nodalConnecIndex->begin());
1753       mcIdType nbOfElem=ToIdType(nodalConnecIndex->getNbOfElems())-1;
1754       if(nbOfElem>0)
1755         for(const mcIdType *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
1756           types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
1757     }
1758 }
1759
1760 /*!
1761  * This method expects that \a this a quadratic 1D, 2D or 3D mesh.
1762  * This method will 'attract' middle points of seg3 (deduced from this by explosion if needed) of cells connected to nodes specified in [\a nodeIdsBg, \a nodeIdsEnd )
1763  * For those selected mid points, their coordinates will be modified by applying a dilation between node in input [\a nodeIdsBg, \a nodeIdsEnd ) and the corresponding mid points using \a ratio input value.
1764  * So this method is non const because coordinates are modified.
1765  * If there is a couple of 2 points in [\a nodeIdsBg, \a nodeIdsEnd ) that are boundaries of a seg3, the corresponding mid point will remain untouched.
1766  *
1767  * \param [in] ratio - ratio of dilation
1768  * \param [in] nodeIdsBg - start (included) of input node list
1769  * \param [in] nodeIdsEnd - end (excluded) of input node list
1770  * \throw if there is a point in [\a nodeIdsBg, \a nodeIdsEnd ) that is a mid point of a seg3
1771  * \warning in case of throw the coordinates may be partially modified before the exception arises
1772  */
1773 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodes(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1774 {
1775   checkFullyDefined();
1776   int mdim(getMeshDimension());
1777   if(mdim==2 || mdim==3)
1778     {
1779       MCAuto<MEDCouplingUMesh> edges;
1780       {
1781         MCAuto<DataArrayIdType> a,b,c,d;
1782         edges=this->explodeIntoEdges(a,b,c,d);
1783       }
1784       return edges->attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1785     }
1786   if(mdim==1)
1787     return attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1788   throw INTERP_KERNEL::Exception("MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : not managed dimension ! Should be in [1,2,3] !");
1789 }
1790
1791 /*!
1792  * \a this is expected to have meshdim==1.
1793  */
1794 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodesUnderground(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1795 {
1796   int spaceDim(getSpaceDimension());
1797   double *coords(getCoords()->getPointer());
1798   auto nbNodes(getNumberOfNodes());
1799   std::size_t nbCells(getNumberOfCells());
1800   std::vector<bool> fastFinder(nbNodes,false);
1801   for(auto work=nodeIdsBg;work!=nodeIdsEnd;work++)
1802     if(*work>=0 && *work<nbNodes)
1803       fastFinder[*work]=true;
1804   MCAuto<DataArrayIdType> cellsIds(getCellIdsLyingOnNodes(nodeIdsBg,nodeIdsEnd,false));
1805   const mcIdType *nc(_nodal_connec->begin()),*nci(_nodal_connec_index->begin());
1806   for(std::size_t cellId=0;cellId<nbCells;cellId++,nci++)
1807     {
1808       const mcIdType *isSelected(std::find_if(nc+nci[0]+1,nc+nci[1],[&fastFinder](mcIdType v) { return fastFinder[v]; }));
1809       if(isSelected!=nc+nci[1])
1810         {
1811           if((INTERP_KERNEL::NormalizedCellType)nc[nci[0]]==INTERP_KERNEL::NORM_SEG3 && nci[1]-nci[0]==4)
1812             {
1813               bool aa(fastFinder[nc[*nci+1]]),bb(fastFinder[nc[*nci+2]]),cc(fastFinder[nc[*nci+3]]);
1814               if(!cc)
1815                 {
1816                   if(aa^bb)
1817                     {
1818                       auto ptToMove(nc[*nci+3]);
1819                       auto attractor(aa?nc[*nci+1]:nc[*nci+2]),endPt(aa?nc[*nci+2]:nc[*nci+1]);
1820                       std::transform(coords+spaceDim*attractor,coords+spaceDim*(attractor+1),coords+spaceDim*endPt,
1821                                      coords+spaceDim*ptToMove,[ratio](const double& stPt, const double& endPt) { return stPt+ratio*(endPt-stPt); });
1822                     }
1823                   else
1824                     continue;//both 2 boundary nodes of current seg3 are un nodeIds input list -> skip it.
1825                 }
1826               else
1827                 {
1828                   std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " has a mid point " << nc[*nci+3] << " ! This node is in input list !";
1829                   throw INTERP_KERNEL::Exception(oss.str());
1830                 }
1831             }
1832           else
1833             {
1834               std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " sharing one of the input nodes list its geo type is NOT SEG3 !";
1835               throw INTERP_KERNEL::Exception(oss.str());
1836             }
1837         }
1838     }
1839 }