1 // Copyright (C) 2007-2014 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 __TRANSFORMED_TRIANGLE_HXX__
21 #define __TRANSFORMED_TRIANGLE_HXX__
23 #include "INTERPKERNELDefines.hxx"
28 // 1 - overview of algorithm + volume result
29 // 2 - algorithm detail
30 // 3 - intersection polygon results detail
31 // 4 - intersection polygon search detail
32 // higher -> misc. gory details of calculation
37 #pragma warning(disable:4251)
42 class TransformedTriangleTest;
43 class TransformedTriangleIntersectTest;
47 namespace INTERP_KERNEL
49 class TetraAffineTransform;
51 /** \class TransformedTriangle
52 * \brief Class representing one of the faces of the triangulated source polyhedron after having been transformed
53 * with the affine transform that takes the target tetrahedron to the unit tetrahedron. It contains the
54 * logic for calculating the volume of intersection between the triangle and the unit tetrahedron.
56 * \see TransformedTriangle.hxx
58 * Reference : J. Grandy, "Conservative Remapping and Region Overlays by Intersecting Arbitrary Polyhedra",
59 * Journal of Computational Physics (1999)
63 /** \file TransformedTriangle.hxx
65 * OVERVIEW of how the class works : (details can be found in the documentation of each method)
68 * The constructor takes as arguments three pointers to double[3] vectors holding the transformed
69 * coordinates of the corners of the triangle. It copies their coordinates and then proceeds to pre-calculating certain
70 * entities used in the intersection calculation : the double products, triple products and the values of the function E
73 * calculateIntersectionVolume() :
74 * This is the only method in the public interface. It calculates the volume under the intersection polygons
75 * between the triangle and the unit tetrahedron, as described in Grandy, pp. 435-447. It does this by first calculating the
76 * intersection polygons A and B, with the method calculateIntersectionPolygons(). It then calculates the barycenter of each
77 * polygon in calculatePolygonBarycenter(), and sorts their points in a circular order around the barycenter in
78 * sortIntersecionPolygon(). The sorting is done with STL sort, using the order defined in the class
79 * ProjectedCentralCircularSortOrder. The volume under each polygon is then calculated with calculateVolumeUnderPolygon(), which
80 * implements formula [34] in Grandy.
82 * calculateIntersectionPolygons() :
83 * This method goes through all the possible ways in which the triangle can intersect the tetrahedron and tests for these
84 * types of intersections in accordance with the formulas described in Grandy. These tests are implemented in the test* - methods.
85 * The formulas in the article are stated for one case each only, while the calculation must take into account all cases.
86 * To this end, a number of tables, implemented as static const arrays of different types, are used. The tables
87 * mainly contain values of the different enumeration types described at the beginning of the class interface. For example,
88 * the formula Grandy gives for the segment-halfstrip intersection tests ([30]) is for use with the halfstrip above the zx edge.
89 * For the other two halfstrips (above the xy and yz edges), other double products are used, which
90 * are stored in the table DP_FOR_HALFSTRIP_INTERSECTION. This allows us to treat
91 * all the edges equally, avoiding switch() - statements. It is the careful choice of order of the enumeration types that makes this
92 * possible. Notably, there is a correspondance between the TetraEdge type and the DoubleProduct type (see Grandy, table III) that
93 * is used throughout the code, permitting statements such as DoubleProduct(some_edge) to work.
94 * When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it
95 * is not known directly. It is then added to the polygon A and/or B as necessary.
98 * If OPTIMIZE is defined, a large number of methods will be prefixed with inline and some optimizations concerning the tests
99 * with zero double products will be used.
101 class INTERPKERNEL_EXPORT TransformedTriangle
107 friend class INTERP_TEST::TransformedTriangleTest;
108 friend class INTERP_TEST::TransformedTriangleIntersectTest;
110 * Enumerations representing the different geometric elements of the unit tetrahedron
111 * and the triangle. The end element, NO_* gives the number of elements in the enumeration
112 * and can be used as end element in loops.
115 /// Corners of tetrahedron
116 enum TetraCorner { O = 0, X, Y, Z, NO_TET_CORNER };
118 /// Edges of tetrahedron
119 enum TetraEdge { OX = 0, OY, OZ, XY, YZ, ZX, H01, H10, NO_TET_EDGE };
121 /// Facets (faces) of tetrahedron
122 enum TetraFacet { OYZ = 0, OZX, OXY, XYZ, NO_TET_FACET };
124 /// Corners of triangle
125 enum TriCorner { P = 0, Q, R, NO_TRI_CORNER };
127 /// Segments (edges) of triangle
128 enum TriSegment { PQ = 0, QR, RP, NO_TRI_SEGMENT };
130 /// Intersection polygons
131 enum IntersectionPolygon{ A = 0, B, NO_INTERSECTION_POLYGONS };
134 /// NB : order corresponds to TetraEdges (Grandy, table III)
135 enum DoubleProduct { C_YZ = 0, C_ZX, C_XY, C_ZH, C_XH, C_YH, C_01, C_10, NO_DP };
137 TransformedTriangle(double* p, double* q, double* r);
138 ~TransformedTriangle();
140 double calculateIntersectionVolume();
141 double calculateIntersectionSurface(TetraAffineTransform* tat);
143 void dumpCoords() const;
145 // Queries of member values used by UnitTetraIntersectionBary
147 const double* getCorner(TriCorner corner) const { return _coords + 5*corner; }
149 const std::vector<double*>& getPolygonA() const { return _polygonA; }
151 double getVolume() const { return _volume; }
155 TransformedTriangle() { }
157 // ----------------------------------------------------------------------------------
158 // High-level methods called directly by calculateIntersectionVolume()
159 // ----------------------------------------------------------------------------------
160 void calculateIntersectionAndProjectionPolygons();
162 void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter);
164 void sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter);
166 double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter);
168 // ----------------------------------------------------------------------------------
169 // High-level methods called directly by calculateIntersectionSurface()
170 // ----------------------------------------------------------------------------------
171 void calculateIntersectionPolygon();
173 double calculateSurfacePolygon();
175 // ----------------------------------------------------------------------------------
176 // Detection of degenerate triangles
177 // ----------------------------------------------------------------------------------
179 bool isTriangleInPlaneOfFacet(const TetraFacet facet) const;
181 bool isTriangleParallelToFacet(const TetraFacet facet) const;
183 int isTriangleInclinedToFacet(const TetraFacet facet) const;
185 bool isTriangleBelowTetraeder() const;
187 // ----------------------------------------------------------------------------------
188 // Intersection test methods and intersection point calculations
189 // ----------------------------------------------------------------------------------
191 inline bool testSurfaceEdgeIntersection(const TetraEdge edge) const;
193 void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const;
195 inline bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const;
197 void calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const;
199 bool testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const;
201 void calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const ;
203 bool testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const ;
205 inline bool testSurfaceRayIntersection(const TetraCorner corner) const;
207 bool testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg);
209 void calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const;
211 bool testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const;
213 inline bool testCornerInTetrahedron(const TriCorner corner) const;
215 inline bool testCornerOnXYZFacet(const TriCorner corner) const;
217 inline bool testCornerAboveXYZFacet(const TriCorner corner) const;
219 // ----------------------------------------------------------------------------------
220 // Utility methods used in intersection tests
221 // ----------------------------------------------------------------------------------
223 bool testTriangleSurroundsEdge(const TetraEdge edge) const;
225 inline bool testEdgeIntersectsTriangle(const TetraEdge edge) const;
227 inline bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const;
229 inline bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const;
231 bool testSegmentIntersectsHPlane(const TriSegment seg) const;
233 bool testSurfaceAboveCorner(const TetraCorner corner) const;
235 bool testTriangleSurroundsRay(const TetraCorner corner) const;
237 // ----------------------------------------------------------------------------------
238 // Double and triple product calculations
239 // ----------------------------------------------------------------------------------
241 void resetNearZeroCoordinates();
243 bool areDoubleProductsConsistent(const TriSegment seg) const;
245 void preCalculateDoubleProducts(void);
247 inline void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
249 double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const;
251 void preCalculateTripleProducts(void);
253 double calculateAngleEdgeTriangle(const TetraEdge edge) const;
255 inline double calcStableC(const TriSegment seg, const DoubleProduct dp) const;
257 inline double calcStableT(const TetraCorner corner) const;
259 inline double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const;
261 double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const;
263 // ----------------------------------------------------------------------------------
265 // ----------------------------------------------------------------------------------
268 /// Array holding the coordinates of the triangle's three corners
270 /// [ p_x, p_y, p_z, p_h, p_H, q_x, q_y, q_z, q_h, q_H, r_x, r_y, r_z, r_h, r_H ]
273 /// Flag showing whether the double products have been calculated yet
274 bool _is_double_products_calculated;
276 /// Flag showing whether the triple products have been calculated yet
277 bool _is_triple_products_calculated;
279 /// Array containing the 24 double products.
280 /// order : c^PQ_YZ, ... ,cPQ_10, ... c^QR_YZ, ... c^RP_YZ
281 /// following order in enumeration DoubleProduct
282 double _doubleProducts[24];
284 /// Array containing the 4 triple products.
285 /// order : t_O, t_X, t_Y, t_Z
286 double _tripleProducts[4];
288 /// Vector holding the points of the intersection polygon A.
289 /// these points are allocated in calculateIntersectionPolygons() and liberated in the destructor
290 std::vector<double*> _polygonA;
292 /// Vector holding the points of the intersection polygon B.
293 /// These points are allocated in calculateIntersectionPolygons() and liberated in the destructor
294 std::vector<double*> _polygonB;
296 /// Array holding the coordinates of the barycenter of the polygon A
297 /// This point is calculated in calculatePolygonBarycenter
298 double _barycenterA[3];
300 /// Array holding the coordinates of the barycenter of the polygon B
301 /// This point is calculated in calculatePolygonBarycenter
302 //double _barycenterB[3];
304 /// Array of flags indicating which of the four triple products have been correctly calculated.
305 /// Used for asserts in debug mode
308 /// calculated volume for use of UnitTetraIntersectionBary
312 * Calls TransformedTriangle::testTriangleSurroundsEdge for edges OX to ZX and stores the result in
313 * member variable array_triangleSurroundsEdgeCache.
316 void preCalculateTriangleSurroundsEdge();
318 /// Array holding results of the test testTriangleSurroundsEdge() for all the edges.
319 /// These are calculated in preCalculateTriangleSurroundsEdge().
320 bool _triangleSurroundsEdgeCache[NO_TET_EDGE];
322 // ----------------------------------------------------------------------------------
324 // ----------------------------------------------------------------------------------
326 // offsets : 0 -> x, 1 -> y, 2 -> z, 3 -> h, 4 -> H
327 // corresponds to order of double products in DoubleProduct
328 // so that offset[C_*] gives the right coordinate
329 static const int DP_OFFSET_1[8];
330 static const int DP_OFFSET_2[8];
332 // the coordinates used in the expansion of triple products by a given row
333 // in constellation (corner, row-1)
334 // (0,1,2,3) <=> (x,y,z,h)
335 static const int COORDINATE_FOR_DETERMINANT_EXPANSION[12];
337 // contains the edge of the double product used when
338 // expanding the triple product determinant associated with each corner
340 static const DoubleProduct DP_FOR_DETERMINANT_EXPANSION[12];
342 // values used to decide how imprecise the double products
343 // should be to set them to 0.0
344 static const long double MACH_EPS; // machine epsilon
345 static const long double MULT_PREC_F; // precision of multiplications (Grandy : f)
346 static const long double THRESHOLD_F; // threshold for zeroing (Grandy : F/f)
348 static const double TRIPLE_PRODUCT_ANGLE_THRESHOLD;
350 // correspondance facet - double product
352 static const DoubleProduct DP_FOR_SEG_FACET_INTERSECTION[12];
354 // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION
355 static const double SIGN_FOR_SEG_FACET_INTERSECTION[12];
357 // coordinates of corners of tetrahedron
358 static const double COORDS_TET_CORNER[12];
360 // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
361 // for the calculation of the coordinates (x,y,z) of the intersection points
362 // for Segment-Facet and Segment-Edge intersections
363 static const int DP_INDEX[12];
365 // correspondance edge - corners
366 static const TetraCorner CORNERS_FOR_EDGE[12];
368 // correspondance edge - facets
369 // facets shared by each edge
370 static const TetraFacet FACET_FOR_EDGE[12];
372 // correspondance edge - corners
373 static const TetraEdge EDGES_FOR_CORNER[12];
375 // double products used in segment-halfstrip test
376 static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12];
378 // double products used in segment - ray test
379 static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21];
383 // include definitions of inline methods
385 #include "TransformedTriangleInline.hxx"