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();
 367
 368     double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
 369     double intersectSourceFace(const NormalizedCellType polyType,
 370                                const int polyNodesNbr,
 371                                const int *const polyNodes,
 372                                const double *const *const polyCoords,
 373                                const double dimCaracteristic,
 374                                const double precision,
 375                                std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
 376                                std::set<TriangleFaceKey>& listOfTetraFacesColinear);
 377
 378     double intersectTetra(const double** tetraCorners);
 379
 380     typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
 381
 382     void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
 383
 384     void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
 385
 386     void clearVolumesCache();
 387
 388   private:
 389     inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const;
 390     inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
 391     inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
 392     inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
 393     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
 394     inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
 395
 396     static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
 397                                 const double *const *const coordsFace, const double precision);
 398     static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
 399                                                                          const double planeConstant,
 400                                                                          const double *const p1, const double *const p2, const double *const p3,
 401                                                                          const double *const p4, const double *const p5, const double *const p6,
 402                                                                          const double dimCaracteristic, const double precision);
 403
 404     /// disallow copying
 405     SplitterTetra(const SplitterTetra& t);
 406
 407     /// disallow assignment
 408     SplitterTetra& operator=(const SplitterTetra& t);
 409
 410     // member variables
 411     /// affine transform associated with this target element
 412     TetraAffineTransform* _t;
 413
 414     /// HashMap relating node numbers to transformed nodes, used for caching
 415     HashMap< int , double* > _nodes;
 416
 417     /// HashMap relating triangular faces to calculated volume contributions, used for caching
 418     HashMap< TriangleFaceKey, double > _volumes;
 419
 420     /// reference to the source mesh
 421     const MyMeshType& _src_mesh;
 422
 423     // node id of the first node in target mesh in C mode.
 424     typename MyMeshType::MyConnType _conn[4];
 425
 426     double _coords[12];
 427
 428     /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
 429     /// 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.
 430     static const double SPARSE_TRUNCATION_LIMIT;
 431   };
 432
 433   /**
 434    * Function used to filter out elements by checking if they belong to one of the halfspaces
 435    * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
 436    * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
 437    * 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.
 438    *
 439    * @param pt        double[3] containing the coordiantes of a transformed point
 440    * @param isOutside bool[8] which indicate the results of earlier checks.
 441    */
 442   template<class MyMeshType>
 443   inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
 444   {
 445     isOutside[0] = isOutside[0] && (pt[0] < errTol);
 446     isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
 447     isOutside[2] = isOutside[2] && (pt[1] < errTol);
 448     isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
 449     isOutside[4] = isOutside[4] && (pt[2] < errTol);
 450     isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
 451     isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
 452     isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
 453   }
 454
 455   template<class MyMeshType>
 456   inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
 457   {
 458     isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
 459     isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
 460     isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
 461     isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
 462     isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
 463     isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
 464     isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
 465     isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
 466   }
 467
 468   /**
 469    * Calculates the transformed node with a given global node number.
 470    * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
 471    * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
 472    * calling.
 473    *
 474    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 475    *
 476    */
 477   template<class MyMeshType>
 478   inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
 479   {
 480     const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
 481     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 482     assert(transformedNode != 0);
 483     _t->apply(transformedNode, node);
 484     _nodes[globalNodeNum] = transformedNode;
 485   }
 486
 487
 488   /**
 489    * Calculates the transformed node with a given global node number.
 490    * Applies TetraAffineTransform * _t to it.
 491    * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
 492    * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
 493    * and are not recalculated in order to optimize the method.
 494    *
 495    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 496    *
 497    */
 498   template<class MyMeshType>
 499   inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
 500   {
 501     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 502     assert(transformedNode != 0);
 503     _t->apply(transformedNode, node);
 504     _nodes[globalNodeNum] = transformedNode;
 505   }
 506
 507   /**
 508    * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
 509    * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
 510    *
 511    * @param tri    triangle for which to calculate the volume contribution
 512    * @param key    key associated with the face
 513    */
 514   template<class MyMeshType>
 515   inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
 516   {
 517     const double vol = tri.calculateIntersectionVolume();
 518     _volumes.insert(std::make_pair(key, vol));
 519   }
 520
 521   /**
 522    * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
 523    * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
 524    *
 525    * @param tri    triangle for which to calculate the surface contribution
 526    * @param key    key associated with the face
 527    */
 528   template<class MyMeshType>
 529   inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
 530   {
 531     const double surf = tri.calculateIntersectionSurface(_t);
 532     _volumes.insert(std::make_pair(key, surf));
 533   }
 534
 535   template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
 536   class SplitterTetra2
 537   {
 538   public:
 539     SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
 540     ~SplitterTetra2();
 541     void releaseArrays();
 542     void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
 543                          typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 544     void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 545     void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 546     void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 547     void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 548     void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 549     void splitConvex(typename MyMeshTypeT::MyConnType                     targetCell,
 550                      typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 551     void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
 552     inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
 553     inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
 554     //template<int n>
 555     inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
 556   private:
 557     const MyMeshTypeT& _target_mesh;
 558     const MyMeshTypeS& _src_mesh;
 559     SplittingPolicy _splitting_pol;
 560     /// vector of pointers to double[3] containing the coordinates of the
 561     /// (sub) - nodes of split target cell
 562     std::vector<const double*> _nodes;
 563     std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
 564   };
 565
 566   /**
 567    * Calculates the barycenter of n (sub) - nodes
 568    *
 569    * @param  n  number of nodes for which to calculate barycenter
 570    * @param  barycenter  pointer to double[3] array in which to store the result
 571    * @param  pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
 572    */
 573   template<class MyMeshTypeT, class MyMeshTypeS>
 574   //template<int n>
 575   inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
 576   {
 577     barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
 578     for(int i = 0; i < n ; ++i)
 579       {
 580        const double* pt = getCoordsOfSubNode(pts[i]);
 581        barycenter[0] += pt[0];
 582        barycenter[1] += pt[1];
 583        barycenter[2] += pt[2];
 584       }
 585
 586     barycenter[0] /= n;
 587     barycenter[1] /= n;
 588     barycenter[2] /= n;
 589   }
 590
 591   /**
 592    * Accessor to the coordinates of a given (sub)-node
 593    *
 594    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 595    * @return pointer to double[3] containing the coordinates of the nodes
 596    */
 597   template<class MyMeshTypeT, class MyMeshTypeS>
 598   inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
 599   {
 600     // replace "at()" with [] for unsafe but faster access
 601     return _nodes.at(node);
 602   }
 603
 604   /**
 605    * Accessor to the coordinates of a given (sub)-node
 606    *
 607    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 608    * @param nodeId is an output that is node id in target whole mesh in C mode.
 609    * @return pointer to double[3] containing the coordinates of the nodes
 610    */
 611   template<class MyMeshTypeT, class MyMeshTypeS>
 612   const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
 613   {
 614     const double *ret=_nodes.at(node);
 615     if(node<8)
 616       nodeId=_node_ids[node];
 617     else
 618       nodeId=-1;
 619     return ret;
 620   }
 621 }
 622
 623 #endif