src/INTERP_KERNEL/SplitterTetra.txx

   1 // Copyright (C) 2007-2019  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #ifndef __SPLITTERTETRA_TXX__
  20 #define __SPLITTERTETRA_TXX__
  21
  22 #include "SplitterTetra.hxx"
  23
  24 #include "TetraAffineTransform.hxx"
  25 #include "TransformedTriangle.hxx"
  26 #include "MeshUtils.hxx"
  27 #include "VectorUtils.hxx"
  28 #include "CellModel.hxx"
  29 #include "Log.hxx"
  30 #include "UnitTetraIntersectionBary.hxx"
  31 #include "VolSurfFormulae.hxx"
  32
  33 #include <cmath>
  34 #include <cassert>
  35 #include <string>
  36 #include <sstream>
  37 #include <vector>
  38
  39 namespace INTERP_KERNEL
  40 {
  41   template<class MyMeshType>
  42   const double SplitterTetra<MyMeshType>::SPARSE_TRUNCATION_LIMIT=1.0e-14;
  43
  44   /*!
  45    * output is expected to be allocated with 24*sizeof(void*) in order to store the 24 tetras.
  46    * These tetras have to be deallocated.
  47    */
  48   template<class MyMeshType>
  49   void SplitterTetra<MyMeshType>::splitIntoDualCells(SplitterTetra<MyMeshType> **output)
  50   {
  51     double tmp[12];
  52     const double *tmp2[4]={tmp,tmp+3,tmp+6,tmp+9};
  53     typename MyMeshType::MyConnType conn[4]={-1,-1,-1,-1};
  54     for(int i=0;i<24;i++)
  55       {
  56         splitMySelfForDual(tmp,i,conn[0]);
  57         output[i]=new SplitterTetra<MyMeshType>(_src_mesh,tmp2,conn);
  58       }
  59   }
  60
  61   /**
  62    * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a
  63    * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
  64    * if a partial computation of a triangle has already been performed (to increase performance).
  65    *
  66    * The \a srcMesh can contain polyhedron cells.
  67    *
  68    *
  69    * Constructor creating object from the four corners of the tetrahedron.
  70    *
  71    * @param srcMesh       mesh containing the source elements
  72    * @param tetraCorners  array of four pointers to double[3] arrays containing the coordinates of the
  73    *                      corners of the tetrahedron
  74    */
  75   template<class MyMeshType>
  76   SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId)
  77     : _t(0), _src_mesh(srcMesh)
  78   {
  79     std::copy(nodesId,nodesId+4,_conn);
  80     _coords[0]=tetraCorners[0][0]; _coords[1]=tetraCorners[0][1]; _coords[2]=tetraCorners[0][2];
  81     _coords[3]=tetraCorners[1][0]; _coords[4]=tetraCorners[1][1]; _coords[5]=tetraCorners[1][2];
  82     _coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2];
  83     _coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2];
  84     // create the affine transform
  85     _t=new TetraAffineTransform(_coords);
  86   }
  87
  88   /**
  89    * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a
  90    * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
  91    * if a partial computation of a triangle has already been performed (to increase performance).
  92    *
  93    * The \a srcMesh can contain polyhedron cells.
  94    *
  95    *
  96    * Constructor creating object from the four corners of the tetrahedron.
  97    *
  98    * \param [in] srcMesh       mesh containing the source elements
  99    * \param [in] tetraCorners  array 4*3 doubles containing corners of input tetrahedron (P0X,P0Y,P0Y,P1X,P1Y,P1Z,P2X,P2Y,P2Z,P3X,P3Y,P3Z).
 100    */
 101   template<class MyMeshType>
 102   SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const int *conn): _t(0),_src_mesh(srcMesh)
 103   {
 104     if(!conn)
 105       { _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3; }
 106     else
 107       { _conn[0]=conn[0]; _conn[1]=conn[1]; _conn[2]=conn[2]; _conn[3]=conn[3]; }
 108     _coords[0]=tetraCorners[0]; _coords[1]=tetraCorners[1]; _coords[2]=tetraCorners[2]; _coords[3]=tetraCorners[3]; _coords[4]=tetraCorners[4]; _coords[5]=tetraCorners[5];
 109     _coords[6]=tetraCorners[6]; _coords[7]=tetraCorners[7]; _coords[8]=tetraCorners[8]; _coords[9]=tetraCorners[9]; _coords[10]=tetraCorners[10]; _coords[11]=tetraCorners[11];
 110     // create the affine transform
 111     _t=new TetraAffineTransform(_coords);
 112   }
 113
 114   /**
 115    * Destructor
 116    *
 117    * Deletes _t and the coordinates (double[3] vectors) in _nodes
 118    *
 119    */
 120   template<class MyMeshType>
 121   SplitterTetra<MyMeshType>::~SplitterTetra()
 122   {
 123     delete _t;
 124     for(HashMap< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
 125       delete[] iter->second;
 126   }
 127
 128   /*!
 129    * \Forget already calculated triangles, which is crucial for calculation of barycenter of intersection
 130    */
 131   template<class MyMeshType>
 132   void SplitterTetra<MyMeshType>::clearVolumesCache()
 133   {
 134     _volumes.clear();
 135   }
 136
 137   /*!
 138    * This method destroys the 4 pointers pointed by tetraCorners[0],tetraCorners[1],tetraCorners[2] and tetraCorners[3]
 139    * @param i is in 0..23 included.
 140    * @param output is expected to be sized of 12 in order to.
 141    */
 142   template<class MyMeshType>
 143   void SplitterTetra<MyMeshType>::splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId)
 144   {
 145     double *tmp[4];
 146     int offset=i/6;
 147     nodeId=_conn[offset];
 148     tmp[0]=_coords+3*offset; tmp[1]=_coords+((offset+1)%4)*3; tmp[2]=_coords+((offset+2)%4)*3; tmp[3]=_coords+((offset+3)%4)*3;
 149     int caseToTreat=i%6;
 150     int case1=caseToTreat/2;
 151     int case2=caseToTreat%2;
 152     const int tab[3][2]={{1,2},{3,2},{1,3}};
 153     const int *curTab=tab[case1];
 154     double pt0[3]; pt0[0]=(tmp[curTab[case2]][0]+tmp[0][0])/2.; pt0[1]=(tmp[curTab[case2]][1]+tmp[0][1])/2.; pt0[2]=(tmp[curTab[case2]][2]+tmp[0][2])/2.;
 155     double pt1[3]; pt1[0]=(tmp[0][0]+tmp[curTab[0]][0]+tmp[curTab[1]][0])/3.; pt1[1]=(tmp[0][1]+tmp[curTab[0]][1]+tmp[curTab[1]][1])/3.; pt1[2]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.;
 156     double pt2[3]; pt2[0]=(tmp[0][0]+tmp[1][0]+tmp[2][0]+tmp[3][0])/4.; pt2[1]=(tmp[0][1]+tmp[1][1]+tmp[2][1]+tmp[3][1])/4.; pt2[2]=(tmp[0][2]+tmp[1][2]+tmp[2][2]+tmp[3][2])/4.;
 157     std::copy(pt1,pt1+3,output+case2*3);
 158     std::copy(pt0,pt0+3,output+(abs(case2-1))*3);
 159     std::copy(pt2,pt2+3,output+2*3);
 160     std::copy(tmp[0],tmp[0]+3,output+3*3);
 161   }
 162
 163   /**
 164    * Calculates the volume of intersection of an element in the source mesh and the target element.
 165    * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
 166    * faces of the source element are triangulated and the calculated transformation is applied
 167    * to each triangle. The algorithm of Grandy, implemented in INTERP_KERNEL::TransformedTriangle is used
 168    * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
 169    * divided by the determinant of the transformation.
 170    *
 171    * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated
 172    * with triangulated faces to avoid having to recalculate these.
 173    *
 174    * @param element      global number of the source element in C mode.
 175    */
 176   template<class MyMeshType>
 177   double SplitterTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element,
 178                                                         double*                         baryCentre)
 179   {
 180     typedef typename MyMeshType::MyConnType ConnType;
 181     const NumberingPolicy numPol=MyMeshType::My_numPol;
 182     //{ could be done on outside?
 183     // check if we have planar tetra element
 184     if(_t->determinant() == 0.0)
 185       {
 186         // tetra is planar
 187         LOG(2, "Planar tetra -- volume 0");
 188         return 0.0;
 189       }
 190
 191     // get type of cell
 192     NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
 193     const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
 194     unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
 195     // halfspace filtering
 196     bool isOutside[8] = {true, true, true, true, true, true, true, true};
 197     bool isTargetOutside = false;
 198
 199     // calculate the coordinates of the nodes
 200     int *cellNodes=new int[nbOfNodes4Type];
 201     for(int i = 0;i<(int)nbOfNodes4Type;++i)
 202       {
 203         // we could store mapping local -> global numbers too, but not sure it is worth it
 204         const int globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
 205         cellNodes[i]=globalNodeNum;
 206         if(_nodes.find(globalNodeNum) == _nodes.end())
 207           {
 208             //for(HashMap< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
 209             //  std::cout << (*iter3).first << " ";
 210             //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
 211             calculateNode(globalNodeNum);
 212           }
 213         CheckIsOutside(_nodes[globalNodeNum], isOutside);
 214       }
 215
 216     // halfspace filtering check
 217     // NB : might not be beneficial for caching of triangles
 218     for(int i = 0; i < 8; ++i)
 219       {
 220         if(isOutside[i])
 221           {
 222             isTargetOutside = true;
 223           }
 224       }
 225
 226     double totalVolume = 0.0;
 227
 228     if(!isTargetOutside)
 229       {
 230         /// calculator of intersection barycentre
 231         UnitTetraIntersectionBary baryCalculator( _t->determinant() < 0.);
 232
 233         // get nb of sons of a cell
 234         const ConnType* rawCellConn = _src_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _src_mesh.getConnectivityIndexPtr()[ element ]);
 235         const int rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
 236         unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
 237
 238         for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
 239           {
 240             // get sons connectivity
 241             NormalizedCellType faceType;
 242             int *faceNodes, nbFaceNodes=-1;
 243             if ( cellModelCell.isDynamic() )
 244               {
 245                 faceNodes=new int[nbOfNodes4Type];
 246                 nbFaceNodes = cellModelCell.fillSonCellNodalConnectivity2(ii,rawCellConn,rawNbCellNodes,faceNodes,faceType);
 247                 for ( int i = 0; i < nbFaceNodes; ++i )
 248                   faceNodes[i] = OTT<ConnType,numPol>::coo2C(faceNodes[i]);
 249               }
 250             else
 251               {
 252                 faceType = cellModelCell.getSonType(ii);
 253                 const CellModel& faceModel=CellModel::GetCellModel(faceType);
 254                 assert(faceModel.getDimension() == 2);
 255                 nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon(ii);
 256                 faceNodes = new int[nbFaceNodes];
 257                 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
 258               }
 259             // intersect a son with the unit tetra
 260             switch(faceType)
 261               {
 262               case NORM_TRI3:
 263                 {
 264                   // create the face key
 265                   TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
 266
 267                   // calculate the triangle if needed
 268                   if(_volumes.find(key) == _volumes.end())
 269                     {
 270                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
 271                       calculateVolume(tri, key);
 272                       totalVolume += _volumes[key];
 273                       if ( baryCentre )
 274                         baryCalculator.addSide( tri );
 275                     } else {
 276                       // count negative as face has reversed orientation
 277                       totalVolume -= _volumes[key];
 278                     }
 279                 }
 280                 break;
 281
 282               case NORM_QUAD4:
 283
 284                 // simple triangulation of faces along a diagonal :
 285                 //
 286                 // 2 ------ 3
 287                 // |      / |
 288                 // |     /  |
 289                 // |    /   |
 290                 // |   /    |
 291                 // |  /     |
 292                 // | /      |
 293                 // 1 ------ 4
 294                 //
 295                 //? not sure if this always works
 296                 {
 297                   // calculate the triangles if needed
 298
 299                   // local nodes 1, 2, 3
 300                   TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
 301                   if(_volumes.find(key1) == _volumes.end())
 302                     {
 303                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
 304                       calculateVolume(tri, key1);
 305                       totalVolume += _volumes[key1];
 306                     } else {
 307                       // count negative as face has reversed orientation
 308                       totalVolume -= _volumes[key1];
 309                     }
 310
 311                   // local nodes 1, 3, 4
 312                   TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]);
 313                   if(_volumes.find(key2) == _volumes.end())
 314                     {
 315                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]);
 316                       calculateVolume(tri, key2);
 317                       totalVolume += _volumes[key2];
 318                     }
 319                   else
 320                     {
 321                       // count negative as face has reversed orientation
 322                       totalVolume -= _volumes[key2];
 323                     }
 324                 }
 325                 break;
 326
 327               case NORM_POLYGON:
 328                 {
 329                   int nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
 330                   for ( int iTri = 0; iTri < nbTria; ++iTri )
 331                     {
 332                       TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]);
 333                       if(_volumes.find(key) == _volumes.end())
 334                         {
 335                           TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1+iTri]], _nodes[faceNodes[2+iTri]]);
 336                           calculateVolume(tri, key);
 337                           totalVolume += _volumes[key];
 338                         }
 339                       else
 340                         {
 341                           totalVolume -= _volumes[key];
 342                         }
 343                     }
 344                 }
 345                 break;
 346
 347               default:
 348                 std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
 349                 assert(false);
 350               }
 351             delete [] faceNodes;
 352           }
 353
 354         if ( baryCentre ) {
 355           baryCalculator.getBary( baryCentre );
 356           _t->reverseApply( baryCentre, baryCentre );
 357         }
 358       }
 359     delete [] cellNodes;
 360     // reset if it is very small to keep the matrix sparse
 361     // is this a good idea?
 362     if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
 363       {
 364         totalVolume = 0.0;
 365       }
 366
 367     LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
 368
 369     // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
 370     // that should be used (which is equivalent to dividing by the determinant)
 371     return std::fabs(1.0 / _t->determinant() * totalVolume) ;
 372   }
 373
 374   /**
 375    * Calculates the intersection surface of two coplanar triangles.
 376    *
 377    * @param palneNormal normal of the plane for the first triangle
 378    * @param planeConstant constant of the equation of the plane for the first triangle
 379    * @param p1 coordinates of the first  node of the first  triangle
 380    * @param p2 coordinates of the second node of the first  triangle
 381    * @param p3 coordinates of the third  node of the first  triangle
 382    * @param p4 coordinates of the first  node of the second triangle
 383    * @param p5 coordinates of the second node of the second triangle
 384    * @param p6 coordinates of the third  node of the second triangle
 385    * @param dimCaracteristic characteristic size of the meshes containing the triangles
 386    * @param precision precision for double float data used for comparison
 387    */
 388   template<class MyMeshType>
 389   double SplitterTetra<MyMeshType>::CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
 390                                                                                     const double planeConstant,
 391                                                                                     const double *const p1, const double *const p2, const double *const p3,
 392                                                                                     const double *const p4, const double *const p5, const double *const p6,
 393                                                                                     const double dimCaracteristic, const double precision)
 394   {
 395     typedef typename MyMeshType::MyConnType ConnType;
 396     typedef double Vect2[2];
 397     typedef double Triangle2[3][2];
 398
 399     const double *const tri0[3] = {p1, p2, p3};
 400     const double *const tri1[3] = {p4, p5, p6};
 401
 402     // Plane of the first triangle defined by the normal of the triangle and the constant
 403     // Project triangles onto coordinate plane most aligned with plane normal
 404     int maxNormal = 0;
 405     double fmax = std::abs(planeNormal[0]);
 406     double absMax = std::abs(planeNormal[1]);
 407     if (absMax > fmax)
 408       {
 409         maxNormal = 1;
 410         fmax = absMax;
 411       }
 412     absMax = std::abs(planeNormal[2]);
 413     if (absMax > fmax)
 414       {
 415         maxNormal = 2;
 416       }
 417
 418     Triangle2 projTri0, projTri1;
 419     int i;
 420
 421     if (maxNormal == 0)
 422       {
 423         // Project onto yz-plane.
 424         for (i = 0; i < 3; ++i)
 425           {
 426             projTri0[i][0] = tri0[i][1];
 427             projTri0[i][1] = tri0[i][2];
 428             projTri1[i][0] = tri1[i][1];
 429             projTri1[i][1] = tri1[i][2];
 430           }
 431       }
 432     else if (maxNormal == 1)
 433       {
 434         // Project onto xz-plane.
 435         for (i = 0; i < 3; ++i)
 436           {
 437             projTri0[i][0] = tri0[i][0];
 438             projTri0[i][1] = tri0[i][2];
 439             projTri1[i][0] = tri1[i][0];
 440             projTri1[i][1] = tri1[i][2];
 441           }
 442       }
 443     else
 444       {
 445         // Project onto xy-plane.
 446         for (i = 0; i < 3; ++i)
 447           {
 448             projTri0[i][0] = tri0[i][0];
 449             projTri0[i][1] = tri0[i][1];
 450             projTri1[i][0] = tri1[i][0];
 451             projTri1[i][1] = tri1[i][1];
 452           }
 453       }
 454
 455     // 2D triangle intersection routines require counterclockwise ordering.
 456     Vect2 save;
 457     Vect2 edge0;
 458     Vect2 edge1;
 459     for (int ii = 0; ii < 2; ++ii)
 460       {
 461         edge0[ii] = projTri0[1][ii] - projTri0[0][ii];
 462         edge1[ii] = projTri0[2][ii] - projTri0[0][ii];
 463       }
 464     if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
 465       {
 466         // Triangle is clockwise, reorder it.
 467         for (int ii = 0; ii < 2; ++ii)
 468           {
 469             save[ii] = projTri0[1][ii];
 470             projTri0[1][ii] = projTri0[2][ii];
 471             projTri0[2][ii] = save[ii];
 472           }
 473       }
 474
 475     for (int ii = 0; ii < 2; ++ii)
 476       {
 477         edge0[ii] = projTri1[1][ii] - projTri1[0][ii];
 478         edge1[ii] = projTri1[2][ii] - projTri1[0][ii];
 479       }
 480     if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
 481       {
 482         // Triangle is clockwise, reorder it.
 483       for (int ii = 0; ii < 2; ++ii)
 484         {
 485           save[ii] = projTri1[1][ii];
 486           projTri1[1][ii] = projTri1[2][ii];
 487           projTri1[2][ii] = save[ii];
 488         }
 489       }
 490
 491     std::vector<double> inter2;
 492     intersec_de_triangle(projTri0[0], projTri0[1], projTri0[2],
 493                          projTri1[0], projTri1[1], projTri1[2],
 494                          inter2,
 495                          dimCaracteristic, precision);
 496     ConnType nb_inter=((ConnType)inter2.size())/2;
 497     double surface = 0.;
 498     if(nb_inter >3) inter2=reconstruct_polygon(inter2);
 499     if (nb_inter > 0)
 500       {
 501         std::vector<double> inter3;
 502         inter3.resize(3 * nb_inter);
 503         // Map 2D intersections back to the 3D triangle space.
 504         if (maxNormal == 0)
 505           {
 506             double invNX = ((double) 1.) / planeNormal[0];
 507             for (i = 0; i < nb_inter; i++)
 508               {
 509                 inter3[3 * i + 1] = inter2[2 * i];
 510                 inter3[3 * i + 2] = inter2[2 * i + 1];
 511                 inter3[3 * i] = invNX * (planeConstant - planeNormal[1] * inter3[3 * i + 1] - planeNormal[2] * inter3[3 * i + 2]);
 512               }
 513           }
 514         else if (maxNormal == 1)
 515           {
 516             double invNY = ((double) 1.) / planeNormal[1];
 517             for (i = 0; i < nb_inter; i++)
 518               {
 519                 inter3[3 * i] = inter2[2 * i];
 520                 inter3[3 * i + 2] = inter2[2 * i + 1];
 521                 inter3[3 * i + 1] = invNY * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[2] * inter3[3 * i + 2]);
 522               }
 523           }
 524         else
 525           {
 526             double invNZ = ((double) 1.) / planeNormal[2];
 527             for (i = 0; i < nb_inter; i++)
 528               {
 529                 inter3[3 * i] = inter2[2 * i];
 530                 inter3[3 * i + 1] = inter2[2 * i + 1];
 531                 inter3[3 * i + 2] = invNZ * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[1] * inter3[3 * i + 1]);
 532               }
 533           }
 534         surface = polygon_area<3>(inter3);
 535       }
 536     return surface;
 537   }
 538
 539   /**
 540    * Determine if a face is coplanar with a triangle.
 541    * The first face is characterized by the equation of her plane
 542    *
 543    * @param palneNormal normal of the plane for the first triangle
 544    * @param planeConstant constant of the equation of the plane for the first triangle
 545    * @param coordsFace coordinates of the triangle face
 546    * @param precision precision for double float data used for comparison
 547    */
 548   template<class MyMeshType>
 549   bool SplitterTetra<MyMeshType>::IsFacesCoplanar(const double *const planeNormal,
 550                                                   const double planeConstant,
 551                                                   const double *const *const coordsFace,
 552                                                   const double precision)
 553   {
 554       // Compute the signed distances of triangle vertices to the plane. Use an epsilon-thick plane test.
 555       // For faces not left
 556       int counter = 0;
 557       for (int i = 0; i < 3; ++i)
 558         {
 559           const double distance = dot(planeNormal, coordsFace[i]) - planeConstant;
 560           if (epsilonEqual(distance, precision))
 561             {
 562               counter++;
 563             }
 564         }
 565       return counter == 3;
 566   }
 567
 568   /**
 569    * Calculates the surface of intersection of a polygon face in the source mesh and a cell of the target mesh.
 570    * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
 571    * faces of the source element are triangulated and the calculated transformation is applied
 572    * to each triangle.
 573    * The algorithm is based on the algorithm of Grandy used in intersectSourceCell to compute
 574    * the volume of intersection of two cell elements.
 575    * The case with a source face colinear to one of the face of tetrahedrons is taking into account:
 576    * the contribution of the face must not be counted two times.
 577    *
 578    * The class will cache the intermediary calculations of transformed nodes of source faces and surfaces associated
 579    * with triangulated faces to avoid having to recalculate these.
 580    *
 581    * @param polyType type of the polygon source face
 582    * @param polyNodesNbr number of the nodes of the polygon source face
 583    * @param polyNodes numbers of the nodes of the polygon source face
 584    * @param polyCoords coordinates of the nodes of the polygon source face
 585    * @param dimCaracteristic characteristic size of the meshes containing the triangles
 586    * @param precision precision for double float data used for comparison
 587    * @param listOfTetraFacesTreated list of tetra faces treated
 588    * @param listOfTetraFacesColinear list of tetra faces colinear with the polygon source faces
 589    */
 590   template<class MyMeshType>
 591   double SplitterTetra<MyMeshType>::intersectSourceFace(const NormalizedCellType polyType,
 592                                                         const int polyNodesNbr,
 593                                                         const int *const polyNodes,
 594                                                         const double *const *const polyCoords,
 595                                                         const double dimCaracteristic,
 596                                                         const double precision,
 597                                                         std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
 598                                                         std::set<TriangleFaceKey>& listOfTetraFacesColinear)
 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     int tmp2[4];
 942     for(std::size_t i=0;i<nbTetras;i++)
 943       {
 944         for(int j=0;j<4;j++)
 945           {
 946             int cellId(tetrasNodalConn[4*i+j]);
 947             tmp2[j]=cellId;
 948             if(cellId>=0)
 949               {
 950                 tmp[j*3+0]=coords[3*cellId+0];
 951                 tmp[j*3+1]=coords[3*cellId+1];
 952                 tmp[j*3+2]=coords[3*cellId+2];
 953               }
 954             else
 955               {
 956                 tmp[j*3+0]=addCoords[3*(-cellId-1)+0];
 957                 tmp[j*3+1]=addCoords[3*(-cellId-1)+1];
 958                 tmp[j*3+2]=addCoords[3*(-cellId-1)+2];
 959               }
 960           }
 961         tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp,tmp2);
 962       }
 963   }
 964
 965   /*!
 966    * @param targetCell in C mode.
 967    * @param tetra is the output result tetra containers.
 968    */
 969   template<class MyMeshTypeT, class MyMeshTypeS>
 970   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell(typename MyMeshTypeT::MyConnType targetCell,
 971                                                                  typename MyMeshTypeT::MyConnType nbOfNodesT,
 972                                                                  typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
 973   {
 974     typedef typename MyMeshTypeT::MyConnType ConnType;
 975     const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
 976     const int numTetra = static_cast<int>(_splitting_pol);
 977     if(nbOfNodesT==4)
 978       {
 979         _nodes.resize(8);
 980         _node_ids.resize(8);
 981         tetra.reserve(1);
 982         const double *nodes[4];
 983         int conn[4];
 984         for(int node = 0; node < 4 ; ++node)
 985           {
 986             nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
 987           }
 988         std::copy(conn,conn+4,_node_ids.begin());
 989         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
 990         tetra.push_back(t);
 991         return ;
 992       }
 993     // Issue 0020634. To pass nbOfNodesT to calculateSubNodes (don't want to add an arg)
 994     _node_ids.resize(nbOfNodesT);
 995
 996     // pre-calculate nodes
 997     calculateSubNodes(_target_mesh, OTT<ConnType,numPol>::indFC(targetCell));
 998
 999     tetra.reserve(numTetra);
1000     _nodes.reserve(30); // we never have more than this
1001
1002     switch ( nbOfNodesT )
1003       {
1004       case 8:
1005         {
1006           switch(_splitting_pol)
1007             {
1008             case PLANAR_FACE_5:
1009               {
1010                 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
1011                 fiveSplit(subZone,tetra);
1012               }
1013               break;
1014
1015             case PLANAR_FACE_6:
1016               {
1017                 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
1018                 sixSplit(subZone,tetra);
1019               }
1020               break;
1021
1022             case GENERAL_24:
1023               {
1024                 calculateGeneral24Tetra(tetra);
1025               }
1026               break;
1027
1028             case GENERAL_48:
1029               {
1030                 calculateGeneral48Tetra(tetra);
1031               }
1032               break;
1033             default:
1034               assert(false);
1035             }
1036           break;
1037         }
1038       case 5:
1039         {
1040           splitPyram5(tetra);
1041           break;
1042         }
1043       default:
1044         {
1045           splitConvex(targetCell, tetra);
1046         }
1047       }
1048   }
1049
1050   /**
1051    * Splits the hexahedron into five tetrahedra.
1052    * This method adds five SplitterTetra objects to the vector tetra.
1053    *
1054    * @param subZone  the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1055    *                 splitting to be reused on the subzones of the GENERAL_* types of splitting
1056    */
1057   template<class MyMeshTypeT, class MyMeshTypeS>
1058   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1059   {
1060     // create tetrahedra
1061     for(int i = 0; i < 5; ++i)
1062       {
1063         const double* nodes[4];
1064         int conn[4];
1065         for(int j = 0; j < 4; ++j)
1066           {
1067             conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ];
1068             nodes[j] = getCoordsOfSubNode(conn[j]);
1069           }
1070         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1071         tetra.push_back(t);
1072       }
1073   }
1074
1075   /**
1076    * Splits the hexahedron into six tetrahedra.
1077    * This method adds six SplitterTetra objects to the vector tetra.
1078    *
1079    * @param subZone  the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1080    *                 splitting to be reused on the subzones of the GENERAL_* types of splitting
1081    */
1082   template<class MyMeshTypeT, class MyMeshTypeS>
1083   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1084   {
1085     for(int i = 0; i < 6; ++i)
1086       {
1087         const double* nodes[4];
1088         int conn[4];
1089         for(int j = 0; j < 4; ++j)
1090           {
1091             conn[j] = subZone[SPLIT_NODES_6[4*i+j]];
1092             nodes[j] = getCoordsOfSubNode(conn[j]);
1093           }
1094         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1095         tetra.push_back(t);
1096       }
1097   }
1098
1099   /**
1100    * Splits the hexahedron into 24 tetrahedra.
1101    * The splitting is done by combining the barycenter of the tetrahedron, the barycenter of each face
1102    * and the nodes of each edge of the face. This creates 6 faces * 4 edges / face = 24 tetrahedra.
1103    * The submesh nodes introduced are the barycenters of the faces and the barycenter of the cell.
1104    *
1105    */
1106   template<class MyMeshTypeT, class MyMeshTypeS>
1107   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1108   {
1109     // The two nodes of the original mesh cell used in each tetrahedron.
1110     // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
1111     // For the correspondence of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
1112
1113     // nodes to use for tetrahedron
1114     const double* nodes[4];
1115     int conn[4];
1116     // get the cell center
1117     conn[0] = 14;
1118     nodes[0] = getCoordsOfSubNode(conn[0]);
1119
1120     for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode)
1121       {
1122         // get the face center
1123         conn[1] = faceCenterNode;
1124         nodes[1] = getCoordsOfSubNode(conn[1]);
1125         for(int j = 0; j < 4; ++j)
1126           {
1127             const int row = 4*(faceCenterNode - 8) + j;
1128             conn[2] = TETRA_EDGES_GENERAL_24[2*row];
1129             conn[3] = TETRA_EDGES_GENERAL_24[2*row + 1];
1130             nodes[2] = getCoordsOfSubNode(conn[2]);
1131             nodes[3] = getCoordsOfSubNode(conn[3]);
1132
1133             SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes, conn);
1134             tetra.push_back(t);
1135           }
1136       }
1137   }
1138
1139
1140   /**
1141    * Splits the hexahedron into 48 tetrahedra.
1142    * The splitting is done by introducing the midpoints of all the edges
1143    * and the barycenter of the element as submesh nodes. The 8 hexahedral subzones thus defined
1144    * are then split into 6 tetrahedra each, as in Grandy, p. 449. The division of the subzones
1145    * is done by calling sixSplit().
1146    *
1147    */
1148   template<class MyMeshTypeT, class MyMeshTypeS>
1149   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1150   {
1151     for(int i = 0; i < 8; ++i)
1152       {
1153         sixSplit(GENERAL_48_SUBZONES+8*i,tetra);
1154       }
1155   }
1156
1157   /**
1158    * Splits the NORM_PYRA5 into 2 tetrahedra.
1159    */
1160   template<class MyMeshTypeT, class MyMeshTypeS>
1161   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1162   {
1163     static const int SPLIT_PYPA5[2][4] =
1164       {
1165         {
1166           0, 1, 2, 4
1167         },
1168         {
1169           0, 2, 3, 4
1170         }
1171       };
1172
1173     // create tetrahedra
1174     const double* nodes[4];
1175     int conn[4];
1176     for(int i = 0; i < 2; ++i)
1177       {
1178         for(int j = 0; j < 4; ++j)
1179           nodes[j] = getCoordsOfSubNode2(SPLIT_PYPA5[i][j],conn[j]);
1180         SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1181         tetra.push_back(t);
1182       }
1183   }
1184
1185   /**
1186    * Splits a convex cell into tetrahedra.
1187    */
1188   template<class MyMeshTypeT, class MyMeshTypeS>
1189   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitConvex(typename MyMeshTypeT::MyConnType targetCell,
1190                                                              typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1191   {
1192     // Each face of a cell is split into triangles and
1193     // each of triangles and a cell barycenter form a tetrahedron.
1194
1195     typedef typename MyMeshTypeT::MyConnType ConnType;
1196     const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
1197
1198     // get type of cell and nb of cell nodes
1199     NormalizedCellType normCellType=_target_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(targetCell));
1200     const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
1201     unsigned nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
1202
1203     // get nb of cell sons (faces)
1204     const ConnType* rawCellConn = _target_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _target_mesh.getConnectivityIndexPtr()[ targetCell ]);
1205     const int rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
1206     unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
1207
1208     // indices of nodes of a son
1209     static std::vector<int> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
1210     while ( allNodeIndices.size() < nbOfCellNodes )
1211       allNodeIndices.push_back( allNodeIndices.size() );
1212     std::vector<int> classicFaceNodes(4);
1213     if(cellModelCell.isQuadratic())
1214       throw INTERP_KERNEL::Exception("SplitterTetra2::splitConvex : quadratic 3D cells are not implemented yet !");
1215     int* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
1216
1217     // nodes of tetrahedron
1218     int conn[4];
1219     const double* nodes[4];
1220     nodes[3] = getCoordsOfSubNode2( nbOfCellNodes,conn[3]); // barycenter
1221
1222     for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
1223       {
1224         // get indices of son's nodes: it's just next portion of allNodeIndices for polyhedron
1225         // and some of allNodeIndices accodring to cell model for a classsic cell
1226         unsigned nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon2(ii, rawCellConn, rawNbCellNodes);
1227         if ( normCellType != NORM_POLYHED )
1228           cellModelCell.fillSonCellNodalConnectivity(ii,&allNodeIndices[0],faceNodes);
1229
1230         int nbTetra = nbFaceNodes - 2; // split polygon into nbTetra triangles
1231
1232         // create tetrahedra
1233         for(int i = 0; i < nbTetra; ++i)
1234           {
1235             nodes[0] = getCoordsOfSubNode2( faceNodes[0],  conn[0]);
1236             nodes[1] = getCoordsOfSubNode2( faceNodes[1+i],conn[1]);
1237             nodes[2] = getCoordsOfSubNode2( faceNodes[2+i],conn[2]);
1238             SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1239             tetra.push_back(t);
1240           }
1241
1242         if ( normCellType == NORM_POLYHED )
1243           faceNodes += nbFaceNodes; // go to the next face
1244       }
1245   }
1246
1247   /**
1248    * Precalculates all the nodes.
1249    * Retrieves the mesh nodes and allocates the necessary sub-mesh
1250    * nodes according to the splitting policy used.
1251    * This method is meant to be called once by the constructor.
1252    *
1253    * @param targetMesh  the target mesh
1254    * @param targetCell  the global number of the cell that the object represents, in targetMesh mode.
1255    * @param policy      the splitting policy of the object
1256    *
1257    */
1258   template<class MyMeshTypeT, class MyMeshTypeS>
1259   void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell)
1260   {
1261     // retrieve real mesh nodes
1262
1263     typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634. _node_ids.resize(8);
1264     for(int node = 0; node < nbOfNodesT ; ++node)
1265       {
1266         // calculate only normal nodes
1267         _nodes.push_back(getCoordsOfNode2(node, targetCell, targetMesh,_node_ids[node]));
1268       }
1269
1270     switch ( nbOfNodesT )
1271       {
1272       case 8:
1273
1274         // create sub-mesh nodes if needed
1275         switch(_splitting_pol)
1276           {
1277           case GENERAL_24:
1278             {
1279               for(int i = 0; i < 7; ++i)
1280                 {
1281                   double* barycenter = new double[3];
1282                   calcBarycenter(4, barycenter, &GENERAL_24_SUB_NODES[4*i]);
1283                   _nodes.push_back(barycenter);
1284                 }
1285             }
1286             break;
1287
1288           case GENERAL_48:
1289             {
1290               for(int i = 0; i < 19; ++i)
1291                 {
1292                   double* barycenter = new double[3];
1293                   calcBarycenter(2, barycenter, &GENERAL_48_SUB_NODES[2*i]);
1294                   _nodes.push_back(barycenter);
1295                 }
1296             }
1297             break;
1298
1299           default:
1300             break;
1301           }
1302
1303       case 5: // NORM_PYRA5
1304         break;
1305
1306       default: // convex 3d cell
1307         {
1308           // add barycenter of a cell
1309           std::vector<int> allIndices(nbOfNodesT);
1310           for ( int i = 0; i < nbOfNodesT; ++i ) allIndices[i] = i;
1311           double* barycenter = new double[3];
1312           calcBarycenter(nbOfNodesT, barycenter, &allIndices[0]);
1313           _nodes.push_back(barycenter);
1314         }
1315       }
1316
1317   }
1318 }
1319
1320 #endif