src/INTERP_KERNEL/SplitterTetra.txx

   1 //  Copyright (C) 2007-2008  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
  31 #include <cmath>
  32 #include <cassert>
  33 #include <string>
  34 #include <sstream>
  35
  36 /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
  37 /// Since the scale is always the same in the transformed space (the target tetrahedron is unitary), this number is independent of the scale of the meshes.
  38 #define SPARSE_TRUNCATION_LIMIT 1.0e-14
  39
  40 namespace INTERP_KERNEL
  41 {
  42
  43   /**
  44    * Constructor creating object from target cell global number
  45    *
  46    * @param srcMesh     mesh containing the source elements
  47    * @param targetMesh  mesh containing the target elements
  48    * @param targetCell  global number of the target cell
  49    *
  50    */
  51   /*template<class MyMeshType>
  52   SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const MyMeshType& targetMesh, typename MyMeshType::MyConnType targetCell)
  53     : _src_mesh(srcMesh), _t(0)
  54   {
  55     // get array of corners of target tetraedron
  56     const double* tetraCorners[4];
  57     for(int i = 0 ; i < 4 ; ++i)
  58       tetraCorners[i] = getCoordsOfNode(i, targetCell, targetMesh);
  59     // create the affine transform
  60     createAffineTransform(tetraCorners);
  61     }*/
  62
  63   /*!
  64    * output is expected to be allocated with 24*sizeof(void*) in order to store the 24 tetras.
  65    * These tetras have to be deallocated.
  66    */
  67   template<class MyMeshType>
  68   void SplitterTetra<MyMeshType>::splitIntoDualCells(SplitterTetra<MyMeshType> **output)
  69   {
  70     double tmp[12];
  71     const double *tmp2[4]={tmp,tmp+3,tmp+6,tmp+9};
  72     typename MyMeshType::MyConnType conn[4]={-1,-1,-1,-1};
  73     for(int i=0;i<24;i++)
  74       {
  75         splitMySelfForDual(tmp,i,conn[0]);
  76         output[i]=new SplitterTetra<MyMeshType>(_src_mesh,tmp2,conn);
  77       }
  78   }
  79
  80   /**
  81    * Constructor creating object from the four corners of the tetrahedron.
  82    *
  83    * @param srcMesh       mesh containing the source elements
  84    * @param tetraCorners  array of four pointers to double[3] arrays containing the coordinates of the
  85    *                      corners of the tetrahedron
  86    */
  87   template<class MyMeshType>
  88   SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId)
  89     : _t(0), _src_mesh(srcMesh)
  90   {
  91     std::copy(nodesId,nodesId+4,_conn);
  92     _coords[0]=tetraCorners[0][0]; _coords[1]=tetraCorners[0][1]; _coords[2]=tetraCorners[0][2];
  93     _coords[3]=tetraCorners[1][0]; _coords[4]=tetraCorners[1][1]; _coords[5]=tetraCorners[1][2];
  94     _coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2];
  95     _coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2];
  96     // create the affine transform
  97     createAffineTransform(tetraCorners);
  98   }
  99
 100   /*!
 101    * This contructor is used to build part of 1/24th dual cell of tetraCorners.
 102    * @param i is in 0..23 included.
 103    * @param nodeId is the id of first node of this in target mesh in C mode.
 104    */
 105   /*template<class MyMeshType>
 106   SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, int i, typename MyMeshType::MyConnType nodeId)
 107     : _t(0), _src_mesh(srcMesh), _conn(nodeId)
 108   {
 109     double *newCoords[4];
 110     splitMySelfForDual(tetraCorners,newCoords,i);
 111     createAffineTransform(newCoords);
 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(hash_map< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
 125       delete[] iter->second;
 126   }
 127
 128   /*!
 129    * This method destroys the 4 pointers pointed by tetraCorners[0],tetraCorners[1],tetraCorners[2] and tetraCorners[3]
 130    * @param i is in 0..23 included.
 131    * @param output is expected to be sized of 12 in order to.
 132    */
 133   template<class MyMeshType>
 134   void SplitterTetra<MyMeshType>::splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId)
 135   {
 136     double *tmp[4];
 137     int offset=i/6;
 138     nodeId=_conn[offset];
 139     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;
 140     int caseToTreat=i%6;
 141     int case1=caseToTreat/2;
 142     int case2=caseToTreat%2;
 143     const int tab[3][2]={{1,2},{3,2},{1,3}};
 144     const int *curTab=tab[case1];
 145     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.;
 146     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[1]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.;
 147     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.;
 148     std::copy(pt1,pt1+3,output+case2*3);
 149     std::copy(pt0,pt0+3,output+(abs(case2-1))*3);
 150     std::copy(pt2,pt2+3,output+2*3);
 151     std::copy(tmp[0],tmp[0]+3,output+3*3);
 152   }
 153
 154   /**
 155    * Calculates the volume of intersection of an element in the source mesh and the target element.
 156    * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
 157    * faces of the source element are triangulated and the calculated transformation is applied
 158    * to each triangle. The algorithm of Grandy, implemented in INTERP_KERNEL::TransformedTriangle is used
 159    * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
 160    * divided by the determinant of the transformation.
 161    *
 162    * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated
 163    * with triangulated faces to avoid having to recalculate these.
 164    *
 165    * @param element      global number of the source element in C mode.
 166    */
 167   template<class MyMeshType>
 168   double SplitterTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element)
 169   {
 170     typedef typename MyMeshType::MyConnType ConnType;
 171     const NumberingPolicy numPol=MyMeshType::My_numPol;
 172     //{ could be done on outside?
 173     // check if we have planar tetra element
 174     if(_t->determinant() == 0.0)
 175       {
 176         // tetra is planar
 177         LOG(2, "Planar tetra -- volume 0");
 178         return 0.0;
 179       }
 180
 181     // get type of cell
 182     NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
 183     const CellModel& cellModelCell=CellModel::getCellModel(normCellType);
 184     unsigned nbOfNodes4Type=cellModelCell.getNumberOfNodes();
 185     // halfspace filtering
 186     bool isOutside[8] = {true, true, true, true, true, true, true, true};
 187     bool isTargetOutside = false;
 188
 189     // calculate the coordinates of the nodes
 190     int *cellNodes=new int[nbOfNodes4Type];
 191     for(int i = 0;i<nbOfNodes4Type;++i)
 192       {
 193         // we could store mapping local -> global numbers too, but not sure it is worth it
 194         const int globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
 195         cellNodes[i]=globalNodeNum;
 196         if(_nodes.find(globalNodeNum) == _nodes.end())
 197           {
 198             //for(hash_map< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
 199             //  std::cout << (*iter3).first << " ";
 200             //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
 201             calculateNode(globalNodeNum);
 202           }
 203
 204         checkIsOutside(_nodes[globalNodeNum], isOutside);
 205       }
 206
 207     // halfspace filtering check
 208     // NB : might not be beneficial for caching of triangles
 209     for(int i = 0; i < 8; ++i)
 210       {
 211         if(isOutside[i])
 212           {
 213             isTargetOutside = true;
 214           }
 215       }
 216
 217     double totalVolume = 0.0;
 218
 219     if(!isTargetOutside)
 220       {
 221         for(unsigned ii = 0 ; ii < cellModelCell.getNumberOfSons() ; ++ii)
 222           {
 223             NormalizedCellType faceType = cellModelCell.getSonType(ii);
 224             const CellModel& faceModel=CellModel::getCellModel(faceType);
 225             assert(faceModel.getDimension() == 2);
 226             int *faceNodes=new int[faceModel.getNumberOfNodes()];
 227             cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
 228             switch(faceType)
 229               {
 230               case NORM_TRI3:
 231                 {
 232                   // create the face key
 233                   TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
 234
 235                   // calculate the triangle if needed
 236                   if(_volumes.find(key) == _volumes.end())
 237                     {
 238                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
 239                       calculateVolume(tri, key);
 240                       totalVolume += _volumes[key];
 241                     } else {
 242                     // count negative as face has reversed orientation
 243                     totalVolume -= _volumes[key];
 244                   }
 245                 }
 246                 break;
 247
 248               case NORM_QUAD4:
 249
 250                 // simple triangulation of faces along a diagonal :
 251                 //
 252                 // 2 ------ 3
 253                 // |      / |
 254                 // |     /  |
 255                 // |    /   |
 256                 // |   /    |
 257                 // |  /     |
 258                 // | /      |
 259                 // 1 ------ 4
 260                 //
 261                 //? not sure if this always works
 262                 {
 263                   // calculate the triangles if needed
 264
 265                   // local nodes 1, 2, 3
 266                   TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
 267                   if(_volumes.find(key1) == _volumes.end())
 268                     {
 269                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
 270                       calculateVolume(tri, key1);
 271                       totalVolume += _volumes[key1];
 272                     } else {
 273                     // count negative as face has reversed orientation
 274                     totalVolume -= _volumes[key1];
 275                   }
 276
 277                   // local nodes 1, 3, 4
 278                   TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]);
 279                   if(_volumes.find(key2) == _volumes.end())
 280                     {
 281                       TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]);
 282                       calculateVolume(tri, key2);
 283                       totalVolume += _volumes[key2];
 284                     }
 285                   else
 286                     {
 287                       // count negative as face has reversed orientation
 288                       totalVolume -= _volumes[key2];
 289                     }
 290                 }
 291                 break;
 292
 293               default:
 294                 std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
 295                 assert(false);
 296               }
 297             delete [] faceNodes;
 298           }
 299       }
 300     delete [] cellNodes;
 301     // reset if it is very small to keep the matrix sparse
 302     // is this a good idea?
 303     if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
 304       {
 305         totalVolume = 0.0;
 306       }
 307
 308     LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
 309
 310     // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
 311     // that should be used (which is equivalent to dividing by the determinant)
 312     return std::fabs(1.0 / _t->determinant() * totalVolume) ;
 313   }
 314
 315   ////////////////////////////////////////////////////////
 316
 317   template<class MyMeshType>
 318   SplitterTetra2<MyMeshType>::SplitterTetra2(const MyMeshType& targetMesh, const MyMeshType& srcMesh, SplittingPolicy policy):_target_mesh(targetMesh),_src_mesh(srcMesh),
 319                                                                                                                               _splitting_pol(policy)
 320   {
 321   }
 322
 323   template<class MyMeshType>
 324   SplitterTetra2<MyMeshType>::~SplitterTetra2()
 325   {
 326     releaseArrays();
 327   }
 328
 329   template<class MyMeshType>
 330   void SplitterTetra2<MyMeshType>::releaseArrays()
 331   {
 332     // free potential sub-mesh nodes that have been allocated
 333     if(_nodes.size()>=8)
 334       {
 335         std::vector<const double*>::iterator iter = _nodes.begin() + 8;
 336         while(iter != _nodes.end())
 337           {
 338             delete[] *iter;
 339             ++iter;
 340           }
 341       }
 342     _nodes.clear();
 343   }
 344
 345   /*!
 346    * @param targetCell in C mode.
 347    * @param tetra is the output result tetra containers.
 348    */
 349   template<class MyMeshType>
 350   void SplitterTetra2<MyMeshType>::splitTargetCell(typename MyMeshType::MyConnType targetCell, typename MyMeshType::MyConnType nbOfNodesT,
 351                                                    typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
 352   {
 353     typedef typename MyMeshType::MyConnType ConnType;
 354     const NumberingPolicy numPol=MyMeshType::My_numPol;
 355     const int numTetra = static_cast<int>(_splitting_pol);
 356     if(nbOfNodesT==4)
 357       {
 358         _nodes.resize(8);
 359         _node_ids.resize(8);
 360         tetra.reserve(1);
 361         const double *nodes[4];
 362         int conn[4];
 363         for(int node = 0; node < 4 ; ++node)
 364           {
 365             nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
 366           }
 367         std::copy(conn,conn+4,_node_ids.begin());
 368         SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes,conn);
 369         tetra.push_back(t);
 370         return ;
 371       }
 372
 373     // pre-calculate nodes
 374     calculateSubNodes(_target_mesh, OTT<ConnType,numPol>::indFC(targetCell));
 375
 376     tetra.reserve(numTetra);
 377     _nodes.reserve(30); // we never have more than this
 378
 379     switch(_splitting_pol)
 380       {
 381       case PLANAR_FACE_5:
 382         {
 383           const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
 384           fiveSplit(subZone,tetra);
 385         }
 386         break;
 387
 388       case PLANAR_FACE_6:
 389         {
 390           const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
 391           sixSplit(subZone,tetra);
 392         }
 393         break;
 394
 395       case GENERAL_24:
 396         {
 397           calculateGeneral24Tetra(tetra);
 398         }
 399         break;
 400
 401       case GENERAL_48:
 402         {
 403           calculateGeneral48Tetra(tetra);
 404         }
 405         break;
 406       default:
 407         assert(false);
 408       }
 409   }
 410
 411   /**
 412    * Splits the hexahedron into five tetrahedra.
 413    * This method adds five SplitterTetra objects to the vector tetra.
 414    *
 415    * @param subZone  the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
 416    *                 splitting to be reused on the subzones of the GENERAL_* types of splitting
 417    */
 418   template<class MyMeshType>
 419   void SplitterTetra2<MyMeshType>::fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
 420   {
 421     // Schema according to which the splitting is performed.
 422     // Each line represents one tetrahedron. The numbering is as follows :
 423     //
 424     //          7 ------ 6
 425     //         /|       /|
 426     //        / |      / |
 427     //       3 ------ 2  |
 428     //       |  |     |  |
 429     //       |  |     |  |
 430     //       |  4-----|- 5
 431     //       | /      | /
 432     //       0 ------ 1
 433
 434
 435     static const int SPLIT_NODES_5[20] =
 436       {
 437         0, 1, 5, 2,
 438         0, 4, 5, 7,
 439         0, 3, 7, 2,
 440         5, 6, 7, 2,
 441         0, 2, 5, 7
 442       };
 443
 444     // create tetrahedra
 445     for(int i = 0; i < 5; ++i)
 446       {
 447         const double* nodes[4];
 448         int conn[4];
 449         for(int j = 0; j < 4; ++j)
 450           {
 451             nodes[j] = getCoordsOfSubNode2(subZone[ SPLIT_NODES_5[4*i+j] ],conn[j]);
 452           }
 453         SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes,conn);
 454         tetra.push_back(t);
 455       }
 456   }
 457
 458   /**
 459    * Splits the hexahedron into six tetrahedra.
 460    * This method adds six SplitterTetra objects to the vector tetra.
 461    *
 462    * @param subZone  the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
 463    *                 splitting to be reused on the subzones of the GENERAL_* types of splitting
 464    */
 465   template<class MyMeshType>
 466   void SplitterTetra2<MyMeshType>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
 467   {
 468     // Schema according to which the splitting is performed.
 469     // Each line represents one tetrahedron. The numbering is as follows :
 470     //
 471     //          7 ------ 6
 472     //         /|       /|
 473     //        / |      / |
 474     //       3 ------ 2  |
 475     //       |  |     |  |
 476     //       |  |     |  |
 477     //       |  4-----|- 5
 478     //       | /      | /
 479     //       0 ------ 1
 480
 481     static const int SPLIT_NODES_6[24] =
 482       {
 483         0, 1, 5, 6,
 484         0, 2, 1, 6,
 485         0, 5, 4, 6,
 486         0, 4, 7, 6,
 487         0, 3, 2, 6,
 488         0, 7, 3, 6
 489       };
 490
 491     for(int i = 0; i < 6; ++i)
 492       {
 493         const double* nodes[4];
 494         int conn[4];
 495         for(int j = 0; j < 4; ++j)
 496           {
 497             nodes[j] = getCoordsOfSubNode2(subZone[ SPLIT_NODES_6[4*i+j] ],conn[j]);
 498           }
 499         SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes,conn);
 500         tetra.push_back(t);
 501       }
 502   }
 503
 504   /**
 505    * Splits the hexahedron into 24 tetrahedra.
 506    * The splitting is done by combining the barycenter of the tetrahedron, the barycenter of each face
 507    * and the nodes of each edge of the face. This creates 6 faces * 4 edges / face = 24 tetrahedra.
 508    * The submesh nodes introduced are the barycenters of the faces and the barycenter of the cell.
 509    *
 510    */
 511   template<class MyMeshType>
 512   void SplitterTetra2<MyMeshType>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
 513   {
 514     // The two nodes of the original mesh cell used in each tetrahedron.
 515     // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
 516     // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
 517     static const int TETRA_EDGES[48] =
 518       {
 519         // face with center 9
 520         0,1,
 521         1,5,
 522         5,4,
 523         4,0,
 524         // face with center 10
 525         0,1,
 526         1,2,
 527         2,3,
 528         3,0,
 529         // face with center 11
 530         0,4,
 531         4,7,
 532         7,3,
 533         3,0,
 534         // face with center 12
 535         1,5,
 536         5,6,
 537         6,2,
 538         2,1,
 539         // face with center 13
 540         5,6,
 541         6,7,
 542         7,4,
 543         4,5,
 544         // face with center 14
 545         2,6,
 546         6,7,
 547         7,3,
 548         3,2
 549       };
 550
 551     // nodes to use for tetrahedron
 552     const double* nodes[4];
 553     int conn[4];
 554     // get the cell center
 555     nodes[0] = getCoordsOfSubNode2(14,conn[0]);
 556
 557     for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode)
 558       {
 559         // get the face center
 560         nodes[1] = getCoordsOfSubNode2(faceCenterNode,conn[1]);
 561         for(int j = 0; j < 4; ++j)
 562           {
 563             const int row = 4*(faceCenterNode - 9) + j;
 564             nodes[2] = getCoordsOfSubNode2(TETRA_EDGES[2*row],conn[2]);
 565             nodes[3] = getCoordsOfSubNode2(TETRA_EDGES[2*row + 1],conn[3]);
 566
 567             SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes, conn);
 568             tetra.push_back(t);
 569           }
 570       }
 571   }
 572
 573
 574   /**
 575    * Splits the hexahedron into 48 tetrahedra.
 576    * The splitting is done by introducing the midpoints of all the edges
 577    * and the barycenter of the element as submesh nodes. The 8 hexahedral subzones thus defined
 578    * are then split into 6 tetrahedra each, as in Grandy, p. 449. The division of the subzones
 579    * is done by calling sixSplit().
 580    *
 581    */
 582   template<class MyMeshType>
 583   void SplitterTetra2<MyMeshType>::calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
 584   {
 585     // Define 8 hexahedral subzones as in Grandy, p449
 586     // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell
 587     // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
 588     static const int subZones[64] =
 589       {
 590         0,8,21,12,9,20,26,22,
 591         8,1,13,21,20,10,23,26,
 592         12,21,16,3,22,26,25,17,
 593         21,13,2,16,26,23,18,25,
 594         9,20,26,22,4,11,24,14,
 595         20,10,23,26,11,5,15,24,
 596         22,26,25,17,14,24,19,7,
 597         26,23,18,25,24,15,6,19
 598       };
 599
 600     for(int i = 0; i < 8; ++i)
 601       {
 602         sixSplit(&subZones[8*i],tetra);
 603       }
 604   }
 605
 606   /**
 607    * Precalculates all the nodes.
 608    * Retrieves the mesh nodes and allocates the necessary sub-mesh
 609    * nodes according to the splitting policy used.
 610    * This method is meant to be called once by the constructor.
 611    *
 612    * @param targetMesh  the target mesh
 613    * @param targetCell  the global number of the cell that the object represents, in targetMesh mode.
 614    * @param policy      the splitting policy of the object
 615    *
 616    */
 617   template<class MyMeshType>
 618   void SplitterTetra2<MyMeshType>::calculateSubNodes(const MyMeshType& targetMesh, typename MyMeshType::MyConnType targetCell)
 619   {
 620     // retrieve real mesh nodes
 621     _node_ids.resize(8);
 622     for(int node = 0; node < 8 ; ++node)
 623       {
 624         // calculate only normal nodes
 625         _nodes.push_back(getCoordsOfNode2(node, targetCell, targetMesh,_node_ids[node]));
 626       }
 627
 628     // create sub-mesh nodes if needed
 629     switch(_splitting_pol)
 630       {
 631       case GENERAL_24:
 632         {
 633           // Each sub-node is the barycenter of 4 other nodes.
 634           // For the faces, these are on the orignal mesh.
 635           // For the barycenter, the four face sub-nodes are used.
 636           static const int GENERAL_24_SUB_NODES[28] =
 637             {
 638               0,1,4,5,// sub-node 9 (face)
 639               0,1,2,3,// sub-node 10 (face)
 640               0,3,4,7,// sub-node 11 (face)
 641               1,2,5,6,// sub-node 12 (face)
 642               4,5,6,7,// sub-node 13 (face)
 643               2,3,6,7,// sub-node 14 (face)
 644               9,10,11,12// sub-node 15 (cell)
 645             };
 646
 647           for(int i = 0; i < 7; ++i)
 648             {
 649               double* barycenter = new double[3];
 650               calcBarycenter<4>(barycenter, &GENERAL_24_SUB_NODES[4*i]);
 651               _nodes.push_back(barycenter);
 652             }
 653         }
 654         break;
 655
 656       case GENERAL_48:
 657         {
 658           // Each sub-node is the barycenter of two other nodes.
 659           // For the edges, these lie on the original mesh.
 660           // For the faces, these are the edge sub-nodes.
 661           // For the cell these are two face sub-nodes.
 662           static const int GENERAL_48_SUB_NODES[38] =
 663             {
 664               0,1,   // sub-node 9 (edge)
 665               0,4,   // sub-node 10 (edge)
 666               1,5,   // sub-node 11 (edge)
 667               4,5,   // sub-node 12 (edge)
 668               0,3,   // sub-node 13 (edge)
 669               1,2,   // sub-node 14 (edge)
 670               4,7,   // sub-node 15 (edge)
 671               5,6,   // sub-node 16 (edge)
 672               2,3,   // sub-node 17 (edge)
 673               3,7,   // sub-node 18 (edge)
 674               2,6,   // sub-node 19 (edge)
 675               6,7,   // sub-node 20 (edge)
 676               8,11,  // sub-node 21 (face)
 677               12,13, // sub-node 22 (face)
 678               9,17,  // sub-node 23 (face)
 679               10,18, // sub-node 24 (face)
 680               14,15, // sub-node 25 (face)
 681               16,19, // sub-node 26 (face)
 682               20,25  // sub-node 27 (cell)
 683             };
 684
 685           for(int i = 0; i < 19; ++i)
 686             {
 687               double* barycenter = new double[3];
 688               calcBarycenter<2>(barycenter, &GENERAL_48_SUB_NODES[2*i]);
 689               _nodes.push_back(barycenter);
 690             }
 691         }
 692         break;
 693
 694       default:
 695         break;
 696       }
 697   }
 698 }
 699
 700 #endif