src/INTERP_KERNEL/SplitterTetra.txx

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