src/INTERP_KERNEL/SplitterTetra.hxx

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