src/INTERP_KERNEL/SplitterTetra.hxx

   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
  20 #ifndef __SPLITTERTETRA_HXX__
  21 #define __SPLITTERTETRA_HXX__
  22
  23 #include "TransformedTriangle.hxx"
  24 #include "TetraAffineTransform.hxx"
  25 #include "InterpolationOptions.hxx"
  26 #include "InterpKernelException.hxx"
  27 #include "InterpKernelHashMap.hxx"
  28 #include "VectorUtils.hxx"
  29
  30 #include <functional>
  31 #include <vector>
  32 #include <cassert>
  33 #include <map>
  34 #include <set>
  35
  36 namespace INTERP_KERNEL
  37 {
  38   // Schema according to which the splitting is performed.
  39     // Each line represents one tetrahedron. The numbering is as follows :
  40     //
  41     //          7 ------ 6
  42     //         /|       /|
  43     //        / |      / |
  44     //       3 ------ 2  |
  45     //       |  |     |  |
  46     //       |  |     |  |
  47     //       |  4-----|- 5
  48     //       | /      | /
  49     //       0 ------ 1
  50
  51   static const int SPLIT_NODES_5[20] = /* WHY not all well oriented ???? */
  52     {
  53       0, 1, 5, 2,
  54       0, 4, 5, 7,
  55       0, 3, 7, 2,
  56       5, 6, 7, 2,
  57       0, 2, 5, 7
  58     };
  59
  60   static const int SPLIT_NODES_5_WO[20] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
  61     {
  62       0, 5, 1, 2,
  63       0, 4, 5, 7,
  64       0, 3, 7, 2,
  65       5, 7, 6, 2,
  66       0, 5, 2, 7
  67     };
  68
  69   static const int SPLIT_NODES_6[24] = /* WHY all badly oriented ???? */
  70     {
  71       0, 1, 5, 6,
  72       0, 2, 1, 6,
  73       0, 5, 4, 6,
  74       0, 4, 7, 6,
  75       0, 3, 2, 6,
  76       0, 7, 3, 6
  77     };
  78
  79   static const int SPLIT_NODES_6_WO[24] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
  80     {
  81       0, 5, 1, 6,
  82       0, 1, 2, 6,
  83       0, 4, 5, 6,
  84       0, 7, 4, 6,
  85       0, 2, 3, 6,
  86       0, 3, 7, 6
  87     };
  88
  89   // Each sub-node is the barycenter of 4 other nodes.
  90   // For the faces, these are on the orignal mesh.
  91   // For the barycenter, the four face sub-nodes are used.
  92   static const int GENERAL_24_SUB_NODES[28] =
  93     {
  94       0,1,4,5,// sub-node 8  (face)
  95       0,1,2,3,// sub-node 9  (face)
  96       0,3,4,7,// sub-node 10 (face)
  97       1,2,5,6,// sub-node 11 (face)
  98       4,5,6,7,// sub-node 12 (face)
  99       2,3,6,7,// sub-node 13 (face)
 100       8,9,10,11// sub-node 14 (cell)
 101     };
 102
 103   static const int GENERAL_24_SUB_NODES_WO[28] =
 104     {
 105       0,4,5,1,// sub-node 8  (face)
 106       0,1,2,3,// sub-node 9  (face)
 107       0,3,7,4,// sub-node 10 (face)
 108       1,5,6,2,// sub-node 11 (face)
 109       4,7,6,5,// sub-node 12 (face)
 110       2,6,7,3,// sub-node 13 (face)
 111       8,9,10,11// sub-node 14 (cell)
 112     };
 113
 114   static const int TETRA_EDGES_GENERAL_24[48] =
 115     {
 116       // face with center 8
 117       0,1,
 118       1,5,
 119       5,4,
 120       4,0,
 121       // face with center 9
 122       0,1,
 123       1,2,
 124       2,3,
 125       3,0,
 126       // face with center 10
 127       0,4,
 128       4,7,
 129       7,3,
 130       3,0,
 131       // face with center 11
 132       1,5,
 133       5,6,
 134       6,2,
 135       2,1,
 136       // face with center 12
 137       5,6,
 138       6,7,
 139       7,4,
 140       4,5,
 141       // face with center 13
 142       2,6,
 143       6,7,
 144       7,3,
 145       3,2
 146     };
 147
 148   // Each sub-node is the barycenter of two other nodes.
 149   // For the edges, these lie on the original mesh.
 150   // For the faces, these are the edge sub-nodes.
 151   // For the cell these are two face sub-nodes.
 152   static const int GENERAL_48_SUB_NODES[38] =
 153     {
 154       0,1,   // sub-node 8 (edge)
 155       0,4,   // sub-node 9 (edge)
 156       1,5,   // sub-node 10 (edge)
 157       4,5,   // sub-node 11 (edge)
 158       0,3,   // sub-node 12 (edge)
 159       1,2,   // sub-node 13 (edge)
 160       4,7,   // sub-node 14 (edge)
 161       5,6,   // sub-node 15 (edge)
 162       2,3,   // sub-node 16 (edge)
 163       3,7,   // sub-node 17 (edge)
 164       2,6,   // sub-node 18 (edge)
 165       6,7,   // sub-node 19 (edge)
 166       8,11,  // sub-node 20 (face)
 167       12,13, // sub-node 21 (face)
 168       9,17,  // sub-node 22 (face)
 169       10,18, // sub-node 23 (face)
 170       14,15, // sub-node 24 (face)
 171       16,19, // sub-node 25 (face)
 172       20,25  // sub-node 26 (cell)
 173     };
 174
 175   // Define 8 hexahedral subzones as in Grandy, p449
 176   // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell
 177   // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
 178   static const int GENERAL_48_SUBZONES[64] =
 179     {
 180       0,8,21,12,9,20,26,22,
 181       8,1,13,21,20,10,23,26,
 182       12,21,16,3,22,26,25,17,
 183       21,13,2,16,26,23,18,25,
 184       9,20,26,22,4,11,24,14,
 185       20,10,23,26,11,5,15,24,
 186       22,26,25,17,14,24,19,7,
 187       26,23,18,25,24,15,6,19
 188     };
 189
 190   static const int GENERAL_48_SUBZONES_2[64] =
 191     {
 192       0,-1,-14,-5,-2,-13,-19,-15,
 193       -1,1,-6,-14,-13,-3,-16,-19,
 194       -5,-14,-9,3,-15,-19,-18,-10,
 195       -14,-6,2,-9,-19,-16,-11,-18,
 196       -2,-13,-19,-15,4,-4,-17,-7,
 197       -13,-3,-16,-19,-4,5,-8,-17,
 198       -15,-19,-18,-10,-7,-17,-12,7,
 199       -19,-16,-11,-18,-17,-8,6,-12};
 200
 201   void SplitHexa8IntoTetras(SplittingPolicy policy, const int *nodalConnBg, const int *nodalConnEnd, const double *coords,
 202                             std::vector<int>& tetrasNodalConn, std::vector<double>& addCoords) throw(INTERP_KERNEL::Exception);
 203
 204   void SplitIntoTetras(SplittingPolicy policy, NormalizedCellType gt, const int *nodalConnBg, const int *nodalConnEnd, const double *coords,
 205                        std::vector<int>& tetrasNodalConn, std::vector<double>& addCoords) throw(INTERP_KERNEL::Exception);
 206
 207   /**
 208    * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra.
 209    *
 210    */
 211   class TriangleFaceKey
 212   {
 213   public:
 214
 215     /**
 216      * Constructor
 217      * Sorts the given nodes (so that the order in which they are passed does not matter) and
 218      * calculates a hash value for the key.
 219      *
 220      * @param node1  global number of the first node of the face
 221      * @param node2  global number of the second node of the face
 222      * @param node3  global number of the third node of the face
 223      */
 224     TriangleFaceKey(int node1, int node2, int node3)
 225     {
 226       Sort3Ints(_nodes, node1, node2, node3);
 227       _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
 228     }
 229
 230     /**
 231      * Equality comparison operator.
 232      * Compares this TriangleFaceKey object to another and determines if they represent the same face.
 233      *
 234      * @param   key  TriangleFaceKey with which to compare
 235      * @return  true if key has the same three nodes as this object, false if not
 236      */
 237     bool operator==(const TriangleFaceKey& key) const
 238     {
 239       return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
 240     }
 241
 242     /**
 243      * Less than operator.
 244      *
 245      * @param   key  TriangleFaceKey with which to compare
 246      * @return  true if this object has the three nodes less than the nodes of the key object, false if not
 247      */
 248     bool operator<(const TriangleFaceKey& key) const
 249     {
 250       for (int i = 0; i < 3; ++i)
 251         {
 252           if (_nodes[i] < key._nodes[i])
 253             {
 254               return true;
 255             }
 256           else if (_nodes[i] > key._nodes[i])
 257             {
 258               return false;
 259             }
 260         }
 261       return false;
 262     }
 263
 264     /**
 265      * Returns a hash value for the object, based on its three nodes.
 266      * This value is not unique for each face.
 267      *
 268      * @return  a hash value for the object
 269      */
 270     int hashVal() const
 271     {
 272       return _hashVal;
 273     }
 274
 275     inline static void Sort3Ints(int* sorted, int node1, int node2, int node3);
 276
 277   private:
 278     /// global numbers of the three nodes, sorted in ascending order
 279     int _nodes[3];
 280
 281     /// hash value for the object, calculated in the constructor
 282     int _hashVal;
 283   };
 284
 285   /**
 286    * Method to sort three integers in ascending order
 287    *
 288    * @param sorted  int[3] array in which to store the result
 289    * @param x1   first integer
 290    * @param x2   second integer
 291    * @param x3   third integer
 292    */
 293   inline void TriangleFaceKey::Sort3Ints(int* sorted, int x1, int x2, int x3)
 294   {
 295     if(x1 < x2)
 296       {
 297         if(x1 < x3)
 298           {
 299             // x1 is min
 300             sorted[0] = x1;
 301             sorted[1] = x2 < x3 ? x2 : x3;
 302             sorted[2] = x2 < x3 ? x3 : x2;
 303           }
 304         else
 305           {
 306             // x3, x1, x2
 307             sorted[0] = x3;
 308             sorted[1] = x1;
 309             sorted[2] = x2;
 310           }
 311       }
 312     else // x2 < x1
 313       {
 314         if(x2 < x3)
 315           {
 316             // x2 is min
 317             sorted[0] = x2;
 318             sorted[1] = x1 < x3 ? x1 : x3;
 319             sorted[2] = x1 < x3 ? x3 : x1;
 320           }
 321         else
 322           {
 323             // x3, x2, x1
 324             sorted[0] = x3;
 325             sorted[1] = x2;
 326             sorted[2] = x1;
 327           }
 328       }
 329   }
 330
 331   /**
 332    * \brief Template specialization of INTERP_KERNEL::hash<T> function object for use with a
 333    * with TriangleFaceKey as key class.
 334    *
 335    */
 336   template<> class hash<INTERP_KERNEL::TriangleFaceKey>
 337   {
 338   public:
 339     /**
 340      * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object.
 341      *
 342      * @param key  a TriangleFaceKey object
 343      * @return an integer hash value for key
 344      */
 345     int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
 346     {
 347       return key.hashVal();
 348     }
 349   };
 350 }
 351
 352 namespace INTERP_KERNEL
 353 {
 354   /**
 355    * \brief Class calculating the volume of intersection between a tetrahedral target element and
 356    * source elements with triangular or quadratilateral faces.
 357    *
 358    */
 359   template<class MyMeshType>
 360   class SplitterTetra
 361   {
 362   public:
 363
 364     SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
 365
 366     SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const int *conn = 0);
 367
 368     ~SplitterTetra();
 369
 370     double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
 371     double intersectSourceFace(const NormalizedCellType polyType,
 372                                const int polyNodesNbr,
 373                                const int *const polyNodes,
 374                                const double *const *const polyCoords,
 375                                const double dimCaracteristic,
 376                                const double precision,
 377                                std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
 378                                std::set<TriangleFaceKey>& listOfTetraFacesColinear);
 379
 380     double intersectTetra(const double** tetraCorners);
 381
 382     typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
 383
 384     void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
 385
 386     void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
 387
 388     void clearVolumesCache();
 389
 390   private:
 391     inline static void CheckIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL);
 392     inline static void CheckIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL);
 393     inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
 394     inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
 395     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
 396     inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
 397
 398     static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
 399                                 const double *const *const coordsFace, const double precision);
 400     static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
 401                                                                          const double planeConstant,
 402                                                                          const double *const p1, const double *const p2, const double *const p3,
 403                                                                          const double *const p4, const double *const p5, const double *const p6,
 404                                                                          const double dimCaracteristic, const double precision);
 405
 406     /// disallow copying
 407     SplitterTetra(const SplitterTetra& t);
 408
 409     /// disallow assignment
 410     SplitterTetra& operator=(const SplitterTetra& t);
 411
 412     // member variables
 413     /// affine transform associated with this target element
 414     TetraAffineTransform* _t;
 415
 416     /// HashMap relating node numbers to transformed nodes, used for caching
 417     HashMap< int , double* > _nodes;
 418
 419     /// HashMap relating triangular faces to calculated volume contributions, used for caching
 420     HashMap< TriangleFaceKey, double > _volumes;
 421
 422     /// reference to the source mesh
 423     const MyMeshType& _src_mesh;
 424
 425     // node id of the first node in target mesh in C mode.
 426     typename MyMeshType::MyConnType _conn[4];
 427
 428     double _coords[12];
 429
 430     /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
 431     /// 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.
 432     static const double SPARSE_TRUNCATION_LIMIT;
 433   };
 434
 435   /**
 436    * Function used to filter out elements by checking if they belong to one of the halfspaces
 437    * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
 438    * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
 439    * the corresponding array element will be true if and only if it was true when the method was called and pt lies in the halfspace.
 440    *
 441    * @param pt        double[3] containing the coordiantes of a transformed point
 442    * @param isOutside bool[8] which indicate the results of earlier checks.
 443    */
 444   template<class MyMeshType>
 445   inline void SplitterTetra<MyMeshType>::CheckIsOutside(const double* pt, bool* isOutside, const double errTol)
 446   {
 447     isOutside[0] = isOutside[0] && (pt[0] < errTol);
 448     isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
 449     isOutside[2] = isOutside[2] && (pt[1] < errTol);
 450     isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
 451     isOutside[4] = isOutside[4] && (pt[2] < errTol);
 452     isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
 453     isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
 454     isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
 455   }
 456
 457   template<class MyMeshType>
 458   inline void SplitterTetra<MyMeshType>::CheckIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol)
 459   {
 460     isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
 461     isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
 462     isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
 463     isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
 464     isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
 465     isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
 466     isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
 467     isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
 468   }
 469
 470   /**
 471    * Calculates the transformed node with a given global node number.
 472    * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
 473    * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
 474    * calling.
 475    *
 476    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 477    *
 478    */
 479   template<class MyMeshType>
 480   inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
 481   {
 482     const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
 483     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 484     assert(transformedNode != 0);
 485     _t->apply(transformedNode, node);
 486     _nodes[globalNodeNum] = transformedNode;
 487   }
 488
 489
 490   /**
 491    * Calculates the transformed node with a given global node number.
 492    * Applies TetraAffineTransform * _t to it.
 493    * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
 494    * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
 495    * and are not recalculated in order to optimize the method.
 496    *
 497    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 498    *
 499    */
 500   template<class MyMeshType>
 501   inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
 502   {
 503     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 504     assert(transformedNode != 0);
 505     _t->apply(transformedNode, node);
 506     _nodes[globalNodeNum] = transformedNode;
 507   }
 508
 509   /**
 510    * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
 511    * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
 512    *
 513    * @param tri    triangle for which to calculate the volume contribution
 514    * @param key    key associated with the face
 515    */
 516   template<class MyMeshType>
 517   inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
 518   {
 519     const double vol = tri.calculateIntersectionVolume();
 520     _volumes.insert(std::make_pair(key, vol));
 521   }
 522
 523   /**
 524    * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
 525    * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
 526    *
 527    * @param tri    triangle for which to calculate the surface contribution
 528    * @param key    key associated with the face
 529    */
 530   template<class MyMeshType>
 531   inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
 532   {
 533     const double surf = tri.calculateIntersectionSurface(_t);
 534     _volumes.insert(std::make_pair(key, surf));
 535   }
 536
 537   template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
 538   class SplitterTetra2
 539   {
 540   public:
 541     SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
 542     ~SplitterTetra2();
 543     void releaseArrays();
 544     void splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 545     void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
 546                          typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
 547     void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
 548     void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
 549     void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
 550     void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
 551     void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
 552     void splitConvex(typename MyMeshTypeT::MyConnType                     targetCell,//to suppress
 553                      typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
 554     void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);//to suppress
 555     inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);//to suppress
 556     inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);//to suppress
 557     //template<int n>
 558     inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);//to suppress
 559   private:
 560     const MyMeshTypeT& _target_mesh;
 561     const MyMeshTypeS& _src_mesh;
 562     SplittingPolicy _splitting_pol;
 563     /// vector of pointers to double[3] containing the coordinates of the
 564     /// (sub) - nodes of split target cell
 565     std::vector<const double*> _nodes;
 566     std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
 567   };
 568
 569   /**
 570    * Calculates the barycenter of n (sub) - nodes
 571    *
 572    * @param  n  number of nodes for which to calculate barycenter
 573    * @param  barycenter  pointer to double[3] array in which to store the result
 574    * @param  pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
 575    */
 576   template<class MyMeshTypeT, class MyMeshTypeS>
 577   //template<int n>
 578   inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
 579   {
 580     barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
 581     for(int i = 0; i < n ; ++i)
 582       {
 583        const double* pt = getCoordsOfSubNode(pts[i]);
 584        barycenter[0] += pt[0];
 585        barycenter[1] += pt[1];
 586        barycenter[2] += pt[2];
 587       }
 588
 589     barycenter[0] /= n;
 590     barycenter[1] /= n;
 591     barycenter[2] /= n;
 592   }
 593
 594   /**
 595    * Accessor to the coordinates of a given (sub)-node
 596    *
 597    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 598    * @return pointer to double[3] containing the coordinates of the nodes
 599    */
 600   template<class MyMeshTypeT, class MyMeshTypeS>
 601   inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
 602   {
 603     // replace "at()" with [] for unsafe but faster access
 604     return _nodes.at(node);
 605   }
 606
 607   /**
 608    * Accessor to the coordinates of a given (sub)-node
 609    *
 610    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 611    * @param nodeId is an output that is node id in target whole mesh in C mode.
 612    * @return pointer to double[3] containing the coordinates of the nodes
 613    */
 614   template<class MyMeshTypeT, class MyMeshTypeS>
 615   const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
 616   {
 617     const double *ret=_nodes.at(node);
 618     if(node<8)
 619       nodeId=_node_ids[node];
 620     else
 621       nodeId=-1;
 622     return ret;
 623   }
 624 }
 625
 626 #endif