src/INTERP_KERNEL/SplitterTetra.hxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 #ifndef __SPLITTERTETRA_HXX__
  20 #define __SPLITTERTETRA_HXX__
  21
  22 #include "TransformedTriangle.hxx"
  23 #include "TetraAffineTransform.hxx"
  24
  25 #include <assert.h>
  26 #include <vector>
  27 #include <functional>
  28 #include <map>
  29 #ifdef WIN32
  30 # include <hash_map>
  31 #else
  32 # include <ext/hash_map>
  33 #endif
  34
  35 #ifndef WIN32
  36 using __gnu_cxx::hash_map;
  37 #else
  38 using stdext::hash_map;
  39 using stdext::hash_compare;
  40 #endif
  41
  42 namespace INTERP_KERNEL
  43 {
  44   /**
  45    * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra.
  46    *
  47    */
  48   class TriangleFaceKey
  49   {
  50   public:
  51
  52     /**
  53      * Constructor
  54      * Sorts the given nodes (so that the order in which they are passed does not matter) and
  55      * calculates a hash value for the key.
  56      *
  57      * @param node1  global number of the first node of the face
  58      * @param node2  global number of the second node of the face
  59      * @param node3  global number of the third node of the face
  60      */
  61     TriangleFaceKey(int node1, int node2, int node3)
  62     {
  63       sort3Ints(_nodes, node1, node2, node3);
  64       _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
  65     }
  66
  67     /**
  68      * Equality comparison operator.
  69      * Compares this TriangleFaceKey object to another and determines if they represent the same face.
  70      *
  71      * @param   key  TriangleFaceKey with which to compare
  72      * @return  true if key has the same three nodes as this object, false if not
  73      */
  74     bool operator==(const TriangleFaceKey& key) const
  75     {
  76       return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
  77     }
  78
  79     /**
  80      * Returns a hash value for the object, based on its three nodes.
  81      * This value is not unique for each face.
  82      *
  83      * @return  a hash value for the object
  84      */
  85     int hashVal() const
  86     {
  87       return _hashVal;
  88     }
  89
  90 #ifdef WIN32
  91     operator size_t () const
  92     {
  93       return _hashVal;
  94     }
  95 #endif
  96
  97     inline void sort3Ints(int* sorted, int node1, int node2, int node3);
  98
  99   private:
 100     /// global numbers of the three nodes, sorted in ascending order
 101     int _nodes[3];
 102
 103     /// hash value for the object, calculated in the constructor
 104     int _hashVal;
 105   };
 106
 107   /**
 108    * Method to sort three integers in ascending order
 109    *
 110    * @param sorted  int[3] array in which to store the result
 111    * @param x1   first integer
 112    * @param x2   second integer
 113    * @param x3   third integer
 114    */
 115   inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3)
 116   {
 117     if(x1 < x2)
 118       {
 119         if(x1 < x3)
 120           {
 121             // x1 is min
 122             sorted[0] = x1;
 123             sorted[1] = x2 < x3 ? x2 : x3;
 124             sorted[2] = x2 < x3 ? x3 : x2;
 125           }
 126         else
 127           {
 128             // x3, x1, x2
 129             sorted[0] = x3;
 130             sorted[1] = x1;
 131             sorted[2] = x2;
 132           }
 133       }
 134     else // x2 < x1
 135       {
 136         if(x2 < x3)
 137           {
 138             // x2 is min
 139             sorted[0] = x2;
 140             sorted[1] = x1 < x3 ? x1 : x3;
 141             sorted[2] = x1 < x3 ? x3 : x1;
 142           }
 143         else
 144           {
 145             // x3, x2, x1
 146             sorted[0] = x3;
 147             sorted[1] = x2;
 148             sorted[2] = x1;
 149           }
 150       }
 151   }
 152 }
 153 #ifndef WIN32
 154 namespace __gnu_cxx
 155 {
 156
 157
 158   /**
 159    * \brief Template specialization of __gnu_cxx::hash<T> function object for use with a __gnu_cxx::hash_map
 160    * with TriangleFaceKey as key class.
 161    *
 162    */
 163   template<>
 164   class hash<INTERP_KERNEL::TriangleFaceKey>
 165
 166   {
 167   public:
 168     /**
 169      * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object.
 170      *
 171      * @param key  a TriangleFaceKey object
 172      * @return an integer hash value for key
 173      */
 174     int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
 175     {
 176       return key.hashVal();
 177     }
 178   };
 179 }
 180 #else
 181   struct TriangleFaceKeyComparator
 182   {
 183     bool operator()(const INTERP_KERNEL::TriangleFaceKey& key1,
 184                     const INTERP_KERNEL::TriangleFaceKey& key2 ) const
 185     {
 186       return key1.hashVal() < key2.hashVal();
 187     }
 188   };
 189 #endif
 190
 191 namespace INTERP_KERNEL
 192 {
 193
 194   /**
 195    * \brief Class calculating the volume of intersection between a tetrahedral target element and
 196    * source elements with triangular or quadratilateral faces.
 197    *
 198    */
 199   template<class MyMeshType>
 200   class SplitterTetra
 201   {
 202   public:
 203
 204     SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
 205
 206     ~SplitterTetra();
 207
 208     double intersectSourceCell(typename MyMeshType::MyConnType srcCell);
 209
 210     typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
 211
 212     void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
 213
 214     void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
 215
 216   private:
 217     // member functions
 218     inline void createAffineTransform(const double** corners);
 219     inline void checkIsOutside(const double* pt, bool* isOutside) const;
 220     inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
 221     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
 222
 223
 224     /// disallow copying
 225     SplitterTetra(const SplitterTetra& t);
 226
 227     /// disallow assignment
 228     SplitterTetra& operator=(const SplitterTetra& t);
 229
 230     // member variables
 231     /// affine transform associated with this target element
 232     TetraAffineTransform* _t;
 233
 234     /// hash_map relating node numbers to transformed nodes, used for caching
 235     hash_map< int , double* > _nodes;
 236
 237     /// hash_map relating triangular faces to calculated volume contributions, used for caching
 238     hash_map< TriangleFaceKey, double
 239 #ifdef WIN32
 240         , hash_compare<TriangleFaceKey,TriangleFaceKeyComparator>
 241 #endif
 242     > _volumes;
 243
 244     /// reference to the source mesh
 245     const MyMeshType& _src_mesh;
 246
 247     // node id of the first node in target mesh in C mode.
 248     typename MyMeshType::MyConnType _conn[4];
 249
 250     double _coords[12];
 251   };
 252
 253   /**
 254    * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors
 255    *
 256    * @param corners  double*[4] array containing pointers to four double[3] arrays with the
 257    *                 coordinates of the corners of the tetrahedron
 258    */
 259   template<class MyMeshType>
 260   inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
 261   {
 262     // create AffineTransform from tetrahedron
 263     _t = new TetraAffineTransform( corners );
 264   }
 265
 266   /**
 267    * Function used to filter out elements by checking if they belong to one of the halfspaces
 268    * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
 269    * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
 270    * 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.
 271    *
 272    * @param pt        double[3] containing the coordiantes of a transformed point
 273    * @param isOutside bool[8] which indicate the results of earlier checks.
 274    */
 275   template<class MyMeshType>
 276   inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside) const
 277   {
 278     isOutside[0] = isOutside[0] && (pt[0] <= 0.0);
 279     isOutside[1] = isOutside[1] && (pt[0] >= 1.0);
 280     isOutside[2] = isOutside[2] && (pt[1] <= 0.0);
 281     isOutside[3] = isOutside[3] && (pt[1] >= 1.0);
 282     isOutside[4] = isOutside[4] && (pt[2] <= 0.0);
 283     isOutside[5] = isOutside[5] && (pt[2] >= 1.0);
 284     isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] <= 0.0);
 285     isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0);
 286   }
 287
 288   /**
 289    * Calculates the transformed node with a given global node number.
 290    * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
 291    * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
 292    * calling.
 293    *
 294    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 295    *
 296    */
 297   template<class MyMeshType>
 298   inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
 299   {
 300     const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
 301     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 302     assert(transformedNode != 0);
 303     _t->apply(transformedNode, node);
 304     _nodes[globalNodeNum] = transformedNode;
 305   }
 306
 307   /**
 308    * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
 309    * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
 310    *
 311    * @param tri    triangle for which to calculate the volume contribution
 312    * @param key    key associated with the face
 313    */
 314   template<class MyMeshType>
 315   inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
 316   {
 317     const double vol = tri.calculateIntersectionVolume();
 318     _volumes.insert(std::make_pair(key, vol));
 319   }
 320
 321   template<class MyMeshType>
 322   class SplitterTetra2
 323   {
 324   public:
 325     SplitterTetra2(const MyMeshType& targetMesh, const MyMeshType& srcMesh, SplittingPolicy policy);
 326     ~SplitterTetra2();
 327     void releaseArrays();
 328     void splitTargetCell(typename MyMeshType::MyConnType targetCell, typename MyMeshType::MyConnType nbOfNodesT,
 329                          typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
 330     void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
 331     void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
 332     void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
 333     void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
 334     void calculateSubNodes(const MyMeshType& targetMesh, typename MyMeshType::MyConnType targetCell);
 335     inline const double* getCoordsOfSubNode(typename MyMeshType::MyConnType node);
 336     inline const double* getCoordsOfSubNode2(typename MyMeshType::MyConnType node, typename MyMeshType::MyConnType& nodeId);
 337     template<int n>
 338     inline void calcBarycenter(double* barycenter, const typename MyMeshType::MyConnType* pts);
 339   private:
 340     const MyMeshType& _target_mesh;
 341     const MyMeshType& _src_mesh;
 342     SplittingPolicy _splitting_pol;
 343     /// vector of pointers to double[3] containing the coordinates of the
 344     /// (sub) - nodes of split target cell
 345     std::vector<const double*> _nodes;
 346     std::vector<typename MyMeshType::MyConnType> _node_ids;
 347   };
 348
 349   /**
 350    * Calculates the barycenter of n (sub) - nodes
 351    *
 352    * @param  n  number of nodes for which to calculate barycenter
 353    * @param  barycenter  pointer to double[3] array in which to store the result
 354    * @param  pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
 355    */
 356   template<class MyMeshType>
 357   template<int n>
 358   inline void SplitterTetra2<MyMeshType>::calcBarycenter(double* barycenter, const typename MyMeshType::MyConnType* pts)
 359   {
 360     barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
 361     for(int i = 0; i < n ; ++i)
 362       {
 363        const double* pt = getCoordsOfSubNode(pts[i]);
 364        barycenter[0] += pt[0];
 365        barycenter[1] += pt[1];
 366        barycenter[2] += pt[2];
 367       }
 368
 369     barycenter[0] /= n;
 370     barycenter[1] /= n;
 371     barycenter[2] /= n;
 372   }
 373
 374   /**
 375    * Accessor to the coordinates of a given (sub)-node
 376    *
 377    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 378    * @return pointer to double[3] containing the coordinates of the nodes
 379    */
 380   template<class MyMeshType>
 381   inline const double* SplitterTetra2<MyMeshType>::getCoordsOfSubNode(typename MyMeshType::MyConnType node)
 382   {
 383     // replace "at()" with [] for unsafe but faster access
 384     return _nodes.at(node);
 385   }
 386
 387   /**
 388    * Accessor to the coordinates of a given (sub)-node
 389    *
 390    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 391    * @param nodeId is an output that is node id in target whole mesh in C mode.
 392    * @return pointer to double[3] containing the coordinates of the nodes
 393    */
 394   template<class MyMeshType>
 395   const double* SplitterTetra2<MyMeshType>::getCoordsOfSubNode2(typename MyMeshType::MyConnType node, typename MyMeshType::MyConnType& nodeId)
 396   {
 397     const double *ret=_nodes.at(node);
 398     if(node<8)
 399       nodeId=_node_ids[node];
 400     else
 401       nodeId=-1;
 402     return ret;
 403   }
 404 }
 405
 406 #endif