src/INTERP_KERNEL/SplitterTetra.hxx

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