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