src/MEDCoupling/MEDCouplingUMesh_internal.cxx

   1 // Copyright (C) 2007-2021  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   if(getNodalConnectivity()->empty())
 514   {
 515     cellIdsKeptArr=cellIdsKept.retn();
 516     return;
 517   }
 518   mcIdType sz=getNodalConnectivity()->getMaxValue(tmp); sz=std::max(sz,ToIdType(0))+1;
 519   std::vector<bool> fastFinder(sz,false);
 520   for(const mcIdType *work=begin;work!=end;work++)
 521     if(*work>=0 && *work<sz)
 522       fastFinder[*work]=true;
 523   mcIdType nbOfCells=getNumberOfCells();
 524   const mcIdType *conn=getNodalConnectivity()->getConstPointer();
 525   const mcIdType *connIndex=getNodalConnectivityIndex()->getConstPointer();
 526   for(mcIdType i=0;i<nbOfCells;i++)
 527     {
 528       mcIdType ref=0,nbOfHit=0;
 529       for(const mcIdType *work2=conn+connIndex[i]+1;work2!=conn+connIndex[i+1];work2++)
 530         if(*work2>=0)
 531           {
 532             ref++;
 533             if(fastFinder[*work2])
 534               nbOfHit++;
 535           }
 536       if((ref==nbOfHit && fullyIn) || (nbOfHit!=0 && !fullyIn))
 537         cellIdsKept->pushBackSilent(i);
 538     }
 539   cellIdsKeptArr=cellIdsKept.retn();
 540 }
 541
 542 /*!
 543  * This method works on a 3D curve linear mesh that is to say (meshDim==1 and spaceDim==3).
 544  * If it is not the case an exception will be thrown.
 545  * This method is non const because the coordinate of \a this can be appended with some new points issued from
 546  * intersection of plane defined by ('origin','vec').
 547  * This method has one in/out parameter : 'cut3DCurve'.
 548  * Param 'cut3DCurve' is expected to be of size 'this->getNumberOfCells()'. For each i in [0,'this->getNumberOfCells()')
 549  * if cut3DCurve[i]==-2, it means that for cell #i in \a this nothing has been detected previously.
 550  * if cut3DCurve[i]==-1, it means that cell#i has been already detected to be fully part of plane defined by ('origin','vec').
 551  * This method will throw an exception if \a this contains a non linear segment.
 552  */
 553 void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double *vec, double eps, std::vector<mcIdType>& cut3DCurve)
 554 {
 555   checkFullyDefined();
 556   if(getMeshDimension()!=1 || getSpaceDimension()!=3)
 557     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane works on umeshes with meshdim equal to 1 and spaceDim equal to 3 !");
 558   mcIdType ncells=getNumberOfCells();
 559   mcIdType nnodes=getNumberOfNodes();
 560   double vec2[3],vec3[3],vec4[3];
 561   double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
 562   if(normm<1e-6)
 563     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : parameter 'vec' should have a norm2 greater than 1e-6 !");
 564   vec2[0]=vec[0]/normm; vec2[1]=vec[1]/normm; vec2[2]=vec[2]/normm;
 565   const mcIdType *conn=_nodal_connec->getConstPointer();
 566   const mcIdType *connI=_nodal_connec_index->getConstPointer();
 567   const double *coo=_coords->getConstPointer();
 568   std::vector<double> addCoo;
 569   for(mcIdType i=0;i<ncells;i++)
 570     {
 571       if(conn[connI[i]]==ToIdType(INTERP_KERNEL::NORM_SEG2))
 572         {
 573           if(cut3DCurve[i]==-2)
 574             {
 575               mcIdType st=conn[connI[i]+1],endd=conn[connI[i]+2];
 576               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];
 577               double normm2=sqrt(vec3[0]*vec3[0]+vec3[1]*vec3[1]+vec3[2]*vec3[2]);
 578               double colin=std::abs((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2])/normm2);
 579               if(colin>eps)//if colin<=eps -> current SEG2 is colinear to the input plane
 580                 {
 581                   const double *st2=coo+3*st;
 582                   vec4[0]=st2[0]-origin[0]; vec4[1]=st2[1]-origin[1]; vec4[2]=st2[2]-origin[2];
 583                   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]));
 584                   if(pos>eps && pos<1-eps)
 585                     {
 586                       mcIdType nNode=ToIdType(addCoo.size())/3;
 587                       vec4[0]=st2[0]+pos*vec3[0]; vec4[1]=st2[1]+pos*vec3[1]; vec4[2]=st2[2]+pos*vec3[2];
 588                       addCoo.insert(addCoo.end(),vec4,vec4+3);
 589                       cut3DCurve[i]=nnodes+nNode;
 590                     }
 591                 }
 592             }
 593         }
 594       else
 595         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : this method is only available for linear cell (NORM_SEG2) !");
 596     }
 597   if(!addCoo.empty())
 598     {
 599       mcIdType newNbOfNodes=nnodes+ToIdType(addCoo.size())/3;
 600       MCAuto<DataArrayDouble> coo2=DataArrayDouble::New();
 601       coo2->alloc(newNbOfNodes,3);
 602       double *tmp=coo2->getPointer();
 603       tmp=std::copy(_coords->begin(),_coords->end(),tmp);
 604       std::copy(addCoo.begin(),addCoo.end(),tmp);
 605       DataArrayDouble::SetArrayIn(coo2,_coords);
 606     }
 607 }
 608
 609 /*!
 610  * This method incarnates the policy 0 for MEDCouplingUMesh::buildExtrudedMesh method.
 611  * \param mesh1D is the input 1D mesh used for translation computation.
 612  * \return newCoords new coords filled by this method.
 613  */
 614 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
 615 {
 616   mcIdType oldNbOfNodes=getNumberOfNodes();
 617   mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
 618   std::size_t spaceDim=getSpaceDimension();
 619   DataArrayDouble *ret=DataArrayDouble::New();
 620   std::vector<bool> isQuads;
 621   mcIdType nbOfLevsInVec=isQuad?2*nbOf1DCells+1:nbOf1DCells+1;
 622   ret->alloc(oldNbOfNodes*nbOfLevsInVec,spaceDim);
 623   double *retPtr=ret->getPointer();
 624   const double *coords=getCoords()->getConstPointer();
 625   double *work=std::copy(coords,coords+spaceDim*oldNbOfNodes,retPtr);
 626   std::vector<mcIdType> v;
 627   std::vector<double> c;
 628   double vec[3];
 629   v.reserve(3);
 630   c.reserve(6);
 631   for(mcIdType i=0;i<nbOf1DCells;i++)
 632     {
 633       v.resize(0);
 634       mesh1D->getNodeIdsOfCell(i,v);
 635       c.resize(0);
 636       mesh1D->getCoordinatesOfNode(v[isQuad?2:1],c);
 637       mesh1D->getCoordinatesOfNode(v[0],c);
 638       std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
 639       for(mcIdType j=0;j<oldNbOfNodes;j++)
 640         work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
 641       if(isQuad)
 642         {
 643           c.resize(0);
 644           mesh1D->getCoordinatesOfNode(v[1],c);
 645           mesh1D->getCoordinatesOfNode(v[0],c);
 646           std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
 647           for(mcIdType j=0;j<oldNbOfNodes;j++)
 648             work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
 649         }
 650     }
 651   ret->copyStringInfoFrom(*getCoords());
 652   return ret;
 653 }
 654
 655 /*!
 656  * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
 657  * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
 658  * \return newCoords new coords filled by this method.
 659  */
 660 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
 661 {
 662   if(mesh1D->getSpaceDimension()==2)
 663     return fillExtCoordsUsingTranslAndAutoRotation2D(mesh1D,isQuad);
 664   if(mesh1D->getSpaceDimension()==3)
 665     return fillExtCoordsUsingTranslAndAutoRotation3D(mesh1D,isQuad);
 666   throw INTERP_KERNEL::Exception("Not implemented rotation and translation alg. for spacedim other than 2 and 3 !");
 667 }
 668
 669 /*!
 670  * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
 671  * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
 672  * \return newCoords new coords filled by this method.
 673  */
 674 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
 675 {
 676   if(isQuad)
 677     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : not implemented for quadratic cells !");
 678   mcIdType oldNbOfNodes=getNumberOfNodes();
 679   mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
 680   if(nbOf1DCells<2)
 681     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
 682   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
 683   mcIdType nbOfLevsInVec=nbOf1DCells+1;
 684   ret->alloc(oldNbOfNodes*nbOfLevsInVec,2);
 685   double *retPtr=ret->getPointer();
 686   retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
 687   MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
 688   MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
 689   tmp->setCoords(tmp2);
 690   const double *coo1D=mesh1D->getCoords()->getConstPointer();
 691   const mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
 692   const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
 693   for(mcIdType i=1;i<nbOfLevsInVec;i++)
 694     {
 695       const double *begin=coo1D+2*conn1D[connI1D[i-1]+1];
 696       const double *end=coo1D+2*conn1D[connI1D[i-1]+2];
 697       const double *third=i+1<nbOfLevsInVec?coo1D+2*conn1D[connI1D[i]+2]:coo1D+2*conn1D[connI1D[i-2]+1];
 698       const double vec[2]={end[0]-begin[0],end[1]-begin[1]};
 699       tmp->translate(vec);
 700       double tmp3[2],radius,alpha,alpha0;
 701       const double *p0=i+1<nbOfLevsInVec?begin:third;
 702       const double *p1=i+1<nbOfLevsInVec?end:begin;
 703       const double *p2=i+1<nbOfLevsInVec?third:end;
 704       INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
 705       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]);
 706       double angle=acos(cosangle/(radius*radius));
 707       tmp->rotate(end,0,angle);
 708       retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
 709     }
 710   return ret.retn();
 711 }
 712
 713 /*!
 714  * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMesh method.
 715  * \param mesh1D is the input 1D mesh used for translation and automatic rotation computation.
 716  * \return newCoords new coords filled by this method.
 717  */
 718 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const
 719 {
 720   if(isQuad)
 721     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : not implemented for quadratic cells !");
 722   mcIdType oldNbOfNodes=getNumberOfNodes();
 723   mcIdType nbOf1DCells=ToIdType(mesh1D->getNumberOfCells());
 724   if(nbOf1DCells<2)
 725     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
 726   MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
 727   mcIdType nbOfLevsInVec=nbOf1DCells+1;
 728   ret->alloc(oldNbOfNodes*nbOfLevsInVec,3);
 729   double *retPtr=ret->getPointer();
 730   retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
 731   MCAuto<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
 732   MCAuto<DataArrayDouble> tmp2=getCoords()->deepCopy();
 733   tmp->setCoords(tmp2);
 734   const double *coo1D=mesh1D->getCoords()->getConstPointer();
 735   const mcIdType *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
 736   const mcIdType *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
 737   for(mcIdType i=1;i<nbOfLevsInVec;i++)
 738     {
 739       const double *begin=coo1D+3*conn1D[connI1D[i-1]+1];
 740       const double *end=coo1D+3*conn1D[connI1D[i-1]+2];
 741       const double *third=i+1<nbOfLevsInVec?coo1D+3*conn1D[connI1D[i]+2]:coo1D+3*conn1D[connI1D[i-2]+1];
 742       const double vec[3]={end[0]-begin[0],end[1]-begin[1],end[2]-begin[2]};
 743       tmp->translate(vec);
 744       double tmp3[2],radius,alpha,alpha0;
 745       const double *p0=i+1<nbOfLevsInVec?begin:third;
 746       const double *p1=i+1<nbOfLevsInVec?end:begin;
 747       const double *p2=i+1<nbOfLevsInVec?third:end;
 748       double vecPlane[3]={
 749         (p1[1]-p0[1])*(p2[2]-p1[2])-(p1[2]-p0[2])*(p2[1]-p1[1]),
 750         (p1[2]-p0[2])*(p2[0]-p1[0])-(p1[0]-p0[0])*(p2[2]-p1[2]),
 751         (p1[0]-p0[0])*(p2[1]-p1[1])-(p1[1]-p0[1])*(p2[0]-p1[0]),
 752       };
 753       double norm=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]+vecPlane[2]*vecPlane[2]);
 754       if(norm>1.e-7)
 755         {
 756           vecPlane[0]/=norm; vecPlane[1]/=norm; vecPlane[2]/=norm;
 757           double norm2=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]);
 758           double vec2[2]={vecPlane[1]/norm2,-vecPlane[0]/norm2};
 759           double s2=norm2;
 760           double c2=cos(asin(s2));
 761           double m[3][3]={
 762             {vec2[0]*vec2[0]*(1-c2)+c2, vec2[0]*vec2[1]*(1-c2), vec2[1]*s2},
 763             {vec2[0]*vec2[1]*(1-c2), vec2[1]*vec2[1]*(1-c2)+c2, -vec2[0]*s2},
 764             {-vec2[1]*s2, vec2[0]*s2, c2}
 765           };
 766           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]};
 767           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]};
 768           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]};
 769           INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0);
 770           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]);
 771           double angle=acos(cosangle/(radius*radius));
 772           tmp->rotate(end,vecPlane,angle);
 773         }
 774       retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
 775     }
 776   return ret.retn();
 777 }
 778
 779 /*!
 780  * This method is private because not easy to use for end user. This method is const contrary to
 781  * MEDCouplingUMesh::buildExtrudedMesh method because this->_coords are expected to contain
 782  * the coords sorted slice by slice.
 783  * \param isQuad specifies presence of quadratic cells.
 784  */
 785 MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(mcIdType nbOfNodesOf1Lev, bool isQuad) const
 786 {
 787   mcIdType nbOf1DCells(getNumberOfNodes()/nbOfNodesOf1Lev-1);
 788   mcIdType nbOf2DCells=getNumberOfCells();
 789   mcIdType nbOf3DCells(nbOf2DCells*nbOf1DCells);
 790   MEDCouplingUMesh *ret(MEDCouplingUMesh::New("Extruded",getMeshDimension()+1));
 791   const mcIdType *conn(_nodal_connec->begin()),*connI(_nodal_connec_index->begin());
 792   MCAuto<DataArrayIdType> newConn(DataArrayIdType::New()),newConnI(DataArrayIdType::New());
 793   newConnI->alloc(nbOf3DCells+1,1);
 794   mcIdType *newConnIPtr(newConnI->getPointer());
 795   *newConnIPtr++=0;
 796   std::vector<mcIdType> newc;
 797   for(mcIdType j=0;j<nbOf2DCells;j++)
 798     {
 799       AppendExtrudedCell(conn+connI[j],conn+connI[j+1],nbOfNodesOf1Lev,isQuad,newc);
 800       *newConnIPtr++=ToIdType(newc.size());
 801     }
 802   newConn->alloc(newc.size()*nbOf1DCells,1);
 803   mcIdType *newConnPtr(newConn->getPointer());
 804   mcIdType deltaPerLev(isQuad?2*nbOfNodesOf1Lev:nbOfNodesOf1Lev);
 805   newConnIPtr=newConnI->getPointer();
 806   for(mcIdType iz=0;iz<nbOf1DCells;iz++)
 807     {
 808       if(iz!=0)
 809         std::transform(newConnIPtr+1,newConnIPtr+1+nbOf2DCells,newConnIPtr+1+iz*nbOf2DCells,std::bind(std::plus<mcIdType>(),std::placeholders::_1,newConnIPtr[iz*nbOf2DCells]));
 810       const mcIdType *posOfTypeOfCell(newConnIPtr);
 811       for(std::vector<mcIdType>::const_iterator iter=newc.begin();iter!=newc.end();iter++,newConnPtr++)
 812         {
 813           mcIdType icell(ToIdType(iter-newc.begin()));//std::distance unfortunately cannot been called here in C++98
 814           if(icell!=*posOfTypeOfCell)
 815             {
 816               if(*iter!=-1)
 817                 *newConnPtr=(*iter)+iz*deltaPerLev;
 818               else
 819                 *newConnPtr=-1;
 820             }
 821           else
 822             {
 823               *newConnPtr=*iter;
 824               posOfTypeOfCell++;
 825             }
 826         }
 827     }
 828   ret->setConnectivity(newConn,newConnI,true);
 829   ret->setCoords(getCoords());
 830   return ret;
 831 }
 832
 833
 834 /*!
 835  * This method find in candidate pool defined by 'candidates' the cells equal following the polycy 'compType'.
 836  * If any true is returned and the results will be put at the end of 'result' output parameter. If not false is returned
 837  * and result remains unchanged.
 838  * The semantic of 'compType' is specified in MEDCouplingPointSet::zipConnectivityTraducer method.
 839  * If in 'candidates' pool -1 value is considered as an empty value.
 840  * WARNING this method returns only ONE set of result !
 841  */
 842 bool MEDCouplingUMesh::AreCellsEqualInPool(const std::vector<mcIdType>& candidates, int compType, const mcIdType *conn, const mcIdType *connI, DataArrayIdType *result)
 843 {
 844   if(candidates.size()<1)
 845     return false;
 846   bool ret=false;
 847   std::vector<mcIdType>::const_iterator iter=candidates.begin();
 848   mcIdType start=(*iter++);
 849   for(;iter!=candidates.end();iter++)
 850     {
 851       int status=AreCellsEqual(conn,connI,start,*iter,compType);
 852       if(status!=0)
 853         {
 854           if(!ret)
 855             {
 856               result->pushBackSilent(start);
 857               ret=true;
 858             }
 859           if(status==1)
 860             result->pushBackSilent(*iter);
 861           else
 862             result->pushBackSilent(status==2?(*iter+1):-(*iter+1));
 863         }
 864     }
 865   return ret;
 866 }
 867
 868 /*!
 869  * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelf.
 870  * Keeps from \a this only cells which constituing point id are in the ids specified by [ \a begin,\a end ).
 871  * The return newly allocated mesh will share the same coordinates as \a this.
 872  */
 873 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const mcIdType *begin, const mcIdType *end) const
 874 {
 875   checkConnectivityFullyDefined();
 876   mcIdType ncell=getNumberOfCells();
 877   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
 878   ret->_mesh_dim=_mesh_dim;
 879   ret->setCoords(_coords);
 880   std::size_t nbOfElemsRet=std::distance(begin,end);
 881   mcIdType *connIndexRet=(mcIdType *)malloc((nbOfElemsRet+1)*sizeof(mcIdType));
 882   connIndexRet[0]=0;
 883   const mcIdType *conn=_nodal_connec->getConstPointer();
 884   const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
 885   mcIdType newNbring=0;
 886   for(const mcIdType *work=begin;work!=end;work++,newNbring++)
 887     {
 888       if(*work>=0 && *work<ncell)
 889         connIndexRet[newNbring+1]=connIndexRet[newNbring]+connIndex[*work+1]-connIndex[*work];
 890       else
 891         {
 892           free(connIndexRet);
 893           std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoords : On pos #" << std::distance(begin,work) << " input cell id =" << *work << " should be in [0," << ncell << ") !";
 894           throw INTERP_KERNEL::Exception(oss.str());
 895         }
 896     }
 897   mcIdType *connRet=(mcIdType *)malloc(connIndexRet[nbOfElemsRet]*sizeof(mcIdType));
 898   mcIdType *connRetWork=connRet;
 899   std::set<INTERP_KERNEL::NormalizedCellType> types;
 900   for(const mcIdType *work=begin;work!=end;work++)
 901     {
 902       types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*work]]);
 903       connRetWork=std::copy(conn+connIndex[*work],conn+connIndex[*work+1],connRetWork);
 904     }
 905   MCAuto<DataArrayIdType> connRetArr=DataArrayIdType::New();
 906   connRetArr->useArray(connRet,true,DeallocType::C_DEALLOC,connIndexRet[nbOfElemsRet],1);
 907   MCAuto<DataArrayIdType> connIndexRetArr=DataArrayIdType::New();
 908   connIndexRetArr->useArray(connIndexRet,true,DeallocType::C_DEALLOC,nbOfElemsRet+1,1);
 909   ret->setConnectivity(connRetArr,connIndexRetArr,false);
 910   ret->_types=types;
 911   ret->copyTinyInfoFrom(this);
 912   return ret.retn();
 913 }
 914
 915 /*!
 916  * This is the low algorithm of MEDCouplingUMesh::buildPartOfMySelfSlice.
 917  * CellIds are given using range specified by a start an end and step.
 918  */
 919 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice(mcIdType start, mcIdType end, mcIdType step) const
 920 {
 921   checkFullyDefined();
 922   mcIdType ncell=getNumberOfCells();
 923   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
 924   ret->_mesh_dim=_mesh_dim;
 925   ret->setCoords(_coords);
 926   mcIdType newNbOfCells=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : ");
 927   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(newNbOfCells+1,1);
 928   mcIdType *newConnIPtr=newConnI->getPointer(); *newConnIPtr=0;
 929   mcIdType work=start;
 930   const mcIdType *conn=_nodal_connec->getConstPointer();
 931   const mcIdType *connIndex=_nodal_connec_index->getConstPointer();
 932   for(mcIdType i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
 933     {
 934       if(work>=0 && work<ncell)
 935         {
 936           newConnIPtr[1]=newConnIPtr[0]+connIndex[work+1]-connIndex[work];
 937         }
 938       else
 939         {
 940           std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoordsSlice : On pos #" << i << " input cell id =" << work << " should be in [0," << ncell << ") !";
 941           throw INTERP_KERNEL::Exception(oss.str());
 942         }
 943     }
 944   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(newConnIPtr[0],1);
 945   mcIdType *newConnPtr=newConn->getPointer();
 946   std::set<INTERP_KERNEL::NormalizedCellType> types;
 947   work=start;
 948   for(mcIdType i=0;i<newNbOfCells;i++,newConnIPtr++,work+=step)
 949     {
 950       types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[work]]);
 951       newConnPtr=std::copy(conn+connIndex[work],conn+connIndex[work+1],newConnPtr);
 952     }
 953   ret->setConnectivity(newConn,newConnI,false);
 954   ret->_types=types;
 955   ret->copyTinyInfoFrom(this);
 956   return ret.retn();
 957 }
 958
 959
 960 mcIdType MEDCouplingFastNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2)
 961 {
 962   return id;
 963 }
 964
 965 mcIdType MEDCouplingOrientationSensitiveNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2)
 966 {
 967   if(!compute)
 968     return id+1;
 969   else
 970     {
 971       if(cm.getOrientationStatus(nb,conn1,conn2))
 972         return id+1;
 973       else
 974         return -(id+1);
 975     }
 976 }
 977
 978
 979 /*!
 980  * Implements \a conversionType 0 for meshes with meshDim = 1, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
 981  * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
 982  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
 983  */
 984 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic1D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
 985 {
 986   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
 987   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
 988   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
 989   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
 990   mcIdType nbOfCells=getNumberOfCells();
 991   mcIdType nbOfNodes=getNumberOfNodes();
 992   const mcIdType *cPtr=_nodal_connec->begin();
 993   const mcIdType *icPtr=_nodal_connec_index->begin();
 994   mcIdType lastVal=0,offset=nbOfNodes;
 995   for(mcIdType i=0;i<nbOfCells;i++,icPtr++)
 996     {
 997       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
 998       if(type==INTERP_KERNEL::NORM_SEG2)
 999         {
1000           types.insert(INTERP_KERNEL::NORM_SEG3);
1001           newConn->pushBackSilent(ToIdType(INTERP_KERNEL::NORM_SEG3));
1002           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[0]+3);
1003           newConn->pushBackSilent(offset++);
1004           lastVal+=4;
1005           newConnI->pushBackSilent(lastVal);
1006           ret->pushBackSilent(i);
1007         }
1008       else
1009         {
1010           types.insert(type);
1011           lastVal+=(icPtr[1]-icPtr[0]);
1012           newConnI->pushBackSilent(lastVal);
1013           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1014         }
1015     }
1016   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1017   coords=DataArrayDouble::Aggregate(getCoords(),tmp); conn=newConn.retn(); connI=newConnI.retn();
1018   return ret.retn();
1019 }
1020
1021 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
1022 {
1023   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1024   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1025   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1026   //
1027   const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1028   DataArrayIdType *conn1D=0,*conn1DI=0;
1029   std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1030   DataArrayDouble *coordsTmp=0;
1031   MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1032   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1033   MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1034   const mcIdType *c1DPtr=conn1D->begin();
1035   const mcIdType *c1DIPtr=conn1DI->begin();
1036   mcIdType nbOfCells=getNumberOfCells();
1037   const mcIdType *cPtr=_nodal_connec->begin();
1038   const mcIdType *icPtr=_nodal_connec_index->begin();
1039   mcIdType lastVal=0;
1040   for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1041     {
1042       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1043       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1044       if(!cm.isQuadratic())
1045         {
1046           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType();
1047           types.insert(typ2); newConn->pushBackSilent(typ2);
1048           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1049           for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1050             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1051           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0]);
1052           newConnI->pushBackSilent(lastVal);
1053           ret->pushBackSilent(i);
1054         }
1055       else
1056         {
1057           types.insert(typ);
1058           lastVal+=(icPtr[1]-icPtr[0]);
1059           newConnI->pushBackSilent(lastVal);
1060           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1061         }
1062     }
1063   conn=newConn.retn(); connI=newConnI.retn(); coords=coordsTmpSafe.retn();
1064   return ret.retn();
1065 }
1066
1067 /*!
1068  * Implements \a conversionType 0 for meshes with meshDim = 2, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1069  * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1070  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1071  */
1072 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1073 {
1074   MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1075   MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1076   return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1077 }
1078
1079 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic2D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1080 {
1081   MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1082   MCAuto<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1083   //
1084   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1085   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1086   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(0,1);
1087   //
1088   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1089   const mcIdType *descPtr(desc->begin()),*descIPtr(descI->begin());
1090   DataArrayIdType *conn1D=0,*conn1DI=0;
1091   std::set<INTERP_KERNEL::NormalizedCellType> types1D;
1092   DataArrayDouble *coordsTmp=0;
1093   MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=0;
1094   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1095   MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1096   const mcIdType *c1DPtr=conn1D->begin();
1097   const mcIdType *c1DIPtr=conn1DI->begin();
1098   mcIdType nbOfCells=getNumberOfCells();
1099   const mcIdType *cPtr=_nodal_connec->begin();
1100   const mcIdType *icPtr=_nodal_connec_index->begin();
1101   mcIdType lastVal=0;
1102   mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1103   for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++)
1104     {
1105       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1106       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1107       if(!cm.isQuadratic())
1108         {
1109           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1110           types.insert(typ2); newConn->pushBackSilent(typ2);
1111           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1112           for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1113             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1114           newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1115           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+1;
1116           newConnI->pushBackSilent(lastVal);
1117           ret->pushBackSilent(i);
1118         }
1119       else
1120         {
1121           types.insert(typ);
1122           lastVal+=(icPtr[1]-icPtr[0]);
1123           newConnI->pushBackSilent(lastVal);
1124           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1125         }
1126     }
1127   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1128   coords=DataArrayDouble::Aggregate(coordsTmpSafe,tmp); conn=newConn.retn(); connI=newConnI.retn();
1129   return ret.retn();
1130 }
1131
1132 /*!
1133  * Implements \a conversionType 0 for meshes with meshDim = 3, of MEDCouplingUMesh::convertLinearCellsToQuadratic method.
1134  * \return a newly created DataArrayIdType instance that the caller should deal with containing cell ids of converted cells.
1135  * \sa MEDCouplingUMesh::convertLinearCellsToQuadratic.
1136  */
1137 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D0(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1138 {
1139   MCAuto<DataArrayIdType> desc(DataArrayIdType::New()),descI(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1140   MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
1141   return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
1142 }
1143
1144 DataArrayIdType *MEDCouplingUMesh::convertLinearCellsToQuadratic3D1(DataArrayIdType *&conn, DataArrayIdType *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
1145 {
1146   MCAuto<DataArrayIdType> desc2(DataArrayIdType::New()),desc2I(DataArrayIdType::New()),tmp2(DataArrayIdType::New()),tmp3(DataArrayIdType::New());
1147   MCAuto<MEDCouplingUMesh> m2D=buildDescendingConnectivityGen<MinusOneSonsGeneratorBiQuadratic>(desc2,desc2I,tmp2,tmp3,MEDCouplingFastNbrer); tmp2=0; tmp3=0;
1148   MCAuto<DataArrayIdType> desc1(DataArrayIdType::New()),desc1I(DataArrayIdType::New()),tmp4(DataArrayIdType::New()),tmp5(DataArrayIdType::New());
1149   MCAuto<MEDCouplingUMesh> m1D=explode3DMeshTo1D(desc1,desc1I,tmp4,tmp5); tmp4=0; tmp5=0;
1150   //
1151   MCAuto<DataArrayIdType> newConn=DataArrayIdType::New(); newConn->alloc(0,1);
1152   MCAuto<DataArrayIdType> newConnI=DataArrayIdType::New(); newConnI->alloc(1,1); newConnI->setIJ(0,0,0);
1153   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(),ret2=DataArrayIdType::New(); ret->alloc(0,1); ret2->alloc(0,1);
1154   //
1155   MCAuto<DataArrayDouble> bary=computeCellCenterOfMass();
1156   const mcIdType *descPtr(desc1->begin()),*descIPtr(desc1I->begin()),*desc2Ptr(desc2->begin()),*desc2IPtr(desc2I->begin());
1157   DataArrayIdType *conn1D=0,*conn1DI=0,*conn2D=0,*conn2DI=0;
1158   std::set<INTERP_KERNEL::NormalizedCellType> types1D,types2D;
1159   DataArrayDouble *coordsTmp=0,*coordsTmp2=0;
1160   MCAuto<DataArrayIdType> ret1D=m1D->convertLinearCellsToQuadratic1D0(conn1D,conn1DI,coordsTmp,types1D); ret1D=DataArrayIdType::New(); ret1D->alloc(0,1);
1161   MCAuto<DataArrayIdType> conn1DSafe(conn1D),conn1DISafe(conn1DI);
1162   MCAuto<DataArrayDouble> coordsTmpSafe(coordsTmp);
1163   MCAuto<DataArrayIdType> ret2D=m2D->convertLinearCellsToQuadratic2D1(conn2D,conn2DI,coordsTmp2,types2D); ret2D=DataArrayIdType::New(); ret2D->alloc(0,1);
1164   MCAuto<DataArrayDouble> coordsTmp2Safe(coordsTmp2);
1165   MCAuto<DataArrayIdType> conn2DSafe(conn2D),conn2DISafe(conn2DI);
1166   const mcIdType *c1DPtr=conn1D->begin(),*c1DIPtr=conn1DI->begin(),*c2DPtr=conn2D->begin(),*c2DIPtr=conn2DI->begin();
1167   mcIdType nbOfCells=getNumberOfCells();
1168   const mcIdType *cPtr=_nodal_connec->begin();
1169   const mcIdType *icPtr=_nodal_connec_index->begin();
1170   mcIdType lastVal=0;
1171   mcIdType offset=coordsTmpSafe->getNumberOfTuples();
1172   for(mcIdType i=0;i<nbOfCells;i++,icPtr++,descIPtr++,desc2IPtr++)
1173     {
1174       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)cPtr[*icPtr];
1175       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
1176       if(!cm.isQuadratic())
1177         {
1178           INTERP_KERNEL::NormalizedCellType typ2=cm.getQuadraticType2();
1179           if(typ2==INTERP_KERNEL::NORM_ERROR)
1180             {
1181               std::ostringstream oss; oss << "MEDCouplingUMesh::convertLinearCellsToQuadratic3D1 : On cell #" << i << " the linear cell type does not support advanced quadratization !";
1182               throw INTERP_KERNEL::Exception(oss.str());
1183             }
1184           types.insert(typ2); newConn->pushBackSilent(typ2);
1185           newConn->pushBackValsSilent(cPtr+icPtr[0]+1,cPtr+icPtr[1]);
1186           for(const mcIdType *d=descPtr+descIPtr[0];d!=descPtr+descIPtr[1];d++)
1187             newConn->pushBackSilent(c1DPtr[c1DIPtr[*d]+3]);
1188           for(const mcIdType *d=desc2Ptr+desc2IPtr[0];d!=desc2Ptr+desc2IPtr[1];d++)
1189             {
1190               mcIdType nodeId2=c2DPtr[c2DIPtr[(*d)+1]-1];
1191               mcIdType tmpPos=newConn->getNumberOfTuples();
1192               newConn->pushBackSilent(nodeId2);
1193               ret2D->pushBackSilent(nodeId2); ret1D->pushBackSilent(tmpPos);
1194             }
1195           newConn->pushBackSilent(offset+ret->getNumberOfTuples());
1196           lastVal+=(icPtr[1]-icPtr[0])+(descIPtr[1]-descIPtr[0])+(desc2IPtr[1]-desc2IPtr[0])+1;
1197           newConnI->pushBackSilent(lastVal);
1198           ret->pushBackSilent(i);
1199         }
1200       else
1201         {
1202           types.insert(typ);
1203           lastVal+=(icPtr[1]-icPtr[0]);
1204           newConnI->pushBackSilent(lastVal);
1205           newConn->pushBackValsSilent(cPtr+icPtr[0],cPtr+icPtr[1]);
1206         }
1207     }
1208   MCAuto<DataArrayIdType> diffRet2D=ret2D->getDifferentValues();
1209   MCAuto<DataArrayIdType> o2nRet2D=diffRet2D->invertArrayN2O2O2N(coordsTmp2Safe->getNumberOfTuples());
1210   coordsTmp2Safe=coordsTmp2Safe->selectByTupleId(diffRet2D->begin(),diffRet2D->end());
1211   MCAuto<DataArrayDouble> tmp=bary->selectByTupleIdSafe(ret->begin(),ret->end());
1212   std::vector<const DataArrayDouble *> v(3); v[0]=coordsTmpSafe; v[1]=coordsTmp2Safe; v[2]=tmp;
1213   mcIdType *c=newConn->getPointer();
1214   const mcIdType *cI(newConnI->begin());
1215   for(const mcIdType *elt=ret1D->begin();elt!=ret1D->end();elt++)
1216     c[*elt]=o2nRet2D->getIJ(c[*elt],0)+offset;
1217   offset=coordsTmp2Safe->getNumberOfTuples();
1218   for(const mcIdType *elt=ret->begin();elt!=ret->end();elt++)
1219     c[cI[(*elt)+1]-1]+=offset;
1220   coords=DataArrayDouble::Aggregate(v); conn=newConn.retn(); connI=newConnI.retn();
1221   return ret.retn();
1222 }
1223
1224 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshLinear(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1225 {
1226   mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1227   const mcIdType *n2oPtr(n2o->begin());
1228   MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1229   skin->getReverseNodalConnectivity(revNodal,revNodalI);
1230   const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1231   const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1232   const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1233   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1234   mcIdType *work(ret->getPointer());  *work++=INTERP_KERNEL::NORM_POLYGON;
1235   if(nbOfNodesExpected<1)
1236     return ret.retn();
1237   mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1238   *work++=n2oPtr[prevNode];
1239   for(mcIdType i=1;i<nbOfNodesExpected;i++)
1240     {
1241       if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==3)
1242         {
1243           std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1244           conn.erase(prevNode);
1245           if(conn.size()==1)
1246             {
1247               mcIdType curNode(*(conn.begin()));
1248               *work++=n2oPtr[curNode];
1249               std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1250               shar.erase(prevCell);
1251               if(shar.size()==1)
1252                 {
1253                   prevCell=*(shar.begin());
1254                   prevNode=curNode;
1255                 }
1256               else
1257                 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 2 !");
1258             }
1259           else
1260             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 1 !");
1261         }
1262       else
1263         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected cell !");
1264     }
1265   return ret.retn();
1266 }
1267
1268 DataArrayIdType *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayIdType *n2o) const
1269 {
1270   mcIdType nbOfNodesExpected(skin->getNumberOfNodes());
1271   mcIdType nbOfTurn(nbOfNodesExpected/2);
1272   const mcIdType *n2oPtr(n2o->begin());
1273   MCAuto<DataArrayIdType> revNodal(DataArrayIdType::New()),revNodalI(DataArrayIdType::New());
1274   skin->getReverseNodalConnectivity(revNodal,revNodalI);
1275   const mcIdType *revNodalPtr(revNodal->begin()),*revNodalIPtr(revNodalI->begin());
1276   const mcIdType *nodalPtr(skin->getNodalConnectivity()->begin());
1277   const mcIdType *nodalIPtr(skin->getNodalConnectivityIndex()->begin());
1278   MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfNodesExpected+1,1);
1279   mcIdType *work(ret->getPointer());  *work++=INTERP_KERNEL::NORM_QPOLYG;
1280   if(nbOfNodesExpected<1)
1281     return ret.retn();
1282   mcIdType prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
1283   *work=n2oPtr[prevNode]; work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[0]+3]]; work++;
1284   for(mcIdType i=1;i<nbOfTurn;i++)
1285     {
1286       if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==4)
1287         {
1288           std::set<mcIdType> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
1289           conn.erase(prevNode);
1290           if(conn.size()==1)
1291             {
1292               mcIdType curNode(*(conn.begin()));
1293               *work=n2oPtr[curNode];
1294               std::set<mcIdType> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
1295               shar.erase(prevCell);
1296               if(shar.size()==1)
1297                 {
1298                   mcIdType curCell(*(shar.begin()));
1299                   work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[curCell]+3]];
1300                   prevCell=curCell;
1301                   prevNode=curNode;
1302                   work++;
1303                 }
1304               else
1305                 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 2 !");
1306             }
1307           else
1308             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 1 !");
1309         }
1310       else
1311         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected cell !");
1312     }
1313   return ret.retn();
1314 }
1315
1316 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(const std::vector<const MEDCouplingUMesh *>& a)
1317 {
1318   if(a.empty())
1319     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : input array must be NON EMPTY !");
1320   std::vector<const MEDCouplingUMesh *>::const_iterator it=a.begin();
1321   int meshDim=(*it)->getMeshDimension();
1322   mcIdType nbOfCells=ToIdType((*it)->getNumberOfCells());
1323   mcIdType meshLgth=(*it++)->getNodalConnectivityArrayLen();
1324   for(;it!=a.end();it++)
1325     {
1326       if(meshDim!=(*it)->getMeshDimension())
1327         throw INTERP_KERNEL::Exception("Mesh dimensions mismatches, MergeUMeshes impossible !");
1328       nbOfCells+=ToIdType((*it)->getNumberOfCells());
1329       meshLgth+=(*it)->getNodalConnectivityArrayLen();
1330     }
1331   std::vector<const MEDCouplingPointSet *> aps(a.size());
1332   std::copy(a.begin(),a.end(),aps.begin());
1333   MCAuto<DataArrayDouble> pts=MergeNodesArray(aps);
1334   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("merge",meshDim);
1335   ret->setCoords(pts);
1336   MCAuto<DataArrayIdType> c=DataArrayIdType::New();
1337   c->alloc(meshLgth,1);
1338   mcIdType *cPtr=c->getPointer();
1339   MCAuto<DataArrayIdType> cI=DataArrayIdType::New();
1340   cI->alloc(nbOfCells+1,1);
1341   mcIdType *cIPtr=cI->getPointer();
1342   *cIPtr++=0;
1343   mcIdType offset=0;
1344   mcIdType offset2=0;
1345   for(it=a.begin();it!=a.end();it++)
1346     {
1347       mcIdType curNbOfCell=ToIdType((*it)->getNumberOfCells());
1348       const mcIdType *curCI=(*it)->_nodal_connec_index->begin();
1349       const mcIdType *curC=(*it)->_nodal_connec->begin();
1350       cIPtr=std::transform(curCI+1,curCI+curNbOfCell+1,cIPtr,std::bind(std::plus<mcIdType>(),std::placeholders::_1,offset));
1351       for(mcIdType j=0;j<curNbOfCell;j++)
1352         {
1353           const mcIdType *src=curC+curCI[j];
1354           *cPtr++=*src++;
1355           for(;src!=curC+curCI[j+1];src++,cPtr++)
1356             {
1357               if(*src!=-1)
1358                 *cPtr=*src+offset2;
1359               else
1360                 *cPtr=-1;
1361             }
1362         }
1363       offset+=curCI[curNbOfCell];
1364       offset2+=(*it)->getNumberOfNodes();
1365     }
1366   //
1367   ret->setConnectivity(c,cI,true);
1368   return ret.retn();
1369 }
1370
1371
1372 /*!
1373  * \param [in] pt the start pointer (included) of the coordinates of the point
1374  * \param [in] cellIdsBg the start pointer (included) of cellIds
1375  * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1376  * \param [in] nc nodal connectivity
1377  * \param [in] ncI nodal connectivity index
1378  * \param [in,out] ret0 the min distance between \a this and the external input point
1379  * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1380  * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1381  */
1382 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)
1383 {
1384   cellId=-1;
1385   ret0=std::numeric_limits<double>::max();
1386   for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1387     {
1388       switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1389       {
1390         case INTERP_KERNEL::NORM_TRI3:
1391           {
1392             double tmp=INTERP_KERNEL::DistanceFromPtToTriInSpaceDim3(pt,coords+3*nc[ncI[*zeCell]+1],coords+3*nc[ncI[*zeCell]+2],coords+3*nc[ncI[*zeCell]+3]);
1393             if(tmp<ret0)
1394               { ret0=tmp; cellId=*zeCell; }
1395             break;
1396           }
1397         case INTERP_KERNEL::NORM_QUAD4:
1398         case INTERP_KERNEL::NORM_POLYGON:
1399           {
1400             double tmp=INTERP_KERNEL::DistanceFromPtToPolygonInSpaceDim3(pt,nc+ncI[*zeCell]+1,nc+ncI[*zeCell+1],coords);
1401             if(tmp<ret0)
1402               { ret0=tmp; cellId=*zeCell; }
1403             break;
1404           }
1405         default:
1406           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint3DSurfAlg : not managed cell type ! Supporting TRI3, QUAD4 and POLYGON !");
1407       }
1408     }
1409 }
1410
1411 /*!
1412  * \param [in] pt the start pointer (included) of the coordinates of the point
1413  * \param [in] cellIdsBg the start pointer (included) of cellIds
1414  * \param [in] cellIdsEnd the end pointer (excluded) of cellIds
1415  * \param [in] nc nodal connectivity
1416  * \param [in] ncI nodal connectivity index
1417  * \param [in,out] ret0 the min distance between \a this and the external input point
1418  * \param [out] cellId that corresponds to minimal distance. If the closer node is not linked to any cell in \a this -1 is returned.
1419  * \sa MEDCouplingUMesh::distanceToPoint, MEDCouplingUMesh::distanceToPoints
1420  */
1421 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)
1422 {
1423   cellId=-1;
1424   ret0=std::numeric_limits<double>::max();
1425   for(const mcIdType *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
1426     {
1427       switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
1428       {
1429         case INTERP_KERNEL::NORM_SEG2:
1430           {
1431             std::size_t uselessEntry=0;
1432             double tmp=INTERP_KERNEL::SquareDistanceFromPtToSegInSpaceDim2(pt,coords+2*nc[ncI[*zeCell]+1],coords+2*nc[ncI[*zeCell]+2],uselessEntry);
1433             tmp=sqrt(tmp);
1434             if(tmp<ret0)
1435               { ret0=tmp; cellId=*zeCell; }
1436             break;
1437           }
1438         default:
1439           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint2DCurveAlg : not managed cell type ! Supporting SEG2 !");
1440       }
1441     }
1442 }
1443 DataArrayIdType *MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg(std::vector<bool>& fetched, const mcIdType *seedBg, const mcIdType *seedEnd, const DataArrayIdType *arrIn, const DataArrayIdType *arrIndxIn, mcIdType nbOfDepthPeeling, mcIdType& nbOfDepthPeelingPerformed)
1444 {
1445   nbOfDepthPeelingPerformed=0;
1446   if(!seedBg || !seedEnd || !arrIn || !arrIndxIn)
1447     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeedAlg : some input pointer is NULL !");
1448   mcIdType nbOfTuples=arrIndxIn->getNumberOfTuples()-1;
1449   std::vector<bool> fetched2(nbOfTuples,false);
1450   int i=0;
1451   for(const mcIdType *seedElt=seedBg;seedElt!=seedEnd;seedElt++,i++)
1452     {
1453       if(*seedElt>=0 && *seedElt<nbOfTuples)
1454         { fetched[*seedElt]=true; fetched2[*seedElt]=true; }
1455       else
1456         { 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()); }
1457     }
1458   const mcIdType *arrInPtr=arrIn->begin();
1459   const mcIdType *arrIndxPtr=arrIndxIn->begin();
1460   mcIdType targetNbOfDepthPeeling=nbOfDepthPeeling!=-1?nbOfDepthPeeling:std::numeric_limits<mcIdType>::max();
1461   std::vector<mcIdType> idsToFetch1(seedBg,seedEnd);
1462   std::vector<mcIdType> idsToFetch2;
1463   std::vector<mcIdType> *idsToFetch=&idsToFetch1;
1464   std::vector<mcIdType> *idsToFetchOther=&idsToFetch2;
1465   while(!idsToFetch->empty() && nbOfDepthPeelingPerformed<targetNbOfDepthPeeling)
1466     {
1467       for(std::vector<mcIdType>::const_iterator it=idsToFetch->begin();it!=idsToFetch->end();it++)
1468         for(const mcIdType *it2=arrInPtr+arrIndxPtr[*it];it2!=arrInPtr+arrIndxPtr[*it+1];it2++)
1469           if(!fetched[*it2])
1470             { fetched[*it2]=true; fetched2[*it2]=true; idsToFetchOther->push_back(*it2); }
1471       std::swap(idsToFetch,idsToFetchOther);
1472       idsToFetchOther->clear();
1473       nbOfDepthPeelingPerformed++;
1474     }
1475   mcIdType lgth=ToIdType(std::count(fetched2.begin(),fetched2.end(),true));
1476   i=0;
1477   MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(lgth,1);
1478   mcIdType *retPtr=ret->getPointer();
1479   for(std::vector<bool>::const_iterator it=fetched2.begin();it!=fetched2.end();it++,i++)
1480     if(*it)
1481       *retPtr++=i;
1482   return ret.retn();
1483 }
1484
1485 /*!
1486  * This method put in zip format into parameter 'zipFrmt' in full interlace mode.
1487  * This format is often asked by INTERP_KERNEL algorithms to avoid many indirections into coordinates array.
1488  */
1489 void MEDCouplingUMesh::FillInCompact3DMode(int spaceDim, mcIdType nbOfNodesInCell, const mcIdType *conn, const double *coo, double *zipFrmt)
1490 {
1491   double *w=zipFrmt;
1492   if(spaceDim==3)
1493     for(int i=0;i<nbOfNodesInCell;i++)
1494       w=std::copy(coo+3*conn[i],coo+3*conn[i]+3,w);
1495   else if(spaceDim==2)
1496     {
1497       for(int i=0;i<nbOfNodesInCell;i++)
1498         {
1499           w=std::copy(coo+2*conn[i],coo+2*conn[i]+2,w);
1500           *w++=0.;
1501         }
1502     }
1503   else
1504     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::FillInCompact3DMode : Invalid spaceDim specified : must be 2 or 3 !");
1505 }
1506
1507 /*!
1508  * 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.
1509  * \param nbOfNodesPerLev in parameter that specifies the number of nodes of one slice of global dataset
1510  * \param isQuad specifies the policy of connectivity.
1511  * @ret in/out parameter in which the result will be append
1512  */
1513 void MEDCouplingUMesh::AppendExtrudedCell(const mcIdType *connBg, const mcIdType *connEnd, mcIdType nbOfNodesPerLev, bool isQuad, std::vector<mcIdType>& ret)
1514 {
1515   INTERP_KERNEL::NormalizedCellType flatType=(INTERP_KERNEL::NormalizedCellType)connBg[0];
1516   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(flatType);
1517   ret.push_back(cm.getExtrudedType());
1518   mcIdType deltaz=isQuad?2*nbOfNodesPerLev:nbOfNodesPerLev;
1519   switch(flatType)
1520   {
1521     case INTERP_KERNEL::NORM_POINT1:
1522       {
1523         ret.push_back(connBg[1]);
1524         ret.push_back(connBg[1]+nbOfNodesPerLev);
1525         break;
1526       }
1527     case INTERP_KERNEL::NORM_SEG2:
1528       {
1529         mcIdType conn[4]={connBg[1],connBg[2],connBg[2]+deltaz,connBg[1]+deltaz};
1530         ret.insert(ret.end(),conn,conn+4);
1531         break;
1532       }
1533     case INTERP_KERNEL::NORM_SEG3:
1534       {
1535         mcIdType conn[8]={connBg[1],connBg[3],connBg[3]+deltaz,connBg[1]+deltaz,connBg[2],connBg[3]+nbOfNodesPerLev,connBg[2]+deltaz,connBg[1]+nbOfNodesPerLev};
1536         ret.insert(ret.end(),conn,conn+8);
1537         break;
1538       }
1539     case INTERP_KERNEL::NORM_QUAD4:
1540       {
1541         mcIdType conn[8]={connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz};
1542         ret.insert(ret.end(),conn,conn+8);
1543         break;
1544       }
1545     case INTERP_KERNEL::NORM_TRI3:
1546       {
1547         mcIdType conn[6]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz};
1548         ret.insert(ret.end(),conn,conn+6);
1549         break;
1550       }
1551     case INTERP_KERNEL::NORM_TRI6:
1552       {
1553         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,
1554           connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
1555         ret.insert(ret.end(),conn,conn+15);
1556         break;
1557       }
1558     case INTERP_KERNEL::NORM_QUAD8:
1559       {
1560         mcIdType conn[20]={
1561           connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz,
1562           connBg[5],connBg[6],connBg[7],connBg[8],connBg[5]+deltaz,connBg[6]+deltaz,connBg[7]+deltaz,connBg[8]+deltaz,
1563           connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev,connBg[4]+nbOfNodesPerLev
1564         };
1565         ret.insert(ret.end(),conn,conn+20);
1566         break;
1567       }
1568     case INTERP_KERNEL::NORM_POLYGON:
1569       {
1570         std::back_insert_iterator< std::vector<mcIdType> > ii(ret);
1571         std::copy(connBg+1,connEnd,ii);
1572         *ii++=-1;
1573         std::reverse_iterator<const mcIdType *> rConnBg(connEnd);
1574         std::reverse_iterator<const mcIdType *> rConnEnd(connBg+1);
1575         std::transform(rConnBg,rConnEnd,ii,std::bind(std::plus<mcIdType>(),std::placeholders::_1,deltaz));
1576         std::size_t nbOfRadFaces=std::distance(connBg+1,connEnd);
1577         for(std::size_t i=0;i<nbOfRadFaces;i++)
1578           {
1579             *ii++=-1;
1580             mcIdType conn[4]={connBg[(i+1)%nbOfRadFaces+1],connBg[i+1],connBg[i+1]+deltaz,connBg[(i+1)%nbOfRadFaces+1]+deltaz};
1581             std::copy(conn,conn+4,ii);
1582           }
1583         break;
1584       }
1585     default:
1586       throw INTERP_KERNEL::Exception("A flat type has been detected that has not its extruded representation !");
1587   }
1588 }
1589
1590
1591 /*!
1592  * 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).
1593  * 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
1594  * with a plane. The result will be put in 'cut3DSuf' out parameter.
1595  * \param [in] cut3DCurve  input parameter that gives for each 3DCurve cell if it owns fully to the plane or partially.
1596  * \param [out] nodesOnPlane returns all the nodes that are on the plane.
1597  * \param [in] nodal3DSurf is the nodal connectivity of 3D surf mesh.
1598  * \param [in] nodalIndx3DSurf is the nodal connectivity index of 3D surf mesh.
1599  * \param [in] nodal3DCurve is the nodal connectivity of 3D curve mesh.
1600  * \param [in] nodalIndx3DCurve is the nodal connectivity index of 3D curve mesh.
1601  * \param [in] desc is the descending connectivity 3DSurf->3DCurve
1602  * \param [in] descIndx is the descending connectivity index 3DSurf->3DCurve
1603  * \param [out] cut3DSurf input/output param.
1604  */
1605 void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<mcIdType>& cut3DCurve, std::vector<mcIdType>& nodesOnPlane, const mcIdType *nodal3DSurf, const mcIdType *nodalIndx3DSurf,
1606                                                    const mcIdType *nodal3DCurve, const mcIdType *nodalIndx3DCurve,
1607                                                    const mcIdType *desc, const mcIdType *descIndx,
1608                                                    std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf)
1609 {
1610   std::set<mcIdType> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
1611   mcIdType nbOf3DSurfCell=ToIdType(cut3DSurf.size());
1612   for(mcIdType i=0;i<nbOf3DSurfCell;i++)
1613     {
1614       std::vector<mcIdType> res;
1615       mcIdType offset=descIndx[i];
1616       mcIdType nbOfSeg=descIndx[i+1]-offset;
1617       for(mcIdType j=0;j<nbOfSeg;j++)
1618         {
1619           mcIdType edgeId=desc[offset+j];
1620           mcIdType status=cut3DCurve[edgeId];
1621           if(status!=-2)
1622             {
1623               if(status>-1)
1624                 res.push_back(status);
1625               else
1626                 {
1627                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
1628                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
1629                 }
1630             }
1631         }
1632       switch(res.size())
1633       {
1634         case 2:
1635           {
1636             cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1637             break;
1638           }
1639         case 1:
1640         case 0:
1641           {
1642             std::set<mcIdType> s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]);
1643             std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector<mcIdType> >(res));
1644             if(res.size()==2)
1645               {
1646                 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
1647               }
1648             else
1649               {
1650                 cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
1651               }
1652             break;
1653           }
1654         default:
1655           {// case when plane is on a multi colinear edge of a polyhedron
1656             if(ToIdType(res.size())==2*nbOfSeg)
1657               {
1658                 cut3DSurf[i].first=-2; cut3DSurf[i].second=i;
1659               }
1660             else
1661               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
1662           }
1663       }
1664     }
1665 }
1666
1667
1668 /*!
1669  * \a this is expected to be a mesh with spaceDim==3 and meshDim==3. If not an exception will be thrown.
1670  * 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).
1671  * 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
1672  * with a plane. The result will be put in 'nodalRes' 'nodalResIndx' and 'cellIds' out parameters.
1673  * \param cut3DSurf  input parameter that gives for each 3DSurf its intersection with plane (result of MEDCouplingUMesh::AssemblyForSplitFrom3DCurve).
1674  * \param desc is the descending connectivity 3D->3DSurf
1675  * \param descIndx is the descending connectivity index 3D->3DSurf
1676  */
1677 void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<mcIdType,mcIdType> >& cut3DSurf,
1678                                                   const mcIdType *desc, const mcIdType *descIndx,
1679                                                   DataArrayIdType *nodalRes, DataArrayIdType *nodalResIndx, DataArrayIdType *cellIds) const
1680 {
1681   checkFullyDefined();
1682   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
1683     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
1684   const mcIdType *nodal3D(_nodal_connec->begin()),*nodalIndx3D(_nodal_connec_index->begin());
1685   mcIdType nbOfCells=getNumberOfCells();
1686   for(mcIdType i=0;i<nbOfCells;i++)
1687     {
1688       std::map<mcIdType, std::set<mcIdType> > m;
1689       mcIdType offset=descIndx[i];
1690       mcIdType nbOfFaces=descIndx[i+1]-offset;
1691       mcIdType start=-1;
1692       mcIdType end=-1;
1693       for(int j=0;j<nbOfFaces;j++)
1694         {
1695           const std::pair<mcIdType,mcIdType>& p=cut3DSurf[desc[offset+j]];
1696           if(p.first!=-1 && p.second!=-1)
1697             {
1698               if(p.first!=-2)
1699                 {
1700                   start=p.first; end=p.second;
1701                   m[p.first].insert(p.second);
1702                   m[p.second].insert(p.first);
1703                 }
1704               else
1705                 {
1706                   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]);
1707                   mcIdType sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
1708                   INTERP_KERNEL::AutoPtr<mcIdType> tmp=new mcIdType[sz];
1709                   INTERP_KERNEL::NormalizedCellType cmsId;
1710                   unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId);
1711                   start=tmp[0]; end=tmp[nbOfNodesSon-1];
1712                   for(unsigned k=0;k<nbOfNodesSon;k++)
1713                     {
1714                       m[tmp[k]].insert(tmp[(k+1)%nbOfNodesSon]);
1715                       m[tmp[(k+1)%nbOfNodesSon]].insert(tmp[k]);
1716                     }
1717                 }
1718             }
1719         }
1720       if(m.empty())
1721         continue;
1722       std::vector<mcIdType> conn(1,ToIdType(INTERP_KERNEL::NORM_POLYGON));
1723       mcIdType prev=end;
1724       while(end!=start)
1725         {
1726           std::map<mcIdType, std::set<mcIdType> >::const_iterator it=m.find(start);
1727           const std::set<mcIdType>& s=(*it).second;
1728           std::set<mcIdType> s2; s2.insert(prev);
1729           std::set<mcIdType> s3;
1730           std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
1731           if(s3.size()==1)
1732             {
1733               mcIdType val=*s3.begin();
1734               conn.push_back(start);
1735               prev=start;
1736               start=val;
1737             }
1738           else
1739             start=end;
1740         }
1741       conn.push_back(end);
1742       if(conn.size()>3)
1743         {
1744           nodalRes->insertAtTheEnd(conn.begin(),conn.end());
1745           nodalResIndx->pushBackSilent(nodalRes->getNumberOfTuples());
1746           cellIds->pushBackSilent(i);
1747         }
1748     }
1749 }
1750
1751
1752 void MEDCouplingUMesh::ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayIdType *nodalConnec, const DataArrayIdType *nodalConnecIndex)
1753 {
1754   if(nodalConnec && nodalConnecIndex)
1755     {
1756       types.clear();
1757       const mcIdType *conn(nodalConnec->begin()),*connIndex(nodalConnecIndex->begin());
1758       mcIdType nbOfElem=ToIdType(nodalConnecIndex->getNbOfElems())-1;
1759       if(nbOfElem>0)
1760         for(const mcIdType *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
1761           types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
1762     }
1763 }
1764
1765 /*!
1766  * This method expects that \a this a quadratic 1D, 2D or 3D mesh.
1767  * 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 )
1768  * 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.
1769  * So this method is non const because coordinates are modified.
1770  * 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.
1771  *
1772  * \param [in] ratio - ratio of dilation
1773  * \param [in] nodeIdsBg - start (included) of input node list
1774  * \param [in] nodeIdsEnd - end (excluded) of input node list
1775  * \throw if there is a point in [\a nodeIdsBg, \a nodeIdsEnd ) that is a mid point of a seg3
1776  * \warning in case of throw the coordinates may be partially modified before the exception arises
1777  */
1778 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodes(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1779 {
1780   checkFullyDefined();
1781   int mdim(getMeshDimension());
1782   if(mdim==2 || mdim==3)
1783     {
1784       MCAuto<MEDCouplingUMesh> edges;
1785       {
1786         MCAuto<DataArrayIdType> a,b,c,d;
1787         edges=this->explodeIntoEdges(a,b,c,d);
1788       }
1789       return edges->attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1790     }
1791   if(mdim==1)
1792     return attractSeg3MidPtsAroundNodesUnderground(ratio,nodeIdsBg,nodeIdsEnd);
1793   throw INTERP_KERNEL::Exception("MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : not managed dimension ! Should be in [1,2,3] !");
1794 }
1795
1796 /*!
1797  * \a this is expected to have meshdim==1.
1798  */
1799 void MEDCouplingUMesh::attractSeg3MidPtsAroundNodesUnderground(double ratio, const mcIdType *nodeIdsBg, const mcIdType *nodeIdsEnd)
1800 {
1801   int spaceDim(getSpaceDimension());
1802   double *coords(getCoords()->getPointer());
1803   auto nbNodes(getNumberOfNodes());
1804   std::size_t nbCells(getNumberOfCells());
1805   std::vector<bool> fastFinder(nbNodes,false);
1806   for(auto work=nodeIdsBg;work!=nodeIdsEnd;work++)
1807     if(*work>=0 && *work<nbNodes)
1808       fastFinder[*work]=true;
1809   MCAuto<DataArrayIdType> cellsIds(getCellIdsLyingOnNodes(nodeIdsBg,nodeIdsEnd,false));
1810   const mcIdType *nc(_nodal_connec->begin()),*nci(_nodal_connec_index->begin());
1811   for(std::size_t cellId=0;cellId<nbCells;cellId++,nci++)
1812     {
1813       const mcIdType *isSelected(std::find_if(nc+nci[0]+1,nc+nci[1],[&fastFinder](mcIdType v) { return fastFinder[v]; }));
1814       if(isSelected!=nc+nci[1])
1815         {
1816           if((INTERP_KERNEL::NormalizedCellType)nc[nci[0]]==INTERP_KERNEL::NORM_SEG3 && nci[1]-nci[0]==4)
1817             {
1818               bool aa(fastFinder[nc[*nci+1]]),bb(fastFinder[nc[*nci+2]]),cc(fastFinder[nc[*nci+3]]);
1819               if(!cc)
1820                 {
1821                   if(aa^bb)
1822                     {
1823                       auto ptToMove(nc[*nci+3]);
1824                       auto attractor(aa?nc[*nci+1]:nc[*nci+2]),endPt(aa?nc[*nci+2]:nc[*nci+1]);
1825                       std::transform(coords+spaceDim*attractor,coords+spaceDim*(attractor+1),coords+spaceDim*endPt,
1826                                      coords+spaceDim*ptToMove,[ratio](const double& stPt2, const double& endPt2) { return stPt2+ratio*(endPt2-stPt2); });
1827                     }
1828                   else
1829                     continue;//both 2 boundary nodes of current seg3 are un nodeIds input list -> skip it.
1830                 }
1831               else
1832                 {
1833                   std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " has a mid point " << nc[*nci+3] << " ! This node is in input list !";
1834                   throw INTERP_KERNEL::Exception(oss.str());
1835                 }
1836             }
1837           else
1838             {
1839               std::ostringstream oss; oss << "MEDCouplingUMesh::attractSeg3MidPtsAroundNodes : cell #" << cellId << " sharing one of the input nodes list its geo type is NOT SEG3 !";
1840               throw INTERP_KERNEL::Exception(oss.str());
1841             }
1842         }
1843     }
1844 }