1 // Copyright (C) 2007-2015 CEA/DEN, EDF R&D
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, or (at your option) any later version.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #ifndef __TRANSFORMEDTRIANGLEINLINE_HXX__
21 #define __TRANSFORMEDTRIANGLEINLINE_HXX__
23 // This file contains inline versions of some of the methods in the TransformedTriangle*.cxx files.
24 // It replaces those methods if OPTIMIZE is defined.
25 // NB : most of these methods have documentation in their corresponding .cxx - file.
27 // ----------------------------------------------------------------------------------
28 // Optimization methods. These are only defined and used if OPTIMIZE is defined.
29 // -----------------------------------------------------------------------------------
32 inline void TransformedTriangle::preCalculateTriangleSurroundsEdge()
34 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
36 _triangleSurroundsEdgeCache[edge] = testTriangleSurroundsEdge(edge);
41 // ----------------------------------------------------------------------------------
42 // TransformedTriangle_math.cxx
43 // ----------------------------------------------------------------------------------
45 inline void TransformedTriangle::resetDoubleProducts(const TriSegment seg, const TetraCorner corner)
47 // set the three corresponding double products to 0.0
48 static const DoubleProduct DOUBLE_PRODUCTS[12] =
50 C_YZ, C_ZX, C_XY, // O
51 C_YZ, C_ZH, C_YH, // X
52 C_ZX, C_ZH, C_XH, // Y
56 for(int i = 0 ; i < 3 ; ++i) {
57 const DoubleProduct dp = DOUBLE_PRODUCTS[3*corner + i];
59 LOG(6, std::endl << "resetting inconsistent dp :" << dp << " for corner " << corner);
60 _doubleProducts[8*seg + dp] = 0.0;
64 inline double TransformedTriangle::calcStableC(const TriSegment seg, const DoubleProduct dp) const
66 return _doubleProducts[8*seg + dp];
69 inline double TransformedTriangle::calcStableT(const TetraCorner corner) const
71 // assert(_isTripleProductsCalculated);
72 // assert(_validTP[corner]);
73 return _tripleProducts[corner];
76 inline double TransformedTriangle::calcUnstableC(const TriSegment seg, const DoubleProduct dp) const
79 // find the points of the triangle
80 // 0 -> P, 1 -> Q, 2 -> R
82 const int pt2 = (seg + 1) % 3;
85 const int off1 = DP_OFFSET_1[dp];
86 const int off2 = DP_OFFSET_2[dp];
88 return _coords[5*pt1 + off1] * _coords[5*pt2 + off2] - _coords[5*pt1 + off2] * _coords[5*pt2 + off1];
91 // ----------------------------------------------------------------------------------
92 // TransformedTriangle_intersect.cxx
93 // ----------------------------------------------------------------------------------
94 inline bool TransformedTriangle::testSurfaceEdgeIntersection(const TetraEdge edge) const
96 return _triangleSurroundsEdgeCache[edge] && testEdgeIntersectsTriangle(edge);
99 inline bool TransformedTriangle::testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const
101 return testFacetSurroundsSegment(seg, facet) && testSegmentIntersectsFacet(seg, facet);
104 inline bool TransformedTriangle::testSurfaceRayIntersection(const TetraCorner corner) const
106 return testTriangleSurroundsRay( corner ) && testSurfaceAboveCorner( corner );
109 inline bool TransformedTriangle::testCornerInTetrahedron(const TriCorner corner) const
113 _coords[5*corner], // x
114 _coords[5*corner + 1], // y
115 _coords[5*corner + 2], // z
116 _coords[5*corner + 3] // z
119 for(int i = 0 ; i < 4 ; ++i)
121 if(pt[i] < 0.0 || pt[i] > 1.0)
129 inline bool TransformedTriangle::testCornerOnXYZFacet(const TriCorner corner) const
134 _coords[5*corner], // x
135 _coords[5*corner + 1], // y
136 _coords[5*corner + 2], // z
137 _coords[5*corner + 3] // h
140 const double* pt = &_coords[5*corner];
147 for(int i = 0 ; i < 3 ; ++i)
149 if(pt[i] < 0.0 || pt[i] > 1.0)
157 inline bool TransformedTriangle::testCornerAboveXYZFacet(const TriCorner corner) const
159 const double x = _coords[5*corner];
160 const double y = _coords[5*corner + 1];
161 const double h = _coords[5*corner + 3];
162 const double H = _coords[5*corner + 4];
164 return h < 0.0 && H >= 0.0 && x >= 0.0 && y >= 0.0;
168 inline bool TransformedTriangle::testEdgeIntersectsTriangle(const TetraEdge edge) const
171 // assert(edge < H01);
173 // correspondance edge - triple products
174 // for edges OX, ..., ZX (Grandy, table III)
175 static const TetraCorner TRIPLE_PRODUCTS[12] =
186 const double t1 = calcStableT(TRIPLE_PRODUCTS[2*edge]);
187 const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]);
189 //? should equality with zero use epsilon?
190 LOG(5, "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 );
191 return (t1*t2 <= 0.0) && (t1 - t2 != 0.0);
194 inline bool TransformedTriangle::testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const
197 const double signs[3] =
199 SIGN_FOR_SEG_FACET_INTERSECTION[3*facet],
200 SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 1],
201 SIGN_FOR_SEG_FACET_INTERSECTION[3*facet + 2]
205 const double* signs = &SIGN_FOR_SEG_FACET_INTERSECTION[3*facet];
206 const double c1 = signs[0]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet]);
207 const double c2 = signs[1]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]);
208 const double c3 = signs[2]*calcStableC(seg, DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]);
210 return (c1*c3 > 0.0) && (c2*c3 > 0.0);
213 inline bool TransformedTriangle::testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const
215 // use correspondance facet a = 0 <=> offset for coordinate a in _coords
216 // and also correspondance segment AB => corner A
217 const double coord1 = _coords[5*seg + facet];
218 const double coord2 = _coords[5*( (seg + 1) % 3) + facet];
220 //? should we use epsilon-equality here in second test?
221 LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
223 return (coord1*coord2 <= 0.0) && (coord1 != coord2);
226 inline bool TransformedTriangle::testSegmentIntersectsHPlane(const TriSegment seg) const
228 // get the H - coordinates
229 const double coord1 = _coords[5*seg + 4];
230 const double coord2 = _coords[5*( (seg + 1) % 3) + 4];
231 //? should we use epsilon-equality here in second test?
232 LOG(5, "coord1 : " << coord1 << " coord2 : " << coord2 );
234 return (coord1*coord2 <= 0.0) && (coord1 != coord2);
237 inline bool TransformedTriangle::testSurfaceAboveCorner(const TetraCorner corner) const
239 // ? There seems to be an error in Grandy -> it should be C_XY instead of C_YZ in [28].
240 // ? I haven't really figured out why, but it seems to work.
241 const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY);
243 LOG(6, "surface above corner " << corner << " : " << "n = " << normal << ", t = [" << calcTByDevelopingRow(corner, 1, false) << ", " << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) );
244 LOG(6, "] - stable : " << calcStableT(corner) );
246 //? we don't care here if the triple product is "invalid", that is, the triangle does not surround one of the
247 // edges going out from the corner (Grandy [53])
248 if(!_validTP[corner])
250 return ( calcTByDevelopingRow(corner, 1, false) * normal ) >= 0.0;
254 return ( calcStableT(corner) * normal ) >= 0.0;
258 inline bool TransformedTriangle::testTriangleSurroundsRay(const TetraCorner corner) const
260 // assert(corner == X || corner == Y || corner == Z);
262 // double products to use for the possible corners
263 static const DoubleProduct DP_FOR_RAY_INTERSECTION[4] =
265 DoubleProduct(0), // O - only here to fill out and make indices match
271 const DoubleProduct dp = DP_FOR_RAY_INTERSECTION[corner];
273 const double cPQ = calcStableC(PQ, dp);
274 const double cQR = calcStableC(QR, dp);
275 const double cRP = calcStableC(RP, dp);
277 //? NB here we have no correction for precision - is this good?
278 // Our authority Grandy says nothing
279 LOG(5, "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" );
281 return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 );