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