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
 365   /**
 366    * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors
 367    *
 368    * @param corners  double*[4] array containing pointers to four double[3] arrays with the
 369    *                 coordinates of the corners of the tetrahedron
 370    */
 371   template<class MyMeshType>
 372   inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
 373   {
 374     // create AffineTransform from tetrahedron
 375     _t = new TetraAffineTransform( corners );
 376   }
 377
 378   /**
 379    * Function used to filter out elements by checking if they belong to one of the halfspaces
 380    * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
 381    * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
 382    * 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.
 383    *
 384    * @param pt        double[3] containing the coordiantes of a transformed point
 385    * @param isOutside bool[8] which indicate the results of earlier checks.
 386    */
 387   template<class MyMeshType>
 388   inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
 389   {
 390     isOutside[0] = isOutside[0] && (pt[0] < errTol);
 391     isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
 392     isOutside[2] = isOutside[2] && (pt[1] < errTol);
 393     isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
 394     isOutside[4] = isOutside[4] && (pt[2] < errTol);
 395     isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
 396     isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
 397     isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
 398   }
 399
 400   template<class MyMeshType>
 401   inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
 402   {
 403     isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
 404     isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
 405     isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
 406     isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
 407     isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
 408     isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
 409     isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
 410     isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
 411   }
 412
 413   /**
 414    * Calculates the transformed node with a given global node number.
 415    * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
 416    * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
 417    * calling.
 418    *
 419    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 420    *
 421    */
 422   template<class MyMeshType>
 423   inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
 424   {
 425     const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
 426     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 427     assert(transformedNode != 0);
 428     _t->apply(transformedNode, node);
 429     _nodes[globalNodeNum] = transformedNode;
 430   }
 431
 432
 433   /**
 434    * Calculates the transformed node with a given global node number.
 435    * Applies TetraAffineTransform * _t to it.
 436    * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
 437    * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
 438    * and are not recalculated in order to optimize the method.
 439    *
 440    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 441    *
 442    */
 443   template<class MyMeshType>
 444   inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
 445   {
 446     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 447     assert(transformedNode != 0);
 448     _t->apply(transformedNode, node);
 449     _nodes[globalNodeNum] = transformedNode;
 450   }
 451
 452   /**
 453    * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
 454    * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
 455    *
 456    * @param tri    triangle for which to calculate the volume contribution
 457    * @param key    key associated with the face
 458    */
 459   template<class MyMeshType>
 460   inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
 461   {
 462     const double vol = tri.calculateIntersectionVolume();
 463     _volumes.insert(std::make_pair(key, vol));
 464   }
 465
 466   /**
 467    * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
 468    * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
 469    *
 470    * @param tri    triangle for which to calculate the surface contribution
 471    * @param key    key associated with the face
 472    */
 473   template<class MyMeshType>
 474   inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
 475   {
 476     const double surf = tri.calculateIntersectionSurface(_t);
 477     _volumes.insert(std::make_pair(key, surf));
 478   }
 479
 480   template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
 481   class SplitterTetra2
 482   {
 483   public:
 484     SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
 485     ~SplitterTetra2();
 486     void releaseArrays();
 487     void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
 488                          typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 489     void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 490     void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 491     void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 492     void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 493     void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 494     void splitConvex(typename MyMeshTypeT::MyConnType                     targetCell,
 495                      typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 496     void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
 497     inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
 498     inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
 499     //template<int n>
 500     inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
 501   private:
 502     const MyMeshTypeT& _target_mesh;
 503     const MyMeshTypeS& _src_mesh;
 504     SplittingPolicy _splitting_pol;
 505     /// vector of pointers to double[3] containing the coordinates of the
 506     /// (sub) - nodes of split target cell
 507     std::vector<const double*> _nodes;
 508     std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
 509   };
 510
 511   /**
 512    * Calculates the barycenter of n (sub) - nodes
 513    *
 514    * @param  n  number of nodes for which to calculate barycenter
 515    * @param  barycenter  pointer to double[3] array in which to store the result
 516    * @param  pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
 517    */
 518   template<class MyMeshTypeT, class MyMeshTypeS>
 519   //template<int n>
 520   inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
 521   {
 522     barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
 523     for(int i = 0; i < n ; ++i)
 524       {
 525        const double* pt = getCoordsOfSubNode(pts[i]);
 526        barycenter[0] += pt[0];
 527        barycenter[1] += pt[1];
 528        barycenter[2] += pt[2];
 529       }
 530
 531     barycenter[0] /= n;
 532     barycenter[1] /= n;
 533     barycenter[2] /= n;
 534   }
 535
 536   /**
 537    * Accessor to the coordinates of a given (sub)-node
 538    *
 539    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 540    * @return pointer to double[3] containing the coordinates of the nodes
 541    */
 542   template<class MyMeshTypeT, class MyMeshTypeS>
 543   inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
 544   {
 545     // replace "at()" with [] for unsafe but faster access
 546     return _nodes.at(node);
 547   }
 548
 549   /**
 550    * Accessor to the coordinates of a given (sub)-node
 551    *
 552    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 553    * @param nodeId is an output that is node id in target whole mesh in C mode.
 554    * @return pointer to double[3] containing the coordinates of the nodes
 555    */
 556   template<class MyMeshTypeT, class MyMeshTypeS>
 557   const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
 558   {
 559     const double *ret=_nodes.at(node);
 560     if(node<8)
 561       nodeId=_node_ids[node];
 562     else
 563       nodeId=-1;
 564     return ret;
 565   }
 566 }
 567
 568 #endif