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