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 static 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     inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const;
 320     inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
 321     inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
 322     inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
 323     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
 324     inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
 325
 326     static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
 327                                 const double *const *const coordsFace, const double precision);
 328     static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
 329                                                                          const double planeConstant,
 330                                                                          const double *const p1, const double *const p2, const double *const p3,
 331                                                                          const double *const p4, const double *const p5, const double *const p6,
 332                                                                          const double dimCaracteristic, const double precision);
 333
 334     /// disallow copying
 335     SplitterTetra(const SplitterTetra& t);
 336
 337     /// disallow assignment
 338     SplitterTetra& operator=(const SplitterTetra& t);
 339
 340     // member variables
 341     /// affine transform associated with this target element
 342     TetraAffineTransform* _t;
 343
 344     /// HashMap relating node numbers to transformed nodes, used for caching
 345     HashMap< int , double* > _nodes;
 346
 347     /// HashMap relating triangular faces to calculated volume contributions, used for caching
 348     HashMap< TriangleFaceKey, double > _volumes;
 349
 350     /// reference to the source mesh
 351     const MyMeshType& _src_mesh;
 352
 353     // node id of the first node in target mesh in C mode.
 354     typename MyMeshType::MyConnType _conn[4];
 355
 356     double _coords[12];
 357
 358     /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
 359     /// 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.
 360     static const double SPARSE_TRUNCATION_LIMIT;
 361   };
 362
 363   /**
 364    * Function used to filter out elements by checking if they belong to one of the halfspaces
 365    * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
 366    * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
 367    * 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.
 368    *
 369    * @param pt        double[3] containing the coordiantes of a transformed point
 370    * @param isOutside bool[8] which indicate the results of earlier checks.
 371    */
 372   template<class MyMeshType>
 373   inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
 374   {
 375     isOutside[0] = isOutside[0] && (pt[0] < errTol);
 376     isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
 377     isOutside[2] = isOutside[2] && (pt[1] < errTol);
 378     isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
 379     isOutside[4] = isOutside[4] && (pt[2] < errTol);
 380     isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
 381     isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
 382     isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
 383   }
 384
 385   template<class MyMeshType>
 386   inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
 387   {
 388     isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
 389     isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
 390     isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
 391     isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
 392     isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
 393     isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
 394     isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
 395     isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
 396   }
 397
 398   /**
 399    * Calculates the transformed node with a given global node number.
 400    * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
 401    * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
 402    * calling.
 403    *
 404    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 405    *
 406    */
 407   template<class MyMeshType>
 408   inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
 409   {
 410     const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
 411     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 412     assert(transformedNode != 0);
 413     _t->apply(transformedNode, node);
 414     _nodes[globalNodeNum] = transformedNode;
 415   }
 416
 417
 418   /**
 419    * Calculates the transformed node with a given global node number.
 420    * Applies TetraAffineTransform * _t to it.
 421    * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
 422    * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
 423    * and are not recalculated in order to optimize the method.
 424    *
 425    * @param globalNodeNum  global node number of the node in the mesh _src_mesh
 426    *
 427    */
 428   template<class MyMeshType>
 429   inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
 430   {
 431     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
 432     assert(transformedNode != 0);
 433     _t->apply(transformedNode, node);
 434     _nodes[globalNodeNum] = transformedNode;
 435   }
 436
 437   /**
 438    * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
 439    * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
 440    *
 441    * @param tri    triangle for which to calculate the volume contribution
 442    * @param key    key associated with the face
 443    */
 444   template<class MyMeshType>
 445   inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
 446   {
 447     const double vol = tri.calculateIntersectionVolume();
 448     _volumes.insert(std::make_pair(key, vol));
 449   }
 450
 451   /**
 452    * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
 453    * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
 454    *
 455    * @param tri    triangle for which to calculate the surface contribution
 456    * @param key    key associated with the face
 457    */
 458   template<class MyMeshType>
 459   inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
 460   {
 461     const double surf = tri.calculateIntersectionSurface(_t);
 462     _volumes.insert(std::make_pair(key, surf));
 463   }
 464
 465   template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
 466   class SplitterTetra2
 467   {
 468   public:
 469     SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
 470     ~SplitterTetra2();
 471     void releaseArrays();
 472     void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
 473                          typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 474     void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 475     void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 476     void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 477     void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 478     void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 479     void splitConvex(typename MyMeshTypeT::MyConnType                     targetCell,
 480                      typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
 481     void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
 482     inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
 483     inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
 484     //template<int n>
 485     inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
 486   private:
 487     const MyMeshTypeT& _target_mesh;
 488     const MyMeshTypeS& _src_mesh;
 489     SplittingPolicy _splitting_pol;
 490     /// vector of pointers to double[3] containing the coordinates of the
 491     /// (sub) - nodes of split target cell
 492     std::vector<const double*> _nodes;
 493     std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
 494   };
 495
 496   /**
 497    * Calculates the barycenter of n (sub) - nodes
 498    *
 499    * @param  n  number of nodes for which to calculate barycenter
 500    * @param  barycenter  pointer to double[3] array in which to store the result
 501    * @param  pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
 502    */
 503   template<class MyMeshTypeT, class MyMeshTypeS>
 504   //template<int n>
 505   inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
 506   {
 507     barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
 508     for(int i = 0; i < n ; ++i)
 509       {
 510        const double* pt = getCoordsOfSubNode(pts[i]);
 511        barycenter[0] += pt[0];
 512        barycenter[1] += pt[1];
 513        barycenter[2] += pt[2];
 514       }
 515
 516     barycenter[0] /= n;
 517     barycenter[1] /= n;
 518     barycenter[2] /= n;
 519   }
 520
 521   /**
 522    * Accessor to the coordinates of a given (sub)-node
 523    *
 524    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 525    * @return pointer to double[3] containing the coordinates of the nodes
 526    */
 527   template<class MyMeshTypeT, class MyMeshTypeS>
 528   inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
 529   {
 530     // replace "at()" with [] for unsafe but faster access
 531     return _nodes.at(node);
 532   }
 533
 534   /**
 535    * Accessor to the coordinates of a given (sub)-node
 536    *
 537    * @param  node  local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
 538    * @param nodeId is an output that is node id in target whole mesh in C mode.
 539    * @return pointer to double[3] containing the coordinates of the nodes
 540    */
 541   template<class MyMeshTypeT, class MyMeshTypeS>
 542   const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
 543   {
 544     const double *ret=_nodes.at(node);
 545     if(node<8)
 546       nodeId=_node_ids[node];
 547     else
 548       nodeId=-1;
 549     return ret;
 550   }
 551 }
 552
 553 #endif