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