1 // Copyright (C) 2007-2008 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.
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
19 #ifndef __TRANSFORMED_TRIANGLE_HXX__
20 #define __TRANSFORMED_TRIANGLE_HXX__
22 #include "INTERPKERNELDefines.hxx"
27 // 1 - overview of algorithm + volume result
28 // 2 - algorithm detail
29 // 3 - intersection polygon results detail
30 // 4 - intersection polygon search detail
31 // higher -> misc. gory details of calculation
37 class TransformedTriangleTest;
38 class TransformedTriangleIntersectTest;
42 namespace INTERP_KERNEL
45 /** \class TransformedTriangle
46 * \brief Class representing one of the faces of the triangulated source polyhedron after having been transformed
47 * with the affine transform that takes the target tetrahedron to the unit tetrahedron. It contains the
48 * logic for calculating the volume of intersection between the triangle and the unit tetrahedron.
50 * \see TransformedTriangle.hxx
52 * Reference : J. Grandy, "Conservative Remapping and Region Overlays by Intersecting Arbitrary Polyhedra",
53 * Journal of Computational Physics (1999)
57 /** \file TransformedTriangle.hxx
59 * OVERVIEW of how the class works : (details can be found in the documentation of each method)
62 * The constructor takes as arguments three pointers to double[3] vectors holding the transformed
63 * coordinates of the corners of the triangle. It copies their coordinates and then proceeds to pre-calculating certain
64 * entities used in the intersection calculation : the double products, triple products and the values of the function E
67 * calculateIntersectionVolume() :
68 * This is the only method in the public interface. It calculates the volume under the intersection polygons
69 * between the triangle and the unit tetrahedron, as described in Grandy, pp. 435-447. It does this by first calculating the
70 * intersection polygons A and B, with the method calculateIntersectionPolygons(). It then calculates the barycenter of each
71 * polygon in calculatePolygonBarycenter(), and sorts their points in a circular order around the barycenter in
72 * sortIntersecionPolygon(). The sorting is done with STL sort, using the order defined in the class
73 * ProjectedCentralCircularSortOrder. The volume under each polygon is then calculated with calculateVolumeUnderPolygon(), which
74 * implements formula [34] in Grandy.
76 * calculateIntersectionPolygons() :
77 * This method goes through all the possible ways in which the triangle can intersect the tetrahedron and tests for these
78 * types of intersections in accordance with the formulas described in Grandy. These tests are implemented in the test* - methods.
79 * The formulas in the article are stated for one case each only, while the calculation must take into account all cases.
80 * To this end, a number of tables, implemented as static const arrays of different types, are used. The tables
81 * mainly contain values of the different enumeration types described at the beginning of the class interface. For example,
82 * the formula Grandy gives for the segment-halfstrip intersection tests ([30]) is for use with the halfstrip above the zx edge.
83 * For the other two halfstrips (above the xy and yz edges), other double products are used, which
84 * are stored in the table DP_FOR_HALFSTRIP_INTERSECTION. This allows us to treat
85 * all the edges equally, avoiding switch() - statements. It is the careful choice of order of the enumeration types that makes this
86 * possible. Notably, there is a correspondance between the TetraEdge type and the DoubleProduct type (see Grandy, table III) that
87 * is used throughout the code, permitting statements such as DoubleProduct(some_edge) to work.
88 * When an intersection point has been detected it is calculated with a corresponding calc* - method in the cases where it
89 * is not known directly. It is then added to the polygon A and/or B as necessary.
92 * If OPTIMIZE is defined, a large number of methods will be prefixed with inline and some optimizations concerning the tests
93 * with zero double products will be used.
95 class INTERPKERNEL_EXPORT TransformedTriangle
101 friend class INTERP_TEST::TransformedTriangleTest;
102 friend class INTERP_TEST::TransformedTriangleIntersectTest;
104 * Enumerations representing the different geometric elements of the unit tetrahedron
105 * and the triangle. The end element, NO_* gives the number of elements in the enumeration
106 * and can be used as end element in loops.
109 /// Corners of tetrahedron
110 enum TetraCorner { O = 0, X, Y, Z, NO_TET_CORNER };
112 /// Edges of tetrahedron
113 enum TetraEdge { OX = 0, OY, OZ, XY, YZ, ZX, H01, H10, NO_TET_EDGE };
115 /// Facets (faces) of tetrahedron
116 enum TetraFacet { OYZ = 0, OZX, OXY, XYZ, NO_TET_FACET };
118 /// Corners of triangle
119 enum TriCorner { P = 0, Q, R, NO_TRI_CORNER };
121 /// Segments (edges) of triangle
122 enum TriSegment { PQ = 0, QR, RP, NO_TRI_SEGMENT };
124 /// Intersection polygons
125 enum IntersectionPolygon{ A = 0, B, NO_INTERSECTION_POLYGONS };
128 /// NB : order corresponds to TetraEdges (Grandy, table III)
129 enum DoubleProduct { C_YZ = 0, C_ZX, C_XY, C_ZH, C_XH, C_YH, C_01, C_10, NO_DP };
131 TransformedTriangle(double* p, double* q, double* r);
132 ~TransformedTriangle();
134 double calculateIntersectionVolume();
136 void dumpCoords() const;
138 // Queries of member values used by UnitTetraIntersectionBary
140 const double* getCorner(TriCorner corner) const { return _coords + 5*corner; }
142 const std::vector<double*>& getPolygonA() const { return _polygonA; }
144 double getVolume() const { return _volume; }
148 TransformedTriangle() { }
150 // ----------------------------------------------------------------------------------
151 // High-level methods called directly by calculateIntersectionVolume()
152 // ----------------------------------------------------------------------------------
153 void calculateIntersectionPolygons();
155 void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter);
157 void sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter);
159 double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter);
161 // ----------------------------------------------------------------------------------
162 // Detection of degenerate triangles
163 // ----------------------------------------------------------------------------------
165 bool isTriangleInPlaneOfFacet(const TetraFacet facet) const;
167 bool isTriangleParallelToFacet(const TetraFacet facet) const;
169 int isTriangleInclinedToFacet(const TetraFacet facet) const;
171 bool isTriangleBelowTetraeder() const;
173 // ----------------------------------------------------------------------------------
174 // Intersection test methods and intersection point calculations
175 // ----------------------------------------------------------------------------------
177 inline bool testSurfaceEdgeIntersection(const TetraEdge edge) const;
179 void calcIntersectionPtSurfaceEdge(const TetraEdge edge, double* pt) const;
181 inline bool testSegmentFacetIntersection(const TriSegment seg, const TetraFacet facet) const;
183 void calcIntersectionPtSegmentFacet(const TriSegment seg, const TetraFacet facet, double* pt) const;
185 bool testSegmentEdgeIntersection(const TriSegment seg, const TetraEdge edge) const;
187 void calcIntersectionPtSegmentEdge(const TriSegment seg, const TetraEdge edge, double* pt) const ;
189 bool testSegmentCornerIntersection(const TriSegment seg, const TetraCorner corner) const ;
191 inline bool testSurfaceRayIntersection(const TetraCorner corner) const;
193 bool testSegmentHalfstripIntersection(const TriSegment seg, const TetraEdge edg);
195 void calcIntersectionPtSegmentHalfstrip(const TriSegment seg, const TetraEdge edge, double* pt) const;
197 bool testSegmentRayIntersection(const TriSegment seg, const TetraCorner corner) const;
199 inline bool testCornerInTetrahedron(const TriCorner corner) const;
201 inline bool testCornerOnXYZFacet(const TriCorner corner) const;
203 inline bool testCornerAboveXYZFacet(const TriCorner corner) const;
205 // ----------------------------------------------------------------------------------
206 // Utility methods used in intersection tests
207 // ----------------------------------------------------------------------------------
209 bool testTriangleSurroundsEdge(const TetraEdge edge) const;
211 inline bool testEdgeIntersectsTriangle(const TetraEdge edge) const;
213 inline bool testFacetSurroundsSegment(const TriSegment seg, const TetraFacet facet) const;
215 inline bool testSegmentIntersectsFacet(const TriSegment seg, const TetraFacet facet) const;
217 bool testSegmentIntersectsHPlane(const TriSegment seg) const;
219 bool testSurfaceAboveCorner(const TetraCorner corner) const;
221 bool testTriangleSurroundsRay(const TetraCorner corner) const;
223 // ----------------------------------------------------------------------------------
224 // Double and triple product calculations
225 // ----------------------------------------------------------------------------------
229 bool areDoubleProductsConsistent(const TriSegment seg) const;
231 void preCalculateDoubleProducts(void);
233 inline void resetDoubleProducts(const TriSegment seg, const TetraCorner corner);
235 double calculateDistanceCornerSegment(const TetraCorner corner, const TriSegment seg) const;
237 void preCalculateTripleProducts(void);
239 double calculateAngleEdgeTriangle(const TetraEdge edge) const;
241 inline double calcStableC(const TriSegment seg, const DoubleProduct dp) const;
243 inline double calcStableT(const TetraCorner corner) const;
245 inline double calcUnstableC(const TriSegment seg, const DoubleProduct dp) const;
247 double calcTByDevelopingRow(const TetraCorner corner, const int row = 1, const bool project = false) const;
249 // ----------------------------------------------------------------------------------
251 // ----------------------------------------------------------------------------------
254 /// Array holding the coordinates of the triangle's three corners
256 /// [ 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 ]
259 /// Flag showing whether the double products have been calculated yet
260 bool _is_double_products_calculated;
262 /// Flag showing whether the triple products have been calculated yet
263 bool _is_triple_products_calculated;
265 /// Array containing the 24 double products.
266 /// order : c^PQ_YZ, ... ,cPQ_10, ... c^QR_YZ, ... c^RP_YZ
267 /// following order in enumeration DoubleProduct
268 double _doubleProducts[24];
270 /// Array containing the 4 triple products.
271 /// order : t_O, t_X, t_Y, t_Z
272 double _tripleProducts[4];
274 /// Vector holding the points of the intersection polygon A.
275 /// these points are allocated in calculateIntersectionPolygons() and liberated in the destructor
276 std::vector<double*> _polygonA;
278 /// Vector holding the points of the intersection polygon B.
279 /// These points are allocated in calculateIntersectionPolygons() and liberated in the destructor
280 std::vector<double*> _polygonB;
282 /// Array holding the coordinates of the barycenter of the polygon A
283 /// This point is calculated in calculatePolygonBarycenter
284 double _barycenterA[3];
286 /// Array holding the coordinates of the barycenter of the polygon B
287 /// This point is calculated in calculatePolygonBarycenter
288 //double _barycenterB[3];
290 /// Array of flags indicating which of the four triple products have been correctly calculated.
291 /// Used for asserts in debug mode
294 /// calculated volume for use of UnitTetraIntersectionBary
298 * Calls TransformedTriangle::testTriangleSurroundsEdge for edges OX to ZX and stores the result in
299 * member variable array_triangleSurroundsEdgeCache.
302 void preCalculateTriangleSurroundsEdge();
304 /// Array holding results of the test testTriangleSurroundsEdge() for all the edges.
305 /// These are calculated in preCalculateTriangleSurroundsEdge().
306 bool _triangleSurroundsEdgeCache[NO_TET_EDGE];
308 // ----------------------------------------------------------------------------------
310 // ----------------------------------------------------------------------------------
312 // offsets : 0 -> x, 1 -> y, 2 -> z, 3 -> h, 4 -> H
313 // corresponds to order of double products in DoubleProduct
314 // so that offset[C_*] gives the right coordinate
315 static const int DP_OFFSET_1[8];
316 static const int DP_OFFSET_2[8];
318 // the coordinates used in the expansion of triple products by a given row
319 // in constellation (corner, row-1)
320 // (0,1,2,3) <=> (x,y,z,h)
321 static const int COORDINATE_FOR_DETERMINANT_EXPANSION[12];
323 // contains the edge of the double product used when
324 // expanding the triple product determinant associated with each corner
326 static const DoubleProduct DP_FOR_DETERMINANT_EXPANSION[12];
328 // values used to decide how imprecise the double products
329 // should be to set them to 0.0
330 static const long double MACH_EPS; // machine epsilon
331 static const long double MULT_PREC_F; // precision of multiplications (Grandy : f)
332 static const long double THRESHOLD_F; // threshold for zeroing (Grandy : F/f)
334 static const double TRIPLE_PRODUCT_ANGLE_THRESHOLD;
336 // correspondance facet - double product
338 static const DoubleProduct DP_FOR_SEG_FACET_INTERSECTION[12];
340 // signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION
341 static const double SIGN_FOR_SEG_FACET_INTERSECTION[12];
343 // coordinates of corners of tetrahedron
344 static const double COORDS_TET_CORNER[12];
346 // indices to use in tables DP_FOR_SEG_FACET_INTERSECTION and SIGN_FOR_SEG_FACET_INTERSECTION
347 // for the calculation of the coordinates (x,y,z) of the intersection points
348 // for Segment-Facet and Segment-Edge intersections
349 static const int DP_INDEX[12];
351 // correspondance edge - corners
352 static const TetraCorner CORNERS_FOR_EDGE[12];
354 // correspondance edge - facets
355 // facets shared by each edge
356 static const TetraFacet FACET_FOR_EDGE[12];
358 // correspondance edge - corners
359 static const TetraEdge EDGES_FOR_CORNER[12];
361 // double products used in segment-halfstrip test
362 static const DoubleProduct DP_FOR_HALFSTRIP_INTERSECTION[12];
364 // double products used in segment - ray test
365 static const DoubleProduct DP_SEGMENT_RAY_INTERSECTION[21];
369 // include definitions of inline methods
371 #include "TransformedTriangleInline.hxx"