src/INTERP_KERNEL/TransformedTriangleInline.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 __TRANSFORMEDTRIANGLEINLINE_HXX__
  20 #define __TRANSFORMEDTRIANGLEINLINE_HXX__
  21
  22 // This file contains inline versions of some of the methods in the TransformedTriangle*.cxx files.
  23 // It replaces those methods if OPTIMIZE is defined.
  24 // NB : most of these methods have documentation in their corresponding .cxx - file.
  25
  26 // ----------------------------------------------------------------------------------
  27 //  Optimization methods. These are only defined and used if OPTIMIZE is defined.
  28 // -----------------------------------------------------------------------------------
  29
  30
  31 inline void TransformedTriangle::preCalculateTriangleSurroundsEdge()
  32 {
  33   for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
  34     {
  35       _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge);
  36     }
  37 }
  38
  39
  40 // ----------------------------------------------------------------------------------
  41 //   TransformedTriangle_math.cxx
  42 // ----------------------------------------------------------------------------------
  43
  44 inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner)
  45 {
  46   // set the three corresponding double products to 0.0
  47   static const DoubleProduct DOUBLE_PRODUCTS[12] =
  48     {
  49       C_YZ, C_ZX, C_XY, // O
  50       C_YZ, C_ZH, C_YH, // X
  51       C_ZX, C_ZH, C_XH, // Y
  52       C_XY, C_YH, C_XH  // Z
  53     };
  54
  55   for(int i = 0 ; i < 3 ; ++i) {
  56     const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i];
  57
  58     LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner);
  59     _doubleProducts[8*seg + dp] = 0.0;
  60   };
  61 }
  62
  63 inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const
  64 {
  65   return _doubleProducts[8*seg + dp];
  66 }
  67
  68 inline double TransformedTriangle::calcStableT(const TetraCorner corner) const
  69 {
  70   //   assert(_isTripleProductsCalculated);
  71   //   assert(_validTP[corner]);
  72   return _tripleProducts[corner];
  73 }
  74
  75 inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const
  76 {
  77
  78   // find the points of the triangle
  79   // 0 -> P, 1 -> Q, 2 -> R
  80   const int pt1 = seg;
  81   const int pt2 = (seg + 1) % 3;
  82
  83   // find offsets
  84   const int off1 = DP_OFFSET_1[dp];
  85   const int off2 = DP_OFFSET_2[dp];
  86
  87   return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
  88 }
  89
  90 // ----------------------------------------------------------------------------------
  91 //  TransformedTriangle_intersect.cxx
  92 // ----------------------------------------------------------------------------------
  93 inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const
  94 {
  95   return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge);
  96 }
  97
  98 inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const
  99 {
 100   return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet);
 101 }
 102
 103 inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const
 104 {
 105   return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner );
 106 }
 107
 108 inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const
 109 {
 110   const double pt[4] =
 111     {
 112       _coords[5*corner],     // x
 113       _coords[5*corner + 1], // y
 114       _coords[5*corner + 2], // z
 115       _coords[5*corner + 3]  // z
 116     };
 117
 118   for(int i = 0 ; i < 4 ; ++i)
 119     {
 120       if(pt[i] < 0.0 || pt[i] > 1.0)
 121         {
 122           return false;
 123         }
 124     }
 125   return true;
 126 }
 127
 128 inline  bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const
 129 {
 130 #if 0
 131   const double pt[4] =
 132     {
 133       _coords[5*corner],     // x
 134       _coords[5*corner + 1], // y
 135       _coords[5*corner + 2], // z
 136       _coords[5*corner + 3]  // h
 137     };
 138 #endif
 139   const double* pt = &_coords[5*corner];
 140
 141   if(pt[3] != 0.0)
 142     {
 143       return false;
 144     }
 145
 146   for(int i = 0 ; i < 3 ; ++i)
 147     {
 148       if(pt[i] < 0.0 || pt[i] > 1.0)
 149         {
 150           return false;
 151         }
 152     }
 153   return true;
 154 }
 155
 156 inline  bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const
 157 {
 158   const double x = _coords[5*corner];
 159   const double y = _coords[5*corner + 1];
 160   const double h = _coords[5*corner + 3];
 161   const double H = _coords[5*corner + 4];
 162
 163   return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
 164
 165 }
 166
 167 inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const
 168 {
 169
 170   //  assert(edge < H01);
 171
 172   // correspondance edge - triple products
 173   // for edges OX, ..., ZX (Grandy, table III)
 174   static const TetraCorner TRIPLE_PRODUCTS[12] =
 175     {
 176       X, O, // OX
 177       Y, O, // OY
 178       Z, O, // OZ
 179       X, Y, // XY
 180       Y, Z, // YZ
 181       Z, X, // ZX
 182     };
 183
 184   // Grandy, [16]
 185   const double t1 = calcStableT(TRIPLE_PRODUCTS[2*edge]);
 186   const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]);
 187
 188   //? should equality with zero use epsilon?
 189   LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 );
 190   return (t1*t2 <= 0.0) && (t1 - t2 != 0.0);
 191 }
 192
 193 inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const
 194 {
 195 #if 0
 196   const double signs[3] =
 197     {
 198       SIGN_FOR_SEG_FACET_INTERSECTION[3*facet],
 199       SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 1],
 200       SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 2]
 201     };
 202 #endif
 203
 204   const double* signs = &SIGN_FOR_SEG_FACET_INTERSECTION[3*facet];
 205   const double c1 = signs[0]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet]);
 206   const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]);
 207   const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]);
 208
 209   return (c1*c3 > 0.0) && (c2*c3 > 0.0);
 210 }
 211
 212 inline bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const
 213 {
 214   // use correspondance facet a = 0 <=> offset for coordinate a in _coords
 215   // and also correspondance segment AB => corner A
 216   const double coord1 = _coords[5*seg + facet];
 217   const double coord2 = _coords[5*( (seg + 1) % 3) + facet];
 218
 219   //? should we use epsilon-equality here in second test?
 220   LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
 221
 222   return (coord1*coord2 <= 0.0) && (coord1 != coord2);
 223 }
 224
 225 inline bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const
 226 {
 227   // get the H - coordinates
 228   const double coord1 = _coords[5*seg + 4];
 229   const double coord2 = _coords[5*( (seg + 1) % 3) + 4];
 230   //? should we use epsilon-equality here in second test?
 231   LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
 232
 233   return (coord1*coord2 <= 0.0) && (coord1 != coord2);
 234 }
 235
 236 inline bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) const
 237 {
 238   // ? There seems to be an error in Grandy -> it should be C_XY instead of C_YZ in [28].
 239   // ? I haven't really figured out why, but it seems to work.
 240   const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY);
 241
 242   LOG(6, "surface above corner " << corner << " : " << "n = " << normal << ", t = [" <<  calcTByDevelopingRow(corner, 1, false) << ", "  << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) );
 243   LOG(6, "] - stable : " << calcStableT(corner)  );
 244
 245   //? we don't care here if the triple product is "invalid", that is, the triangle does not surround one of the
 246   // edges going out from the corner (Grandy [53])
 247   if(!_validTP[corner])
 248     {
 249       return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0;
 250     }
 251   else
 252     {
 253       return ( calcStableT(corner) * normal ) >= 0.0;
 254     }
 255 }
 256
 257 inline bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const
 258 {
 259   //  assert(corner == X || corner == Y || corner == Z);
 260
 261   // double products to use for the possible corners
 262   static const DoubleProduct DP_FOR_RAY_INTERSECTION[4] =
 263     {
 264       DoubleProduct(0),        // O - only here to fill out and make indices match
 265       C_10,     // X
 266       C_01,     // Y
 267       C_XY      // Z
 268     };
 269
 270   const DoubleProduct dp = DP_FOR_RAY_INTERSECTION[corner];
 271
 272   const double cPQ = calcStableC(PQ, dp);
 273   const double cQR = calcStableC(QR, dp);
 274   const double cRP = calcStableC(RP, dp);
 275
 276   //? NB here we have no correction for precision - is this good?
 277   // Our authority Grandy says nothing
 278   LOG(5, "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
 279
 280   return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 );
 281
 282 }
 283 #endif