src/INTERP_KERNEL/SplitterTetra.hxx

   1 // Copyright (C) 2007-2012  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 "InterpKernelHashMap.hxx"
  27 #include "VectorUtils.hxx"
  28
  29 #include <assert.h>
  30 #include <vector>
  31 #include <functional>
  32 #include <map>
  33 #include <set>
  34
  35 namespace INTERP_KERNEL
  36 {
  37   /**
  38    * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra.
  39    *
  40    */
  41   class TriangleFaceKey
  42   {
  43   public:
  44
  45     /**
  46      * Constructor
  47      * Sorts the given nodes (so that the order in which they are passed does not matter) and
  48      * calculates a hash value for the key.
  49      *
  50      * @param node1  global number of the first node of the face
  51      * @param node2  global number of the second node of the face
  52      * @param node3  global number of the third node of the face
  53      */
  54     TriangleFaceKey(int node1, int node2, int node3)
  55     {
  56       sort3Ints(_nodes, node1, node2, node3);
  57       _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
  58     }
  59
  60     /**
  61      * Equality comparison operator.
  62      * Compares this TriangleFaceKey object to another and determines if they represent the same face.
  63      *
  64      * @param   key  TriangleFaceKey with which to compare
  65      * @return  true if key has the same three nodes as this object, false if not
  66      */
  67     bool operator==(const TriangleFaceKey& key) const
  68     {
  69       return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
  70     }
  71
  72     /**
  73      * Less than operator.
  74      *
  75      * @param   key  TriangleFaceKey with which to compare
  76      * @return  true if this object has the three nodes less than the nodes of the key object, false if not
  77      */
  78     bool operator<(const TriangleFaceKey& key) const
  79     {
  80       for (int i = 0; i < 3; ++i)
  81         {
  82           if (_nodes[i] < key._nodes[i])
  83             {
  84               return true;
  85             }
  86           else if (_nodes[i] > key._nodes[i])
  87             {
  88               return false;
  89             }
  90         }
  91       return false;
  92     }
  93
  94     /**
  95      * Returns a hash value for the object, based on its three nodes.
  96      * This value is not unique for each face.
  97      *
  98      * @return  a hash value for the object
  99      */
 100     int hashVal() const
 101     {
 102       return _hashVal;
 103     }
 104
 105     inline void sort3Ints(int* sorted, int node1, int node2, int node3);
 106
 107   private:
 108     /// global numbers of the three nodes, sorted in ascending order
 109     int _nodes[3];
 110
 111     /// hash value for the object, calculated in the constructor
 112     int _hashVal;
 113   };
 114
 115   /**
 116    * Method to sort three integers in ascending order
 117    *
 118    * @param sorted  int[3] array in which to store the result
 119    * @param x1   first integer
 120    * @param x2   second integer
 121    * @param x3   third integer
 122    */
 123   inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3)
 124   {
 125     if(x1 < x2)
 126       {
 127         if(x1 < x3)
 128           {
 129             // x1 is min
 130             sorted[0] = x1;
 131             sorted[1] = x2 < x3 ? x2 : x3;
 132             sorted[2] = x2 < x3 ? x3 : x2;
 133           }
 134         else
 135           {
 136             // x3, x1, x2
 137             sorted[0] = x3;
 138             sorted[1] = x1;
 139             sorted[2] = x2;
 140           }
 141       }
 142     else // x2 < x1
 143       {
 144         if(x2 < x3)
 145           {
 146             // x2 is min
 147             sorted[0] = x2;
 148             sorted[1] = x1 < x3 ? x1 : x3;
 149             sorted[2] = x1 < x3 ? x3 : x1;
 150           }
 151         else
 152           {
 153             // x3, x2, x1
 154             sorted[0] = x3;
 155             sorted[1] = x2;
 156             sorted[2] = x1;
 157           }
 158       }
 159   }
 160
 161   /**
 162    * \brief Template specialization of INTERP_KERNEL::hash<T> function object for use with a
 163    * with TriangleFaceKey as key class.
 164    *
 165    */
 166   template<> class hash<INTERP_KERNEL::TriangleFaceKey>
 167   {
 168   public:
 169     /**
 170      * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object.
 171      *
 172      * @param key  a TriangleFaceKey object
 173      * @return an integer hash value for key
 174      */
 175     int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
 176     {
 177       return key.hashVal();
 178     }
 179   };
 180 }
 181
 182 namespace INTERP_KERNEL
 183 {
 184
 185   /**
 186    * \brief Class calculating the volume of intersection between a tetrahedral target element and
 187    * source elements with triangular or quadratilateral faces.
 188    *
 189    */
 190   template<class MyMeshType>
 191   class SplitterTetra
 192   {
 193   public:
 194
 195     SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
 196
 197     ~SplitterTetra();
 198
 199     double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
 200     double intersectSourceFace(const NormalizedCellType polyType,
 201                                const int polyNodesNbr,
 202                                const int *const polyNodes,
 203                                const double *const *const polyCoords,
 204                                const double dimCaracteristic,
 205                                const double precision,
 206                                std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
 207                                std::set<TriangleFaceKey>& listOfTetraFacesColinear);
 208
 209     double intersectTetra(const double** tetraCorners);
 210
 211     typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
 212
 213     void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
 214
 215     void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
 216
 217     void clearVolumesCache();
 218
 219   private:
 220     // member functions
 221     inline void createAffineTransform(const double** corners);
 222     inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const;
 223     inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
 224     inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
 225     inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
 226     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
 227     inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
 228
 229     static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
 230                                 const double *const *const coordsFace, const double precision);
 231     static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
 232                                                                          const double planeConstant,
 233                                                                          const double *const p1, const double *const p2, const double *const p3,
 234                                                                          const double *const p4, const double *const p5, const double *const p6,
 235                                                                          const double dimCaracteristic, const double precision);
 236
 237     /// disallow copying
 238     SplitterTetra(const SplitterTetra& t);
 239
 240     /// disallow assignment
 241     SplitterTetra& operator=(const SplitterTetra& t);
 242
 243     // member variables
 244     /// affine transform associated with this target element
 245     TetraAffineTransform* _t;
 246
 247     /// HashMap relating node numbers to transformed nodes, used for caching
 248     HashMap< int , double* > _nodes;
 249
 250     /// HashMap relating triangular faces to calculated volume contributions, used for caching
 251     HashMap< TriangleFaceKey, double
 252 // #ifdef WIN32
 253 //         , hash_compare<TriangleFaceKey,TriangleFaceKeyComparator>
 254 // #endif
 255     > _volumes;
 256
 257     /// reference to the source mesh
 258     const MyMeshType& _src_mesh;
 259
 260     // node id of the first node in target mesh in C mode.
 261     typename MyMeshType::MyConnType _conn[4];
 262
 263     double _coords[12];
 264   };
 265
 266   /**
 267    * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors
 268    *
 269    * @param corners  double*[4] array containing pointers to four double[3] arrays with the
 270    *                 coordinates of the corners of the tetrahedron
 271    */
 272   template<class MyMeshType>
 273   inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
 274   {
 275     // create AffineTransform from tetrahedron
 276     _t = new TetraAffineTransform( corners );
 277   }
 278
 279   /**
 280    * Function used to filter out elements by checking if they belong to one of the halfspaces
 281    * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
 282    * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
 283    * 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.
 284    *
 285    * @param pt        double[3] containing the coordiantes of a transformed point
 286    * @param isOutside bool[8] which indicate the results of earlier checks.
 287    */
 288   template<class MyMeshType>
 289   inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
 290   {
 291     isOutside[0] = isOutside[0] && (pt[0] < errTol);
 292     isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
 293     isOutside[2] = isOutside[2] && (pt[1] < errTol);
 294     isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
 295     isOutside[4] = isOutside[4] && (pt[2] < errTol);
 296     isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
 297     isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
 298     isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
 299   }
 300
 301   template<class MyMeshType>
 302   inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
 303   {
 304     isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
 305     isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
 306     isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
 307     isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
 308     isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
 309     isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
 310     isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
 311     isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
 312   }
 313
 314   /**
 315    * Calculates the transformed node with a given global node number.
 316    * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
 317    * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
 318    * calling.
 319    *
 320    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 321    *
 322    */
 323   template<class MyMeshType>
 324   inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
 325   {
 326     const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
 327     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 328     assert(transformedNode != 0);
 329     _t->apply(transformedNode, node);
 330     _nodes[globalNodeNum] = transformedNode;
 331   }
 332
 333
 334   /**
 335    * Calculates the transformed node with a given global node number.
 336    * Applies TetraAffineTransform * _t to it.
 337    * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
 338    * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
 339    * and are not recalculated in order to optimize the method.
 340    *
 341    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 342    *
 343    */
 344   template<class MyMeshType>
 345   inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
 346   {
 347     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 348     assert(transformedNode != 0);
 349     _t->apply(transformedNode, node);
 350     _nodes[globalNodeNum] = transformedNode;
 351   }
 352
 353   /**
 354    * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
 355    * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
 356    *
 357    * @param tri    triangle for which to calculate the volume contribution
 358    * @param key    key associated with the face
 359    */
 360   template<class MyMeshType>
 361   inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
 362   {
 363     const double vol = tri.calculateIntersectionVolume();
 364     _volumes.insert(std::make_pair(key, vol));
 365   }
 366
 367   /**
 368    * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
 369    * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
 370    *
 371    * @param tri    triangle for which to calculate the surface contribution
 372    * @param key    key associated with the face
 373    */
 374   template<class MyMeshType>
 375   inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
 376   {
 377     const double surf = tri.calculateIntersectionSurface(_t);
 378     _volumes.insert(std::make_pair(key, surf));
 379   }
 380
 381   template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
 382   class SplitterTetra2
 383   {
 384   public:
 385     SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
 386     ~SplitterTetra2();
 387     void releaseArrays();
 388     void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
 389                          typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 390     void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 391     void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 392     void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 393     void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 394     void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 395     void splitConvex(typename MyMeshTypeT::MyConnType                     targetCell,
 396                      typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 397     void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
 398     inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
 399     inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
 400     //template<int n>
 401     inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
 402   private:
 403     const MyMeshTypeT& _target_mesh;
 404     const MyMeshTypeS& _src_mesh;
 405     SplittingPolicy _splitting_pol;
 406     /// vector of pointers to double[3] containing the coordinates of the
 407     /// (sub) - nodes of split target cell
 408     std::vector<const double*> _nodes;
 409     std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
 410   };
 411
 412   /**
 413    * Calculates the barycenter of n (sub) - nodes
 414    *
 415    * @param  n  number of nodes for which to calculate barycenter
 416    * @param  barycenter  pointer to double[3] array in which to store the result
 417    * @param  pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
 418    */
 419   template<class MyMeshTypeT, class MyMeshTypeS>
 420   //template<int n>
 421   inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
 422   {
 423     barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
 424     for(int i = 0; i < n ; ++i)
 425       {
 426        const double* pt = getCoordsOfSubNode(pts[i]);
 427        barycenter[0] += pt[0];
 428        barycenter[1] += pt[1];
 429        barycenter[2] += pt[2];
 430       }
 431
 432     barycenter[0] /= n;
 433     barycenter[1] /= n;
 434     barycenter[2] /= n;
 435   }
 436
 437   /**
 438    * Accessor to the coordinates of a given (sub)-node
 439    *
 440    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 441    * @return pointer to double[3] containing the coordinates of the nodes
 442    */
 443   template<class MyMeshTypeT, class MyMeshTypeS>
 444   inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
 445   {
 446     // replace "at()" with [] for unsafe but faster access
 447     return _nodes.at(node);
 448   }
 449
 450   /**
 451    * Accessor to the coordinates of a given (sub)-node
 452    *
 453    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 454    * @param nodeId is an output that is node id in target whole mesh in C mode.
 455    * @return pointer to double[3] containing the coordinates of the nodes
 456    */
 457   template<class MyMeshTypeT, class MyMeshTypeS>
 458   const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
 459   {
 460     const double *ret=_nodes.at(node);
 461     if(node<8)
 462       nodeId=_node_ids[node];
 463     else
 464       nodeId=-1;
 465     return ret;
 466   }
 467 }
 468
 469 #endif