src/INTERP_KERNEL/SplitterTetra.txx

   1 // Copyright (C) 2007-2024  CEA, EDF
   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 #ifndef __SPLITTERTETRA_TXX__
  20 #define __SPLITTERTETRA_TXX__
  21
  22 #include "SplitterTetra.hxx"
  23
  24 #include "TetraAffineTransform.hxx"
  25 #include "TransformedTriangle.hxx"
  26 #include "MeshUtils.hxx"
  27 #include "VectorUtils.hxx"
  28 #include "CellModel.hxx"
  29 #include "Log.hxx"
  30 #include "UnitTetraIntersectionBary.hxx"
  31 #include "VolSurfFormulae.hxx"
  32
  33 #include <cmath>
  34 #include <cassert>
  35 #include <string>
  36 #include <sstream>
  37 #include <vector>
  38
  39 namespace INTERP_KERNEL
  40 {
  41   template<class MyMeshType>
  42   const double SplitterTetra<MyMeshType>::SPARSE_TRUNCATION_LIMIT=1.0e-14;
  43
  44   /*!
  45    * output is expected to be allocated with 24*sizeof(void*) in order to store the 24 tetras.
  46    * These tetras have to be deallocated.
  47    */
  48   template<class MyMeshType>
  49   void SplitterTetra<MyMeshType>::splitIntoDualCells(SplitterTetra<MyMeshType> **output)
  50   {
  51     double tmp[12];
  52     const double *tmp2[4]={tmp,tmp+3,tmp+6,tmp+9};
  53     typename MyMeshType::MyConnType conn[4]={-1,-1,-1,-1};
  54     for(int i=0;i<24;i++)
  55       {
  56         splitMySelfForDual(tmp,i,conn[0]);
  57         output[i]=new SplitterTetra<MyMeshType>(_src_mesh,tmp2,conn);
  58       }
  59   }
  60
  61   /**
  62    * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a
  63    * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
  64    * if a partial computation of a triangle has already been performed (to increase performance).
  65    *
  66    * The \a srcMesh can contain polyhedron cells.
  67    *
  68    *
  69    * Constructor creating object from the four corners of the tetrahedron.
  70    *
  71    * @param srcMesh       mesh containing the source elements
  72    * @param tetraCorners  array of four pointers to double[3] arrays containing the coordinates of the
  73    *                      corners of the tetrahedron
  74    */
  75   template<class MyMeshType>
  76   SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId)
  77     : _t(0), _src_mesh(srcMesh)
  78   {
  79     std::copy(nodesId,nodesId+4,_conn);
  80     _coords[0]=tetraCorners[0][0]; _coords[1]=tetraCorners[0][1]; _coords[2]=tetraCorners[0][2];
  81     _coords[3]=tetraCorners[1][0]; _coords[4]=tetraCorners[1][1]; _coords[5]=tetraCorners[1][2];
  82     _coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2];
  83     _coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2];
  84     // create the affine transform
  85     _t=new TetraAffineTransform(_coords);
  86   }
  87
  88   /**
  89    * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a
  90    * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
  91    * if a partial computation of a triangle has already been performed (to increase performance).
  92    *
  93    * The \a srcMesh can contain polyhedron cells.
  94    *
  95    *
  96    * Constructor creating object from the four corners of the tetrahedron.
  97    *
  98    * \param [in] srcMesh       mesh containing the source elements
  99    * \param [in] tetraCorners  array 4*3 doubles containing corners of input tetrahedron (P0X,P0Y,P0Y,P1X,P1Y,P1Z,P2X,P2Y,P2Z,P3X,P3Y,P3Z).
 100    */
 101   template<class MyMeshType>
 102   SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const ConnType *conn): _t(0),_src_mesh(srcMesh)
 103   {
 104     if(!conn)
 105       { _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3; }
 106     else
 107       { _conn[0]=conn[0]; _conn[1]=conn[1]; _conn[2]=conn[2]; _conn[3]=conn[3]; }
 108     _coords[0]=tetraCorners[0]; _coords[1]=tetraCorners[1]; _coords[2]=tetraCorners[2]; _coords[3]=tetraCorners[3]; _coords[4]=tetraCorners[4]; _coords[5]=tetraCorners[5];
 109     _coords[6]=tetraCorners[6]; _coords[7]=tetraCorners[7]; _coords[8]=tetraCorners[8]; _coords[9]=tetraCorners[9]; _coords[10]=tetraCorners[10]; _coords[11]=tetraCorners[11];
 110     // create the affine transform
 111     _t=new TetraAffineTransform(_coords);
 112   }
 113
 114   /**
 115    * Destructor
 116    *
 117    * Deletes _t and the coordinates (double[3] vectors) in _nodes
 118    *
 119    */
 120   template<class MyMeshType>
 121   SplitterTetra<MyMeshType>::~SplitterTetra()
 122   {
 123     delete _t;
 124     for(typename HashMap< ConnType, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
 125       delete[] iter->second;
 126   }
 127
 128   /*!
 129    * \Forget already calculated triangles, which is crucial for calculation of barycenter of intersection
 130    */
 131   template<class MyMeshType>
 132   void SplitterTetra<MyMeshType>::clearVolumesCache()
 133   {
 134     _volumes.clear();
 135   }
 136
 137   /*!
 138    * This method destroys the 4 pointers pointed by tetraCorners[0],tetraCorners[1],tetraCorners[2] and tetraCorners[3]
 139    * @param i is in 0..23 included.
 140    * @param output is expected to be sized of 12 in order to.
 141    */
 142   template<class MyMeshType>
 143   void SplitterTetra<MyMeshType>::splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId)
 144   {
 145     double *tmp[4];
 146     int offset=i/6;
 147     nodeId=_conn[offset];
 148     tmp[0]=_coords+3*offset; tmp[1]=_coords+((offset+1)%4)*3; tmp[2]=_coords+((offset+2)%4)*3; tmp[3]=_coords+((offset+3)%4)*3;
 149     int caseToTreat=i%6;
 150     int case1=caseToTreat/2;
 151     int case2=caseToTreat%2;
 152     const int tab[3][2]={{1,2},{3,2},{1,3}};
 153     const int *curTab=tab[case1];
 154     double pt0[3]; pt0[0]=(tmp[curTab[case2]][0]+tmp[0][0])/2.; pt0[1]=(tmp[curTab[case2]][1]+tmp[0][1])/2.; pt0[2]=(tmp[curTab[case2]][2]+tmp[0][2])/2.;
 155     double pt1[3]; pt1[0]=(tmp[0][0]+tmp[curTab[0]][0]+tmp[curTab[1]][0])/3.; pt1[1]=(tmp[0][1]+tmp[curTab[0]][1]+tmp[curTab[1]][1])/3.; pt1[2]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.;
 156     double pt2[3]; pt2[0]=(tmp[0][0]+tmp[1][0]+tmp[2][0]+tmp[3][0])/4.; pt2[1]=(tmp[0][1]+tmp[1][1]+tmp[2][1]+tmp[3][1])/4.; pt2[2]=(tmp[0][2]+tmp[1][2]+tmp[2][2]+tmp[3][2])/4.;
 157     std::copy(pt1,pt1+3,output+case2*3);
 158     std::copy(pt0,pt0+3,output+(abs(case2-1))*3);
 159     std::copy(pt2,pt2+3,output+2*3);
 160     std::copy(tmp[0],tmp[0]+3,output+3*3);
 161   }
 162
 163   /**
 164    * Calculates the volume of intersection of an element in the source mesh and the target element.
 165    * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
 166    * faces of the source element are triangulated and the calculated transformation is applied
 167    * to each triangle. The algorithm of Grandy, implemented in INTERP_KERNEL::TransformedTriangle is used
 168    * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
 169    * divided by the determinant of the transformation.
 170    *
 171    * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated
 172    * with triangulated faces to avoid having to recalculate these.
 173    *
 174    * @param element      global number of the source element in C mode.
 175    */
 176   template<class MyMeshType>
 177   double SplitterTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element,
 178                                                         double*                         baryCentre)
 179   {
 180     typedef typename MyMeshType::MyConnType ConnType;
 181     const NumberingPolicy numPol=MyMeshType::My_numPol;
 182     //{ could be done on outside?
 183     // check if we have planar tetra element
 184     if(_t->determinant() == 0.0)
 185       {
 186         // tetra is planar
 187         LOG(2, "Planar tetra -- volume 0");
 188         return 0.0;
 189       }
 190
 191     // get type of cell
 192     NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
 193     const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
 194     ConnType nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
 195     // halfspace filtering
 196     bool isOutside[8] = {true, true, true, true, true, true, true, true};
 197     bool isTargetOutside = false;
 198
 199     // calculate the coordinates of the nodes
 200     ConnType *cellNodes=new ConnType[nbOfNodes4Type];
 201     for(ConnType i = 0;i<nbOfNodes4Type;++i)
 202       {
 203         // we could store mapping local -> global numbers too, but not sure it is worth it
 204         const ConnType globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
 205         cellNodes[i]=globalNodeNum;
 206         if(_nodes.find(globalNodeNum) == _nodes.end())
 207           {
 208             //for(HashMap< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
 209             //  std::cout << (*iter3).first << " ";
 210             //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
 211             calculateNode(globalNodeNum);
 212           }
 213         CheckIsOutside(_nodes[globalNodeNum], isOutside);
 214       }
 215
 216     // halfspace filtering check
 217     // NB : might not be beneficial for caching of triangles
 218     for(int i = 0; i < 8; ++i)
 219       if(isOutside[i])
 220         isTargetOutside = true;
 221
 222     double totalVolume = 0.0;
 223
 224     if(!isTargetOutside)
 225       {
 226         /// calculator of intersection barycentre
 227         UnitTetraIntersectionBary baryCalculator( _t->determinant() < 0.);
 228
 229         // get nb of sons of a cell
 230         const ConnType* rawCellConn = _src_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _src_mesh.getConnectivityIndexPtr()[ element ]);
 231         const ConnType rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
 232         unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
 233
 234         for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
 235           {
 236             // get sons connectivity
 237             NormalizedCellType faceType;
 238             ConnType *faceNodes, nbFaceNodes=-1;
 239             if ( cellModelCell.isDynamic() )
 240               {
 241                 faceNodes=new ConnType[nbOfNodes4Type];
 242                 nbFaceNodes = cellModelCell.fillSonCellNodalConnectivity2(ii,rawCellConn,rawNbCellNodes,faceNodes,faceType);
 243                 for ( ConnType i = 0; i < nbFaceNodes; ++i )
 244                   faceNodes[i] = OTT<ConnType,numPol>::coo2C(faceNodes[i]);
 245               }
 246             else
 247               {
 248                 faceType = cellModelCell.getSonType(ii);
 249                 assert(CellModel::GetCellModel(faceType).getDimension() == 2);
 250                 nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon(ii);
 251                 faceNodes = new ConnType[nbFaceNodes];
 252                 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
 253               }
 254             // intersect a son with the unit tetra
 255             switch(faceType)
 256               {
 257               case NORM_TRI3:
 258                 {
 259                   // create the face key
 260                   TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
 261
 262                   // calculate the triangle if needed
 263                   if(_volumes.find(key) == _volumes.end())
 264                     {
 265                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
 266                       calculateVolume(tri, key);
 267                       totalVolume += _volumes[key];
 268                       if ( baryCentre )
 269                         baryCalculator.addSide( tri );
 270                     } else {
 271                       // count negative as face has reversed orientation
 272                       totalVolume -= _volumes[key];
 273                     }
 274                 }
 275                 break;
 276
 277               case NORM_QUAD4:
 278
 279                 // simple triangulation of faces along a diagonal :
 280                 //
 281                 // 2 ------ 3
 282                 // |      / |
 283                 // |     /  |
 284                 // |    /   |
 285                 // |   /    |
 286                 // |  /     |
 287                 // | /      |
 288                 // 1 ------ 4
 289                 //
 290                 //? not sure if this always works
 291                 {
 292                   // calculate the triangles if needed
 293
 294                   // local nodes 1, 2, 3
 295                   TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
 296                   if(_volumes.find(key1) == _volumes.end())
 297                     {
 298                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
 299                       calculateVolume(tri, key1);
 300                       totalVolume += _volumes[key1];
 301                     } else {
 302                       // count negative as face has reversed orientation
 303                       totalVolume -= _volumes[key1];
 304                     }
 305
 306                   // local nodes 1, 3, 4
 307                   TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]);
 308                   if(_volumes.find(key2) == _volumes.end())
 309                     {
 310                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]);
 311                       calculateVolume(tri, key2);
 312                       totalVolume += _volumes[key2];
 313                     }
 314                   else
 315                     {
 316                       // count negative as face has reversed orientation
 317                       totalVolume -= _volumes[key2];
 318                     }
 319                 }
 320                 break;
 321
 322               case NORM_POLYGON:
 323                 {
 324                   ConnType nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
 325                   for ( ConnType iTri = 0; iTri < nbTria; ++iTri )
 326                     {
 327                       TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]);
 328                       if(_volumes.find(key) == _volumes.end())
 329                         {
 330                           TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1+iTri]], _nodes[faceNodes[2+iTri]]);
 331                           calculateVolume(tri, key);
 332                           totalVolume += _volumes[key];
 333                         }
 334                       else
 335                         {
 336                           totalVolume -= _volumes[key];
 337                         }
 338                     }
 339                 }
 340                 break;
 341
 342               default:
 343                 std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
 344                 assert(false);
 345               }
 346             delete [] faceNodes;
 347           }
 348
 349         if ( baryCentre ) {
 350           baryCalculator.getBary( baryCentre );
 351           _t->reverseApply( baryCentre, baryCentre );
 352         }
 353       }
 354     delete [] cellNodes;
 355     // reset if it is very small to keep the matrix sparse
 356     // is this a good idea?
 357     if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
 358       {
 359         totalVolume = 0.0;
 360       }
 361
 362     LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
 363
 364     // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
 365     // that should be used (which is equivalent to dividing by the determinant)
 366     return std::fabs(1.0 / _t->determinant() * totalVolume) ;
 367   }
 368
 369   /**
 370    * Calculates the intersection surface of two coplanar triangles.
 371    *
 372    * @param planeNormal normal of the plane for the first triangle
 373    * @param planeConstant constant of the equation of the plane for the first triangle
 374    * @param p1 coordinates of the first  node of the first  triangle
 375    * @param p2 coordinates of the second node of the first  triangle
 376    * @param p3 coordinates of the third  node of the first  triangle
 377    * @param p4 coordinates of the first  node of the second triangle
 378    * @param p5 coordinates of the second node of the second triangle
 379    * @param p6 coordinates of the third  node of the second triangle
 380    * @param dimCaracteristic characteristic size of the meshes containing the triangles
 381    * @param precision precision for double float data used for comparison
 382    */
 383   template<class MyMeshType>
 384   double SplitterTetra<MyMeshType>::CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
 385                                                                                     const double planeConstant,
 386                                                                                     const double *const p1, const double *const p2, const double *const p3,
 387                                                                                     const double *const p4, const double *const p5, const double *const p6,
 388                                                                                     const double dimCaracteristic, const double precision)
 389   {
 390     typedef typename MyMeshType::MyConnType ConnType;
 391     typedef double Vect2[2];
 392     typedef double Triangle2[3][2];
 393
 394     const double *const tri0[3] = {p1, p2, p3};
 395     const double *const tri1[3] = {p4, p5, p6};
 396
 397     // Plane of the first triangle defined by the normal of the triangle and the constant
 398     // Project triangles onto coordinate plane most aligned with plane normal
 399     int maxNormal = 0;
 400     double fmax = std::abs(planeNormal[0]);
 401     double absMax = std::abs(planeNormal[1]);
 402     if (absMax > fmax)
 403       {
 404         maxNormal = 1;
 405         fmax = absMax;
 406       }
 407     absMax = std::abs(planeNormal[2]);
 408     if (absMax > fmax)
 409       {
 410         maxNormal = 2;
 411       }
 412
 413     Triangle2 projTri0, projTri1;
 414     int i;
 415
 416     if (maxNormal == 0)
 417       {
 418         // Project onto yz-plane.
 419         for (i = 0; i < 3; ++i)
 420           {
 421             projTri0[i][0] = tri0[i][1];
 422             projTri0[i][1] = tri0[i][2];
 423             projTri1[i][0] = tri1[i][1];
 424             projTri1[i][1] = tri1[i][2];
 425           }
 426       }
 427     else if (maxNormal == 1)
 428       {
 429         // Project onto xz-plane.
 430         for (i = 0; i < 3; ++i)
 431           {
 432             projTri0[i][0] = tri0[i][0];
 433             projTri0[i][1] = tri0[i][2];
 434             projTri1[i][0] = tri1[i][0];
 435             projTri1[i][1] = tri1[i][2];
 436           }
 437       }
 438     else
 439       {
 440         // Project onto xy-plane.
 441         for (i = 0; i < 3; ++i)
 442           {
 443             projTri0[i][0] = tri0[i][0];
 444             projTri0[i][1] = tri0[i][1];
 445             projTri1[i][0] = tri1[i][0];
 446             projTri1[i][1] = tri1[i][1];
 447           }
 448       }
 449
 450     // 2D triangle intersection routines require counterclockwise ordering.
 451     Vect2 save;
 452     Vect2 edge0;
 453     Vect2 edge1;
 454     for (int ii = 0; ii < 2; ++ii)
 455       {
 456         edge0[ii] = projTri0[1][ii] - projTri0[0][ii];
 457         edge1[ii] = projTri0[2][ii] - projTri0[0][ii];
 458       }
 459     if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
 460       {
 461         // Triangle is clockwise, reorder it.
 462         for (int ii = 0; ii < 2; ++ii)
 463           {
 464             save[ii] = projTri0[1][ii];
 465             projTri0[1][ii] = projTri0[2][ii];
 466             projTri0[2][ii] = save[ii];
 467           }
 468       }
 469
 470     for (int ii = 0; ii < 2; ++ii)
 471       {
 472         edge0[ii] = projTri1[1][ii] - projTri1[0][ii];
 473         edge1[ii] = projTri1[2][ii] - projTri1[0][ii];
 474       }
 475     if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
 476       {
 477         // Triangle is clockwise, reorder it.
 478       for (int ii = 0; ii < 2; ++ii)
 479         {
 480           save[ii] = projTri1[1][ii];
 481           projTri1[1][ii] = projTri1[2][ii];
 482           projTri1[2][ii] = save[ii];
 483         }
 484       }
 485
 486     std::vector<double> inter2;
 487     intersec_de_triangle(projTri0[0], projTri0[1], projTri0[2],
 488                          projTri1[0], projTri1[1], projTri1[2],
 489                          inter2,
 490                          dimCaracteristic, precision);
 491     ConnType nb_inter=((ConnType)inter2.size())/2;
 492     double surface = 0.;
 493     if(nb_inter >3) inter2=reconstruct_polygon(inter2);
 494     if (nb_inter > 0)
 495       {
 496         std::vector<double> inter3;
 497         inter3.resize(3 * nb_inter);
 498         // Map 2D intersections back to the 3D triangle space.
 499         if (maxNormal == 0)
 500           {
 501             double invNX = ((double) 1.) / planeNormal[0];
 502             for (i = 0; i < nb_inter; i++)
 503               {
 504                 inter3[3 * i + 1] = inter2[2 * i];
 505                 inter3[3 * i + 2] = inter2[2 * i + 1];
 506                 inter3[3 * i] = invNX * (planeConstant - planeNormal[1] * inter3[3 * i + 1] - planeNormal[2] * inter3[3 * i + 2]);
 507               }
 508           }
 509         else if (maxNormal == 1)
 510           {
 511             double invNY = ((double) 1.) / planeNormal[1];
 512             for (i = 0; i < nb_inter; i++)
 513               {
 514                 inter3[3 * i] = inter2[2 * i];
 515                 inter3[3 * i + 2] = inter2[2 * i + 1];
 516                 inter3[3 * i + 1] = invNY * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[2] * inter3[3 * i + 2]);
 517               }
 518           }
 519         else
 520           {
 521             double invNZ = ((double) 1.) / planeNormal[2];
 522             for (i = 0; i < nb_inter; i++)
 523               {
 524                 inter3[3 * i] = inter2[2 * i];
 525                 inter3[3 * i + 1] = inter2[2 * i + 1];
 526                 inter3[3 * i + 2] = invNZ * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[1] * inter3[3 * i + 1]);
 527               }
 528           }
 529         surface = polygon_area<3>(inter3);
 530       }
 531     return surface;
 532   }
 533
 534   /**
 535    * Determine if a face is coplanar with a triangle.
 536    * The first face is characterized by the equation of her plane
 537    *
 538    * @param planeNormal normal of the plane for the first triangle
 539    * @param planeConstant constant of the equation of the plane for the first triangle
 540    * @param coordsFace coordinates of the triangle face
 541    * @param precision precision for double float data used for comparison
 542    */
 543   template<class MyMeshType>
 544   bool SplitterTetra<MyMeshType>::IsFacesCoplanar(const double *const planeNormal,
 545                                                   const double planeConstant,
 546                                                   const double *const *const coordsFace,
 547                                                   const double precision)
 548   {
 549       // Compute the signed distances of triangle vertices to the plane. Use an epsilon-thick plane test.
 550       // For faces not left
 551       int counter = 0;
 552       for (int i = 0; i < 3; ++i)
 553         {
 554           const double distance = dot(planeNormal, coordsFace[i]) - planeConstant;
 555           if (epsilonEqual(distance, precision))
 556             {
 557               counter++;
 558             }
 559         }
 560       return counter == 3;
 561   }
 562
 563   /**
 564    * Calculates the surface of intersection of a polygon face in the source mesh and a cell of the target mesh.
 565    * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
 566    * faces of the source element are triangulated and the calculated transformation is applied
 567    * to each triangle.
 568    * The algorithm is based on the algorithm of Grandy used in intersectSourceCell to compute
 569    * the volume of intersection of two cell elements.
 570    * The case with a source face colinear to one of the face of tetrahedrons is taking into account:
 571    * the contribution of the face must not be counted two times.
 572    *
 573    * The class will cache the intermediary calculations of transformed nodes of source faces and surfaces associated
 574    * with triangulated faces to avoid having to recalculate these.
 575    *
 576    * @param polyType type of the polygon source face
 577    * @param polyNodesNbr number of the nodes of the polygon source face
 578    * @param polyNodes numbers of the nodes of the polygon source face
 579    * @param polyCoords coordinates of the nodes of the polygon source face
 580    * @param dimCaracteristic characteristic size of the meshes containing the triangles
 581    * @param precision precision for double float data used for comparison
 582    * @param listOfTetraFacesTreated list of tetra faces treated
 583    * @param listOfTetraFacesColinear list of tetra faces colinear with the polygon source faces
 584    */
 585   template<class MyMeshType>
 586   double SplitterTetra<MyMeshType>::intersectSourceFace(const NormalizedCellType polyType,
 587                                                         const ConnType polyNodesNbr,
 588                                                         const ConnType *const polyNodes,
 589                                                         const double *const *const polyCoords,
 590                                                         const double dimCaracteristic,
 591                                                         const double precision,
 592                                                         std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
 593                                                         std::set<TriangleFaceKey>& listOfTetraFacesColinear)
 594   {
 595     double totalSurface = 0.0;
 596
 597     // check if we have planar tetra element
 598     if(_t->determinant() == 0.0)
 599       {
 600         // tetra is planar
 601         LOG(2, "Planar tetra -- volume 0");
 602         return 0.0;
 603       }
 604
 605     // halfspace filtering
 606     bool isOutside[8] = {true, true, true, true, true, true, true, true};
 607     bool isStrictlyOutside[8] = {true, true, true, true, true, true, true, true};
 608     bool isTargetStrictlyOutside = false;
 609     bool isTargetOutside = false;
 610
 611     // calculate the coordinates of the nodes
 612     for(ConnType i = 0;i<polyNodesNbr;++i)
 613       {
 614         const ConnType globalNodeNum = polyNodes[i];
 615         if(_nodes.find(globalNodeNum) == _nodes.end())
 616           {
 617             calculateNode2(globalNodeNum, polyCoords[i]);
 618           }
 619
 620         CheckIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside, precision);
 621         CheckIsOutside(_nodes[globalNodeNum], isOutside, precision);
 622       }
 623
 624     // halfspace filtering check
 625     // NB : might not be beneficial for caching of triangles
 626     for(int i = 0; i < 8; ++i)
 627       {
 628         if(isStrictlyOutside[i])
 629           {
 630             isTargetStrictlyOutside = true;
 631             break;
 632           }
 633         else if (isOutside[i])
 634           {
 635             isTargetOutside = true;
 636           }
 637       }
 638
 639     if (!isTargetStrictlyOutside)
 640       {
 641
 642         if (isTargetOutside)
 643           {
 644             // Faces are parallel
 645             const int tetraFacesNodesConn[4][3] = {
 646                 { 0, 1, 2 },
 647                 { 0, 2, 3 },
 648                 { 0, 3, 1 },
 649                 { 1, 2, 3 } };
 650             double planeNormal[3];
 651             for (int iTetraFace = 0; iTetraFace < 4; ++iTetraFace)
 652               {
 653                 const int * const tetraFaceNodesConn = tetraFacesNodesConn[iTetraFace];
 654                 TriangleFaceKey key = TriangleFaceKey(_conn[tetraFaceNodesConn[0]],
 655                                                       _conn[tetraFaceNodesConn[1]],
 656                                                       _conn[tetraFaceNodesConn[2]]);
 657                 if (listOfTetraFacesTreated.find(key) == listOfTetraFacesTreated.end())
 658                   {
 659                     const double * const coordsTetraTriNode1 = _coords + tetraFaceNodesConn[0] * MyMeshType::MY_SPACEDIM;
 660                     const double * const coordsTetraTriNode2 = _coords + tetraFaceNodesConn[1] * MyMeshType::MY_SPACEDIM;
 661                     const double * const coordsTetraTriNode3 = _coords + tetraFaceNodesConn[2] * MyMeshType::MY_SPACEDIM;
 662                     calculateNormalForTria(coordsTetraTriNode1, coordsTetraTriNode2, coordsTetraTriNode3, planeNormal);
 663                     const double normOfTetraTriNormal = norm(planeNormal);
 664                     if (epsilonEqual(normOfTetraTriNormal, 0.))
 665                       {
 666                         for (int i = 0; i < 3; ++i)
 667                           {
 668                             planeNormal[i] = 0.;
 669                           }
 670                       }
 671                     else
 672                       {
 673                         const double invNormOfTetraTriNormal = 1. / normOfTetraTriNormal;
 674                         for (int i = 0; i < 3; ++i)
 675                           {
 676                             planeNormal[i] *= invNormOfTetraTriNormal;
 677                           }
 678                       }
 679                     double planeConstant = dot(planeNormal, coordsTetraTriNode1);
 680                     if (IsFacesCoplanar(planeNormal, planeConstant, polyCoords, precision))
 681                       {
 682                         ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
 683                         for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
 684                           {
 685                             double volume = CalculateIntersectionSurfaceOfCoplanarTriangles(planeNormal,
 686                                                                                             planeConstant,
 687                                                                                             polyCoords[0],
 688                                                                                             polyCoords[1 + iTri],
 689                                                                                             polyCoords[2 + iTri],
 690                                                                                             coordsTetraTriNode1,
 691                                                                                             coordsTetraTriNode2,
 692                                                                                             coordsTetraTriNode3,
 693                                                                                             dimCaracteristic,
 694                                                                                             precision);
 695                             if (!epsilonEqual(volume, 0.))
 696                               {
 697                                 totalSurface += volume;
 698                                 listOfTetraFacesColinear.insert(key);
 699                               }
 700                           }
 701                       }
 702                   }
 703                 listOfTetraFacesTreated.insert(key);
 704               }
 705           }
 706         else
 707           {
 708               // intersect a son with the unit tetra
 709               switch (polyType)
 710                 {
 711                 case NORM_TRI3:
 712                   {
 713                     // create the face key
 714                     TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]);
 715
 716                     // calculate the triangle if needed
 717                     if (_volumes.find(key) == _volumes.end())
 718                       {
 719                         TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
 720                         calculateSurface(tri, key);
 721                         totalSurface += _volumes[key];
 722                       }
 723                     else
 724                       {
 725                         // count negative as face has reversed orientation
 726                         totalSurface -= _volumes[key];
 727                       }
 728                   }
 729                   break;
 730
 731                 case NORM_QUAD4:
 732
 733                   // simple triangulation of faces along a diagonal :
 734                   //
 735                   // 2 ------ 3
 736                   // |      / |
 737                   // |     /  |
 738                   // |    /   |
 739                   // |   /    |
 740                   // |  /     |
 741                   // | /      |
 742                   // 1 ------ 4
 743                   //
 744                   //? not sure if this always works
 745                   {
 746                     // calculate the triangles if needed
 747
 748                     // local nodes 1, 2, 3
 749                     TriangleFaceKey key1 = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]);
 750                     if (_volumes.find(key1) == _volumes.end())
 751                       {
 752                         TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
 753                         calculateSurface(tri, key1);
 754                         totalSurface += _volumes[key1];
 755                       }
 756                     else
 757                       {
 758                         // count negative as face has reversed orientation
 759                         totalSurface -= _volumes[key1];
 760                       }
 761
 762                     // local nodes 1, 3, 4
 763                     TriangleFaceKey key2 = TriangleFaceKey(polyNodes[0], polyNodes[2], polyNodes[3]);
 764                     if (_volumes.find(key2) == _volumes.end())
 765                       {
 766                         TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[2]], _nodes[polyNodes[3]]);
 767                         calculateSurface(tri, key2);
 768                         totalSurface += _volumes[key2];
 769                       }
 770                     else
 771                       {
 772                         // count negative as face has reversed orientation
 773                         totalSurface -= _volumes[key2];
 774                       }
 775                   }
 776                   break;
 777
 778                 case NORM_POLYGON:
 779                   {
 780                     ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
 781                     for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
 782                       {
 783                         TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1 + iTri], polyNodes[2 + iTri]);
 784                         if (_volumes.find(key) == _volumes.end())
 785                           {
 786                             TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1 + iTri]],
 787                                 _nodes[polyNodes[2 + iTri]]);
 788                             calculateSurface(tri, key);
 789                             totalSurface += _volumes[key];
 790                           }
 791                         else
 792                           {
 793                             totalSurface -= _volumes[key];
 794                           }
 795                       }
 796                   }
 797                   break;
 798
 799                 default:
 800                   std::cout
 801                       << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment."
 802                       << std::endl;
 803                   assert(false);
 804                 }
 805
 806           }
 807       }
 808
 809     // reset if it is very small to keep the matrix sparse
 810     // is this a good idea?
 811     if(epsilonEqual(totalSurface, 0.0, SPARSE_TRUNCATION_LIMIT))
 812       {
 813         totalSurface = 0.0;
 814       }
 815
 816     LOG(2, "Volume = " << totalSurface << ", det= " << _t->determinant());
 817
 818     return totalSurface;
 819   }
 820
 821   /**
 822    * Calculates the volume of intersection of this tetrahedron with another one.
 823    */
 824   template<class MyMeshType>
 825   double SplitterTetra<MyMeshType>::intersectTetra(const double** tetraCorners)
 826   {
 827     //{ could be done on outside?
 828     // check if we have planar tetra element
 829     if(_t->determinant() == 0.0)
 830     {
 831       // tetra is planar
 832       LOG(2, "Planar tetra -- volume 0");
 833       return 0.0;
 834     }
 835
 836     const unsigned nbOfNodes4Type=4;
 837     // halfspace filtering
 838     bool isOutside[8] = {true, true, true, true, true, true, true, true};
 839     bool isTargetOutside = false;
 840
 841     // calculate the transformed coordinates of the nodes
 842     double nodes[nbOfNodes4Type][3];
 843     for(int i = 0;i<(int)nbOfNodes4Type;++i)
 844     {
 845       _t->apply(nodes[i], tetraCorners[i]);
 846       CheckIsOutside(nodes[i], isOutside);
 847     }
 848
 849     // halfspace filtering check
 850     // NB : might not be beneficial for caching of triangles
 851     for(int i = 0; i < 8; ++i)
 852     {
 853       if(isOutside[i])
 854       {
 855         isTargetOutside = true;
 856       }
 857     }
 858
 859     double totalVolume = 0.0;
 860
 861     if(!isTargetOutside)
 862     {
 863       const CellModel& cellModelCell=CellModel::GetCellModel(NORM_TETRA4);
 864       ConnType cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3];
 865
 866       for(unsigned ii = 0 ; ii < 4 ; ++ii)
 867       {
 868         cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
 869
 870         TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
 871         double vol = tri.calculateIntersectionVolume();
 872         LOG(1, "ii = " << ii << " Volume=" << vol)
 873         totalVolume += vol;
 874       }
 875
 876       // reset if it is very small to keep the matrix sparse
 877       // is this a good idea?
 878       if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
 879       {
 880         totalVolume = 0.0;
 881       }
 882     }
 883     LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
 884
 885     // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
 886     // that should be used (which is equivalent to dividing by the determinant)
 887     return std::fabs(1.0 / _t->determinant() * totalVolume) ;
 888   }
 889
 890   ////////////////////////////////////////////////////////
 891
 892   template<class MyMeshTypeT, class MyMeshTypeS>
 893   SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy)
 894     :_target_mesh(targetMesh),_src_mesh(srcMesh),_splitting_pol(policy)
 895   {
 896   }
 897
 898   template<class MyMeshTypeT, class MyMeshTypeS>
 899   SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::~SplitterTetra2()
 900   {
 901     releaseArrays();
 902   }
 903
 904   template<class MyMeshTypeT, class MyMeshTypeS>
 905   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::releaseArrays()
 906   {
 907     // free potential sub-mesh nodes that have been allocated
 908     if(_nodes.size()>=/*8*/_node_ids.size())
 909       {
 910         typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());
 911         std::vector<const double*>::iterator iter = _nodes.begin() + /*8*/nbOfNodesT;
 912         while(iter != _nodes.end())
 913           {
 914             delete[] *iter;
 915             ++iter;
 916           }
 917       }
 918     _nodes.clear();
 919   }
 920
 921   /*!
 922    * \param [in] targetCell in C mode.
 923    * \param [out] tetra is the output result tetra containers.
 924    */
 925   template<class MyMeshTypeT, class MyMeshTypeS>
 926   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
 927   {
 928     typedef typename MyMeshTypeT::MyConnType TConnType;
 929     const TConnType *refConn(_target_mesh.getConnectivityPtr());
 930     const TConnType *cellConn(refConn+_target_mesh.getConnectivityIndexPtr()[targetCell]);
 931     INTERP_KERNEL::NormalizedCellType gt(_target_mesh.getTypeOfElement(targetCell));
 932     std::vector<TConnType> tetrasNodalConn;
 933     std::vector<double> addCoords;
 934     const double *coords(_target_mesh.getCoordinatesPtr());
 935     SplitIntoTetras(_splitting_pol,gt,cellConn,refConn+_target_mesh.getConnectivityIndexPtr()[targetCell+1],coords,tetrasNodalConn,addCoords);
 936     std::size_t nbTetras(tetrasNodalConn.size()/4); tetra.resize(nbTetras);
 937     double tmp[12];
 938     typename MyMeshTypeS::MyConnType tmp2[4];
 939     for(std::size_t i=0;i<nbTetras;i++)
 940       {
 941         for(int j=0;j<4;j++)
 942           {
 943             TConnType cellId(tetrasNodalConn[4*i+j]);
 944             tmp2[j]=cellId;
 945             if(cellId>=0)
 946               {
 947                 tmp[j*3+0]=coords[3*cellId+0];
 948                 tmp[j*3+1]=coords[3*cellId+1];
 949                 tmp[j*3+2]=coords[3*cellId+2];
 950               }
 951             else
 952               {
 953                 tmp[j*3+0]=addCoords[3*(-cellId-1)+0];
 954                 tmp[j*3+1]=addCoords[3*(-cellId-1)+1];
 955                 tmp[j*3+2]=addCoords[3*(-cellId-1)+2];
 956               }
 957           }
 958         tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp,tmp2);
 959       }
 960   }
 961
 962   /*!
 963    * @param targetCell in C mode.
 964    * @param tetra is the output result tetra containers.
 965    */
 966   template<class MyMeshTypeT, class MyMeshTypeS>
 967   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell(typename MyMeshTypeT::MyConnType targetCell,
 968                                                                  typename MyMeshTypeT::MyConnType nbOfNodesT,
 969                                                                  typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
 970   {
 971     typedef typename MyMeshTypeT::MyConnType ConnType;
 972     const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
 973     const int numTetra = static_cast<int>(_splitting_pol);
 974     if(nbOfNodesT==4)
 975       {
 976         _nodes.resize(8);
 977         _node_ids.resize(8);
 978         tetra.reserve(1);
 979         const double *nodes[4];
 980         ConnType conn[4];
 981         for(int node = 0; node < 4 ; ++node)
 982           {
 983             nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
 984           }
 985         std::copy(conn,conn+4,_node_ids.begin());
 986         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
 987         tetra.push_back(t);
 988         return ;
 989       }
 990     // Issue 0020634. To pass nbOfNodesT to calculateSubNodes (don't want to add an arg)
 991     _node_ids.resize(nbOfNodesT);
 992
 993     // pre-calculate nodes
 994     calculateSubNodes(_target_mesh, OTT<ConnType,numPol>::indFC(targetCell));
 995
 996     tetra.reserve(numTetra);
 997     _nodes.reserve(30); // we never have more than this
 998
 999     switch ( nbOfNodesT )
1000       {
1001       case 8:
1002         {
1003           switch(_splitting_pol)
1004             {
1005             case PLANAR_FACE_5:
1006               {
1007                 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
1008                 fiveSplit(subZone,tetra);
1009               }
1010               break;
1011
1012             case PLANAR_FACE_6:
1013               {
1014                 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
1015                 sixSplit(subZone,tetra);
1016               }
1017               break;
1018
1019             case GENERAL_24:
1020               {
1021                 calculateGeneral24TetraOld(tetra);
1022               }
1023               break;
1024
1025             case GENERAL_48:
1026               {
1027                 calculateGeneral48Tetra(tetra);
1028               }
1029               break;
1030             default:
1031               assert(false);
1032             }
1033           break;
1034         }
1035       case 5:
1036         {
1037           splitPyram5(tetra);
1038           break;
1039         }
1040       default:
1041         {
1042           splitConvex(targetCell, tetra);
1043         }
1044       }
1045   }
1046
1047   /**
1048    * Splits the hexahedron into five tetrahedra.
1049    * This method adds five SplitterTetra objects to the vector tetra.
1050    *
1051    * @param subZone  the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1052    *                 splitting to be reused on the subzones of the GENERAL_* types of splitting
1053    */
1054   template<class MyMeshTypeT, class MyMeshTypeS>
1055   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1056   {
1057     // create tetrahedra
1058     for(int i = 0; i < 5; ++i)
1059       {
1060         const double* nodes[4];
1061         typename MyMeshTypeS::MyConnType conn[4];
1062         for(int j = 0; j < 4; ++j)
1063           {
1064             conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ];
1065             typename MyMeshTypeS::MyConnType realConn;
1066             nodes[j] = getCoordsOfSubNode2(conn[j],realConn);
1067             conn[j] = realConn;
1068           }
1069         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1070         tetra.push_back(t);
1071       }
1072   }
1073
1074   template<class MyMeshTypeT, class MyMeshTypeS>
1075   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1076   {
1077     sixSplitGen(subZone,tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType& conn, const double *&coords)
1078     {
1079       typename MyMeshTypeS::MyConnType realConn;
1080       coords = obj.getCoordsOfSubNode2(conn,realConn);
1081       conn = realConn;
1082     });
1083   }
1084
1085   /**
1086    * Splits the hexahedron into six tetrahedra.
1087    * This method adds six SplitterTetra objects to the vector tetra.
1088    *
1089    * @param subZone  the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1090    *                 splitting to be reused on the subzones of the GENERAL_* types of splitting
1091    */
1092   template<class MyMeshTypeT, class MyMeshTypeS>
1093   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplitGen(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType&, const double*&)> func)
1094   {
1095     for(int i = 0; i < 6; ++i)
1096       {
1097         const double* nodes[4];
1098         typename MyMeshTypeS::MyConnType conn[4];
1099         for(int j = 0; j < 4; ++j)
1100           {
1101             conn[j] = subZone[SPLIT_NODES_6[4*i+j]];
1102             func(*this,conn[j],nodes[j]);
1103           }
1104         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1105         tetra.push_back(t);
1106       }
1107   }
1108
1109   /**
1110    * Version of calculateGeneral24Tetra connectivity aware for P1P0 and P0P1
1111    */
1112   template<class MyMeshTypeT, class MyMeshTypeS>
1113   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1114   {
1115     calculateGeneral24TetraGen(tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType conn[4], const double* nodes[4]) {
1116       typename MyMeshTypeS::MyConnType realConn;
1117       nodes[2] = obj.getCoordsOfSubNode2(conn[2],realConn); conn[2] = realConn;
1118       nodes[3] = obj.getCoordsOfSubNode2(conn[3],realConn); conn[3] = realConn;
1119     });
1120   }
1121
1122   /*!
1123    * Version for 3D2DP0P0
1124    */
1125   template<class MyMeshTypeT, class MyMeshTypeS>
1126   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24TetraOld(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1127   {
1128     calculateGeneral24TetraGen(tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType conn[4], const double* nodes[4]) {
1129       nodes[2] = obj.getCoordsOfSubNode(conn[2]);
1130       nodes[3] = obj.getCoordsOfSubNode(conn[3]);
1131     });
1132   }
1133
1134   /**
1135    * Splits the hexahedron into 24 tetrahedra.
1136    * The splitting is done by combining the barycenter of the tetrahedron, the barycenter of each face
1137    * and the nodes of each edge of the face. This creates 6 faces * 4 edges / face = 24 tetrahedra.
1138    * The submesh nodes introduced are the barycenters of the faces and the barycenter of the cell.
1139    */
1140   template<class MyMeshTypeT, class MyMeshTypeS>
1141   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24TetraGen(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType[4], const double*[4])> func)
1142   {
1143     // The two nodes of the original mesh cell used in each tetrahedron.
1144     // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
1145     // For the correspondence of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
1146
1147     // nodes to use for tetrahedron
1148     const double* nodes[4];
1149     typename MyMeshTypeS::MyConnType conn[4];
1150     // get the cell center
1151     conn[0] = 14;
1152     nodes[0] = getCoordsOfSubNode(conn[0]);
1153
1154     for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode)
1155       {
1156         // get the face center
1157         conn[1] = faceCenterNode;
1158         nodes[1] = getCoordsOfSubNode(conn[1]);
1159         for(int j = 0; j < 4; ++j)
1160           {
1161             const int row = 4*(faceCenterNode - 8) + j;
1162             conn[2] = TETRA_EDGES_GENERAL_24[2*row];
1163             conn[3] = TETRA_EDGES_GENERAL_24[2*row + 1];
1164             func(*this,conn,nodes);
1165             SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes, conn);
1166             tetra.push_back(t);
1167           }
1168       }
1169   }
1170
1171   /**
1172    * Splits the hexahedron into 48 tetrahedra.
1173    * The splitting is done by introducing the midpoints of all the edges
1174    * and the barycenter of the element as submesh nodes. The 8 hexahedral subzones thus defined
1175    * are then split into 6 tetrahedra each, as in Grandy, p. 449. The division of the subzones
1176    * is done by calling sixSplit().
1177    *
1178    */
1179   template<class MyMeshTypeT, class MyMeshTypeS>
1180   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1181   {
1182     for(int i = 0; i < 8; ++i)
1183       {
1184         sixSplitGen(GENERAL_48_SUBZONES+8*i,tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType& conn, const double *&coords){
1185           coords = obj.getCoordsOfSubNode(conn);
1186         });
1187       }
1188   }
1189
1190   /**
1191    * Splits the NORM_PYRA5 into 2 tetrahedra.
1192    */
1193   template<class MyMeshTypeT, class MyMeshTypeS>
1194   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1195   {
1196     static const int SPLIT_PYPA5[2][4] =
1197       {
1198         {
1199           0, 1, 2, 4
1200         },
1201         {
1202           0, 2, 3, 4
1203         }
1204       };
1205
1206     // create tetrahedra
1207     const double* nodes[4];
1208     typename MyMeshTypeS::MyConnType conn[4];
1209     for(int i = 0; i < 2; ++i)
1210       {
1211         for(int j = 0; j < 4; ++j)
1212           nodes[j] = getCoordsOfSubNode2(SPLIT_PYPA5[i][j],conn[j]);
1213         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1214         tetra.push_back(t);
1215       }
1216   }
1217
1218   /**
1219    * Splits a convex cell into tetrahedra.
1220    */
1221   template<class MyMeshTypeT, class MyMeshTypeS>
1222   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitConvex(typename MyMeshTypeT::MyConnType targetCell,
1223                                                              typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1224   {
1225     // Each face of a cell is split into triangles and
1226     // each of triangles and a cell barycenter form a tetrahedron.
1227
1228     typedef typename MyMeshTypeT::MyConnType ConnType;
1229     const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
1230
1231     // get type of cell and nb of cell nodes
1232     NormalizedCellType normCellType=_target_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(targetCell));
1233     const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
1234     ConnType nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
1235
1236     // get nb of cell sons (faces)
1237     const ConnType* rawCellConn = _target_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _target_mesh.getConnectivityIndexPtr()[ targetCell ]);
1238     const ConnType rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
1239     unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
1240
1241     // indices of nodes of a son
1242     static std::vector<ConnType> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
1243     while ( allNodeIndices.size() < (std::size_t)nbOfCellNodes )
1244       allNodeIndices.push_back( static_cast<ConnType>(allNodeIndices.size()) );
1245     std::vector<ConnType> classicFaceNodes(4);
1246     if(cellModelCell.isQuadratic())
1247       throw INTERP_KERNEL::Exception("SplitterTetra2::splitConvex : quadratic 3D cells are not implemented yet !");
1248     ConnType* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
1249
1250     // nodes of tetrahedron
1251     typename MyMeshTypeS::MyConnType conn[4];
1252     const double* nodes[4];
1253     nodes[3] = getCoordsOfSubNode2( nbOfCellNodes,conn[3]); // barycenter
1254
1255     for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
1256       {
1257         // get indices of son's nodes: it's just next portion of allNodeIndices for polyhedron
1258         // and some of allNodeIndices accodring to cell model for a classsic cell
1259         unsigned nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon2(ii, rawCellConn, rawNbCellNodes);
1260         if ( normCellType != NORM_POLYHED )
1261           cellModelCell.fillSonCellNodalConnectivity(ii,&allNodeIndices[0],faceNodes);
1262
1263         int nbTetra = nbFaceNodes - 2; // split polygon into nbTetra triangles
1264
1265         // create tetrahedra
1266         for(int i = 0; i < nbTetra; ++i)
1267           {
1268             nodes[0] = getCoordsOfSubNode2( faceNodes[0],  conn[0]);
1269             nodes[1] = getCoordsOfSubNode2( faceNodes[1+i],conn[1]);
1270             nodes[2] = getCoordsOfSubNode2( faceNodes[2+i],conn[2]);
1271             SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1272             tetra.push_back(t);
1273           }
1274
1275         if ( normCellType == NORM_POLYHED )
1276           faceNodes += nbFaceNodes; // go to the next face
1277       }
1278   }
1279
1280   /**
1281    * Precalculates all the nodes.
1282    * Retrieves the mesh nodes and allocates the necessary sub-mesh
1283    * nodes according to the splitting policy used.
1284    * This method is meant to be called once by the constructor.
1285    *
1286    * @param targetMesh  the target mesh
1287    * @param targetCell  the global number of the cell that the object represents, in targetMesh mode.
1288    *
1289    */
1290   template<class MyMeshTypeT, class MyMeshTypeS>
1291   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell)
1292   {
1293     // retrieve real mesh nodes
1294
1295     typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());// Issue 0020634. _node_ids.resize(8);
1296     for(int node = 0; node < nbOfNodesT ; ++node)
1297       {
1298         // calculate only normal nodes
1299         _nodes.push_back(getCoordsOfNode2(node, targetCell, targetMesh,_node_ids[node]));
1300       }
1301
1302     switch ( nbOfNodesT )
1303       {
1304       case 8:
1305
1306         // create sub-mesh nodes if needed
1307         switch(_splitting_pol)
1308           {
1309           case GENERAL_24:
1310             {
1311               for(int i = 0; i < 7; ++i)
1312                 {
1313                   double* barycenter = new double[3];
1314                   calcBarycenter(4, barycenter, &GENERAL_24_SUB_NODES[4*i]);
1315                   _nodes.push_back(barycenter);
1316                 }
1317             }
1318             break;
1319
1320           case GENERAL_48:
1321             {
1322               for(int i = 0; i < 19; ++i)
1323                 {
1324                   double* barycenter = new double[3];
1325                   calcBarycenter(2, barycenter, &GENERAL_48_SUB_NODES[2*i]);
1326                   _nodes.push_back(barycenter);
1327                 }
1328             }
1329             break;
1330
1331           default:
1332             break;
1333           }
1334         break;
1335
1336       case 5: // NORM_PYRA5
1337         break;
1338
1339       default: // convex 3d cell
1340         {
1341           // add barycenter of a cell
1342           std::vector<int> allIndices(nbOfNodesT);
1343           for ( int i = 0; i < nbOfNodesT; ++i ) allIndices[i] = i;
1344           double* barycenter = new double[3];
1345           calcBarycenter(nbOfNodesT, barycenter, &allIndices[0]);
1346           _nodes.push_back(barycenter);
1347         }
1348       }
1349
1350   }
1351 }
1352
1353 #endif