1 // Copyright (C) 2007-2021 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)
64 * OVERVIEW of how the class works : (details can be found in the documentation of each method)
67 * The constructor takes as arguments three pointers to double[3] vectors holding the transformed
68 * coordinates of the corners of the triangle. It copies their coordinates and then proceeds to pre-calculating certain
69 * entities used in the intersection calculation : the double products, triple products and the values of the function E
72 * calculateIntersectionVolume() :
73 * This is the only method in the public interface. It calculates the volume under the intersection polygons
74 * between the triangle and the unit tetrahedron, as described in Grandy, pp. 435-447. It does this by first calculating the
75 * intersection polygons A and B, with the method calculateIntersectionPolygons(). It then calculates the barycenter of each
76 * polygon in calculatePolygonBarycenter(), and sorts their points in a circular order around the barycenter in
77 * sortIntersecionPolygon(). The sorting is done with STL sort, using the order defined in the class
78 * ProjectedCentralCircularSortOrder. The volume under each polygon is then calculated with calculateVolumeUnderPolygon(), which
79 * implements formula [34] in Grandy.
81 * calculateIntersectionPolygons() :
82 * This method goes through all the possible ways in which the triangle can intersect the tetrahedron and tests for these
83 * types of intersections in accordance with the formulas described in Grandy. These tests are implemented in the test* - methods.
84 * The formulas in the article are stated for one case each only, while the calculation must take into account all cases.
85 * To this end, a number of tables, implemented as static const arrays of different types, are used. The tables
86 * mainly contain values of the different enumeration types described at the beginning of the class interface. For example,
87 * the formula Grandy gives for the segment-halfstrip intersection tests ([30]) is for use with the halfstrip above the zx edge.
88 * For the other two halfstrips (above the xy and yz edges), other double products are used, which
89 * are stored in the table DP_FOR_HALFSTRIP_INTERSECTION. This allows us to treat
90 * all the edges equally, avoiding switch() - statements. It is the careful choice of order of the enumeration types that makes this
91 * possible. Notably, there is a correspondence between the TetraEdge type and the DoubleProduct type (see Grandy, table III) that
92 * is used throughout the code, permitting statements such as DoubleProduct(some_edge) to work.
93 * When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it
94 * is not known directly. It is then added to the polygon A and/or B as necessary.
97 * If OPTIMIZE is defined, a large number of methods will be prefixed with inline and some optimizations concerning the tests
98 * with zero double products will be used.
100 class INTERPKERNEL_EXPORT TransformedTriangle
106 friend class INTERP_TEST::TransformedTriangleTest;
107 friend class INTERP_TEST::TransformedTriangleIntersectTest;
109 * Enumerations representing the different geometric elements of the unit tetrahedron
110 * and the triangle. The end element, NO_* gives the number of elements in the enumeration
111 * and can be used as end element in loops.
114 /// Corners of tetrahedron
115 enum TetraCorner { O = 0, X, Y, Z, NO_TET_CORNER };
117 /// Edges of tetrahedron
118 enum TetraEdge { OX = 0, OY, OZ, XY, YZ, ZX, H01, H10, NO_TET_EDGE };
120 /// Facets (faces) of tetrahedron
121 enum TetraFacet { OYZ = 0, OZX, OXY, XYZ, NO_TET_FACET };
123 /// Corners of triangle
124 enum TriCorner { P = 0, Q, R, NO_TRI_CORNER };
126 /// Segments (edges) of triangle
127 enum TriSegment { PQ = 0, QR, RP, NO_TRI_SEGMENT };
129 /// Intersection polygons
130 enum IntersectionPolygon{ A = 0, B, NO_INTERSECTION_POLYGONS };
133 /// NB : order corresponds to TetraEdges (Grandy, table III)
134 enum DoubleProduct { C_YZ = 0, C_ZX, C_XY, C_ZH, C_XH, C_YH, C_01, C_10, NO_DP };
136 TransformedTriangle(double* p, double* q, double* r);
137 ~TransformedTriangle();
139 double calculateIntersectionVolume();
140 double calculateIntersectionSurface(TetraAffineTransform* tat);
142 void dumpCoords() const;
144 // Queries of member values used by UnitTetraIntersectionBary
146 const double* getCorner(TriCorner corner) const { return _coords + 5*corner; }
148 const std::vector<double*>& getPolygonA() const { return _polygonA; }
150 double getVolume() const { return _volume; }
154 TransformedTriangle() { }
156 // ----------------------------------------------------------------------------------
157 // High-level methods called directly by calculateIntersectionVolume()
158 // ----------------------------------------------------------------------------------
159 void calculateIntersectionAndProjectionPolygons();
161 void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter);
163 void sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter);
165 double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter);
167 // ----------------------------------------------------------------------------------
168 // High-level methods called directly by calculateIntersectionSurface()
169 // ----------------------------------------------------------------------------------
170 void calculateIntersectionPolygon();
172 double calculateSurfacePolygon();
174 // ----------------------------------------------------------------------------------
175 // Detection of degenerate triangles
176 // ----------------------------------------------------------------------------------
178 bool isTriangleInPlaneOfFacet(const TetraFacet facet) const;
180 bool isTriangleParallelToFacet(const TetraFacet facet) const;
182 int isTriangleInclinedToFacet(const TetraFacet facet) const;
184 bool isTriangleBelowTetraeder() const;
186 // ----------------------------------------------------------------------------------
187 // Intersection test methods and intersection point calculations
188 // ----------------------------------------------------------------------------------
190 inline bool testSurfaceEdgeIntersection(const TetraEdge edge) const;
192 void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const;
194 inline bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const;
196 void calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const;
198 bool testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const;
200 void calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const ;
202 bool testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const ;
204 inline bool testSurfaceRayIntersection(const TetraCorner corner) const;
206 bool testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg);
208 void calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const;
210 bool testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const;
212 inline bool testCornerInTetrahedron(const TriCorner corner) const;
214 inline bool testCornerOnXYZFacet(const TriCorner corner) const;
216 inline bool testCornerAboveXYZFacet(const TriCorner corner) const;
218 // ----------------------------------------------------------------------------------
219 // Utility methods used in intersection tests
220 // ----------------------------------------------------------------------------------
222 bool testTriangleSurroundsEdge(const TetraEdge edge) const;
224 inline bool testEdgeIntersectsTriangle(const TetraEdge edge) const;
226 inline bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const;
228 inline bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const;
230 bool testSegmentIntersectsHPlane(const TriSegment seg) const;
232 bool testSurfaceAboveCorner(const TetraCorner corner) const;
234 bool testTriangleSurroundsRay(const TetraCorner corner) const;
236 // ----------------------------------------------------------------------------------
237 // Double and triple product calculations
238 // ----------------------------------------------------------------------------------
240 void resetNearZeroCoordinates();
242 bool areDoubleProductsConsistent(const TriSegment seg) const;
244 void preCalculateDoubleProducts(void);
246 inline void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
248 double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const;
250 void preCalculateTripleProducts(void);
252 double calculateAngleEdgeTriangle(const TetraEdge edge) const;
254 inline double calcStableC(const TriSegment seg, const DoubleProduct dp) const;
256 inline double calcStableT(const TetraCorner corner) const;
258 inline double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const;
260 double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const;
262 // ----------------------------------------------------------------------------------
264 // ----------------------------------------------------------------------------------
267 /// Array holding the coordinates of the triangle's three corners
269 /// [ 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 ]
272 /// Flag showing whether the double products have been calculated yet
273 bool _is_double_products_calculated;
275 /// Flag showing whether the triple products have been calculated yet
276 bool _is_triple_products_calculated;
278 /// Array containing the 24 double products.
279 /// order : c^PQ_YZ, ... ,cPQ_10, ... c^QR_YZ, ... c^RP_YZ
280 /// following order in enumeration DoubleProduct
281 double _doubleProducts[24];
283 /// Array containing the 4 triple products.
284 /// order : t_O, t_X, t_Y, t_Z
285 double _tripleProducts[4];
287 /// Vector holding the points of the intersection polygon A.
288 /// these points are allocated in calculateIntersectionPolygons() and liberated in the destructor
289 std::vector<double*> _polygonA;
291 /// Vector holding the points of the intersection polygon B.
292 /// These points are allocated in calculateIntersectionPolygons() and liberated in the destructor
293 std::vector<double*> _polygonB;
295 /// Array holding the coordinates of the barycenter of the polygon A
296 /// This point is calculated in calculatePolygonBarycenter
297 double _barycenterA[3];
299 /// Array holding the coordinates of the barycenter of the polygon B
300 /// This point is calculated in calculatePolygonBarycenter
301 //double _barycenterB[3];
303 /// Array of flags indicating which of the four triple products have been correctly calculated.
304 /// Used for asserts in debug mode
307 /// calculated volume for use of UnitTetraIntersectionBary
311 * Calls TransformedTriangle::testTriangleSurroundsEdge for edges OX to ZX and stores the result in
312 * member variable array_triangleSurroundsEdgeCache.
315 void preCalculateTriangleSurroundsEdge();
317 /// Array holding results of the test testTriangleSurroundsEdge() for all the edges.
318 /// These are calculated in preCalculateTriangleSurroundsEdge().
319 bool _triangleSurroundsEdgeCache[NO_TET_EDGE];
321 // ----------------------------------------------------------------------------------
323 // ----------------------------------------------------------------------------------
325 // offsets : 0 -> x, 1 -> y, 2 -> z, 3 -> h, 4 -> H
326 // corresponds to order of double products in DoubleProduct
327 // so that offset[C_*] gives the right coordinate
328 static const int DP_OFFSET_1[8];
329 static const int DP_OFFSET_2[8];
331 // the coordinates used in the expansion of triple products by a given row
332 // in constellation (corner, row-1)
333 // (0,1,2,3) <=> (x,y,z,h)
334 static const int COORDINATE_FOR_DETERMINANT_EXPANSION[12];
336 // contains the edge of the double product used when
337 // expanding the triple product determinant associated with each corner
339 static const DoubleProduct DP_FOR_DETERMINANT_EXPANSION[12];
341 // values used to decide how imprecise the double products
342 // should be to set them to 0.0
343 static const long double MACH_EPS; // machine epsilon
344 static const long double MULT_PREC_F; // precision of multiplications (Grandy : f)
345 static const long double THRESHOLD_F; // threshold for zeroing (Grandy : F/f)
347 static const double TRIPLE_PRODUCT_ANGLE_THRESHOLD;
349 // correspondence facet - double product
351 static const DoubleProduct DP_FOR_SEG_FACET_INTERSECTION[12];
353 // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION
354 static const double SIGN_FOR_SEG_FACET_INTERSECTION[12];
356 // coordinates of corners of tetrahedron
357 static const double COORDS_TET_CORNER[12];
359 // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
360 // for the calculation of the coordinates (x,y,z) of the intersection points
361 // for Segment-Facet and Segment-Edge intersections
362 static const int DP_INDEX[12];
364 // correspondence edge - corners
365 static const TetraCorner CORNERS_FOR_EDGE[12];
367 // correspondence edge - facets
368 // facets shared by each edge
369 static const TetraFacet FACET_FOR_EDGE[12];
371 // correspondence edge - corners
372 static const TetraEdge EDGES_FOR_CORNER[12];
374 // double products used in segment-halfstrip test
375 static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12];
377 // double products used in segment - ray test
378 static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21];
382 // include definitions of inline methods
384 #include "TransformedTriangleInline.hxx"