src/INTERP_KERNEL/SplitterTetra.txx

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