1 // Copyright (C) 2007-2013 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
20 #include "TransformedTriangle.hxx"
21 #include "VectorUtils.hxx"
22 #include "TetraAffineTransform.hxx"
32 namespace INTERP_KERNEL
36 * \brief Class representing a circular order of a set of points around their barycenter.
37 * It is used with the STL sort() algorithm to sort the point of the two polygons
40 class ProjectedCentralCircularSortOrder
44 /// Enumeration of different planes to project on when calculating order
45 enum CoordType { XY, XZ, YZ };
50 * @param barycenter double[3] containing the barycenter of the points to be compared
51 * @param type plane to project on when comparing. The comparison will not work if all the points are in a plane perpendicular
52 * to the plane being projected on
54 ProjectedCentralCircularSortOrder(const double* barycenter, const CoordType type)
55 : _aIdx((type == YZ) ? 1 : 0),
56 _bIdx((type == XY) ? 1 : 2),
57 _a(barycenter[_aIdx]),
63 * Comparison operator.
64 * Compares the relative position between two points in their ordering around the barycenter.
66 * @param pt1 a double[3] representing a point
67 * @param pt2 a double[3] representing a point
68 * @return true if the angle of the difference vector between pt1 and the barycenter is greater than that
69 * of the difference vector between pt2 and the barycenter.
71 bool operator()(const double* pt1, const double* pt2)
73 // calculate angles with the axis
74 const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b);
75 const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b);
81 /// index corresponding to first coordinate of plane on which points are projected
84 /// index corresponding to second coordinate of plane on which points are projected
87 /// value of first projected coordinate of the barycenter
90 /// value of second projected coordinate of the barycenter
94 // ----------------------------------------------------------------------------------
95 // TransformedTriangle PUBLIC
96 // ----------------------------------------------------------------------------------
101 * The coordinates are copied to the internal member variables
103 * @param p array of three doubles containing coordinates of P
104 * @param q array of three doubles containing coordinates of Q
105 * @param r array of three doubles containing coordinates of R
107 TransformedTriangle::TransformedTriangle(double* p, double* q, double* r)
108 : _is_double_products_calculated(false), _is_triple_products_calculated(false), _volume(0)
111 for(int i = 0 ; i < 3 ; ++i)
114 _coords[5*P + i] = p[i];
115 _coords[5*Q + i] = q[i];
116 _coords[5*R + i] = r[i];
121 _coords[5*P + 3] = 1 - p[0] - p[1] - p[2];
122 _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
123 _coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
126 _coords[5*P + 4] = 1 - p[0] - p[1];
127 _coords[5*Q + 4] = 1 - q[0] - q[1];
128 _coords[5*R + 4] = 1 - r[0] - r[1];
130 resetNearZeroCoordinates();
132 // initialise rest of data
133 preCalculateDoubleProducts();
135 preCalculateTriangleSurroundsEdge();
137 preCalculateTripleProducts();
144 * Deallocates the memory used to store the points of the polygons.
145 * This memory is allocated in calculateIntersectionAndProjectionPolygons().
147 TransformedTriangle::~TransformedTriangle()
149 // delete elements of polygons
150 for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
154 for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
161 * Calculates the volume of intersection between the triangle and the
164 * @return volume of intersection of this triangle with unit tetrahedron,
165 * as described in Grandy
168 double TransformedTriangle::calculateIntersectionVolume()
170 // check first that we are not below z - plane
171 if(isTriangleBelowTetraeder())
173 LOG(2, " --- Triangle is below tetraeder - V = 0.0");
177 // get the sign of the volume - equal to the sign of the z-component of the normal
178 // of the triangle, u_x * v_y - u_y * v_x, where u = q - p and v = r - p
179 // if it is zero, the triangle is perpendicular to the z - plane and so the volume is zero
180 // const double uv_xy[4] =
182 // _coords[5*Q] - _coords[5*P], _coords[5*Q + 1] - _coords[5*P + 1], // u_x, u_y
183 // _coords[5*R] - _coords[5*P], _coords[5*R + 1] - _coords[5*P + 1] // v_x, v_y
186 // double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
187 int sign = isTriangleInclinedToFacet( OXY );
191 LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0");
192 return _volume = 0.0;
197 //sign = sign > 0.0 ? 1.0 : -1.0;
199 LOG(2, "-- Calculating intersection polygons ... ");
200 calculateIntersectionAndProjectionPolygons();
202 double barycenter[3];
204 // calculate volume under A
206 if(_polygonA.size() > 2)
208 LOG(2, "---- Treating polygon A ... ");
209 calculatePolygonBarycenter(A, barycenter);
210 sortIntersectionPolygon(A, barycenter);
211 volA = calculateVolumeUnderPolygon(A, barycenter);
212 LOG(2, "Volume is " << sign * volA);
216 // if triangle is not in h = 0 plane, calculate volume under B
217 if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
219 LOG(2, "---- Treating polygon B ... ");
221 calculatePolygonBarycenter(B, barycenter);
222 sortIntersectionPolygon(B, barycenter);
223 volB = calculateVolumeUnderPolygon(B, barycenter);
224 LOG(2, "Volume is " << sign * volB);
227 LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
229 return _volume = sign * (volA + volB);
234 * Calculates the volume of intersection between the triangle and the
237 * @return volume of intersection of this triangle with unit tetrahedron,
238 * as described in Grandy
241 double TransformedTriangle::calculateIntersectionSurface(TetraAffineTransform* tat)
243 // check first that we are not below z - plane
244 if(isTriangleBelowTetraeder())
246 LOG(2, " --- Triangle is below tetraeder - V = 0.0");
250 LOG(2, "-- Calculating intersection polygon ... ");
251 calculateIntersectionPolygon();
254 if(_polygonA.size() > 2) {
255 double barycenter[3];
256 calculatePolygonBarycenter(A, barycenter);
257 sortIntersectionPolygon(A, barycenter);
258 const std::size_t nbPoints = _polygonA.size();
259 for(std::size_t i = 0 ; i < nbPoints ; ++i)
260 tat->reverseApply(_polygonA[i], _polygonA[i]);
261 _volume = calculateSurfacePolygon();
267 // ----------------------------------------------------------------------------------
268 // TransformedTriangle PRIVATE
269 // ----------------------------------------------------------------------------------
272 * Calculates the intersection polygons A and B, performing the intersection tests
273 * and storing the corresponding points in the vectors _polygonA and _polygonB.
275 * @post _polygonA contains the intersection polygon A and _polygonB contains the
276 * intersection polygon B.
279 void TransformedTriangle::calculateIntersectionAndProjectionPolygons()
281 assert(_polygonA.size() == 0);
282 assert(_polygonB.size() == 0);
283 // avoid reallocations in push_back() by pre-allocating enough memory
284 // we should never have more than 20 points
285 _polygonA.reserve(20);
286 _polygonB.reserve(20);
287 // -- surface intersections
289 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
291 if(testSurfaceEdgeIntersection(edge))
293 double* ptA = new double[3];
294 calcIntersectionPtSurfaceEdge(edge, ptA);
295 _polygonA.push_back(ptA);
296 LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
299 double* ptB = new double[3];
300 copyVector3(ptA, ptB);
301 _polygonB.push_back(ptB);
302 LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B ");
309 for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
311 if(testSurfaceRayIntersection(corner))
313 double* ptB = new double[3];
314 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
315 _polygonB.push_back(ptB);
316 LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
320 // -- segment intersections
321 for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
326 // check beforehand which double-products are zero
327 for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
329 isZero[dp] = (calcStableC(seg, dp) == 0.0);
333 for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
335 // is this test worth it?
337 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
338 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
339 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
341 if(doTest && testSegmentFacetIntersection(seg, facet))
343 double* ptA = new double[3];
344 calcIntersectionPtSegmentFacet(seg, facet, ptA);
345 _polygonA.push_back(ptA);
346 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
349 double* ptB = new double[3];
350 copyVector3(ptA, ptB);
351 _polygonB.push_back(ptB);
352 LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
359 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
361 const DoubleProduct edge_dp = DoubleProduct(edge);
363 if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
365 double* ptA = new double[3];
366 calcIntersectionPtSegmentEdge(seg, edge, ptA);
367 _polygonA.push_back(ptA);
368 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
371 double* ptB = new double[3];
372 copyVector3(ptA, ptB);
373 _polygonB.push_back(ptB);
379 for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
382 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
383 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
384 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
386 if(doTest && testSegmentCornerIntersection(seg, corner))
388 double* ptA = new double[3];
389 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
390 _polygonA.push_back(ptA);
391 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
394 double* ptB = new double[3];
395 _polygonB.push_back(ptB);
396 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
397 LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
403 for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
405 if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
407 double* ptB = new double[3];
408 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
409 _polygonB.push_back(ptB);
410 LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
414 // segment - halfstrip
415 for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
419 const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
421 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
422 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
425 if(doTest && testSegmentHalfstripIntersection(seg, edge))
427 if(testSegmentHalfstripIntersection(seg, edge))
429 double* ptB = new double[3];
430 calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
431 _polygonB.push_back(ptB);
432 LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
438 for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
440 // { XYZ - inclusion only possible if in Tetrahedron?
442 if(testCornerInTetrahedron(corner))
444 double* ptA = new double[3];
445 copyVector3(&_coords[5*corner], ptA);
446 _polygonA.push_back(ptA);
447 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
451 if(testCornerOnXYZFacet(corner))
453 double* ptB = new double[3];
454 copyVector3(&_coords[5*corner], ptB);
455 _polygonB.push_back(ptB);
456 LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
459 // projection on XYZ - facet
460 if(testCornerAboveXYZFacet(corner))
462 double* ptB = new double[3];
463 copyVector3(&_coords[5*corner], ptB);
464 ptB[2] = 1 - ptB[0] - ptB[1];
465 assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
466 _polygonB.push_back(ptB);
467 LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
475 * Calculates the intersection polygon A, performing the intersection tests
476 * and storing the corresponding point in the vector _polygonA.
478 * @post _polygonA contains the intersection polygon A.
481 void TransformedTriangle::calculateIntersectionPolygon()
483 assert(_polygonA.size() == 0);
484 // avoid reallocations in push_back() by pre-allocating enough memory
485 // we should never have more than 20 points
486 _polygonA.reserve(20);
487 // -- surface intersections
489 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
491 if(testSurfaceEdgeIntersection(edge))
493 double* ptA = new double[3];
494 calcIntersectionPtSurfaceEdge(edge, ptA);
495 _polygonA.push_back(ptA);
496 LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
500 // -- segment intersections
501 for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
506 // check beforehand which double-products are zero
507 for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
509 isZero[dp] = (calcStableC(seg, dp) == 0.0);
513 for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
515 // is this test worth it?
517 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
518 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
519 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
521 if(doTest && testSegmentFacetIntersection(seg, facet))
523 double* ptA = new double[3];
524 calcIntersectionPtSegmentFacet(seg, facet, ptA);
525 _polygonA.push_back(ptA);
526 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
531 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
533 const DoubleProduct edge_dp = DoubleProduct(edge);
535 if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
537 double* ptA = new double[3];
538 calcIntersectionPtSegmentEdge(seg, edge, ptA);
539 _polygonA.push_back(ptA);
540 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
545 for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
548 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
549 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
550 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
552 if(doTest && testSegmentCornerIntersection(seg, corner))
554 double* ptA = new double[3];
555 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
556 _polygonA.push_back(ptA);
557 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
564 for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
566 // { XYZ - inclusion only possible if in Tetrahedron?
568 if(testCornerInTetrahedron(corner))
570 double* ptA = new double[3];
571 copyVector3(&_coords[5*corner], ptA);
572 _polygonA.push_back(ptA);
573 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
582 * Returns the surface of polygon A.
584 * @return the surface of polygon A.
586 double TransformedTriangle::calculateSurfacePolygon()
588 const std::size_t nbPoints = _polygonA.size();
590 double sum[3] = {0., 0., 0.};
592 for(std::size_t i = 0 ; i < nbPoints ; ++i)
594 const double *const ptCurr = _polygonA[i]; // pt "i"
595 const double *const ptNext = _polygonA[(i + 1) % nbPoints]; // pt "i+1" (pt nbPoints == pt 0)
597 cross(ptCurr, ptNext, pdt);
601 const double surface = norm(sum) * 0.5;
602 LOG(2,"Surface is " << surface);
607 * Calculates the barycenters of the given intersection polygon.
609 * @pre the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons()
611 * @param poly one of the two intersection polygons
612 * @param barycenter array of three doubles where barycenter is stored
615 void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
617 LOG(3,"--- Calculating polygon barycenter");
619 // get the polygon points
620 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
622 // calculate barycenter
623 const std::size_t m = polygon.size();
625 for(int j = 0 ; j < 3 ; ++j)
632 for(std::size_t i = 0 ; i < m ; ++i)
634 const double* pt = polygon[i];
635 for(int j = 0 ; j < 3 ; ++j)
637 barycenter[j] += pt[j] / double(m);
641 LOG(3,"Barycenter is " << vToStr(barycenter));
645 * Sorts the given intersection polygon in circular order around its barycenter.
646 * @pre the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons()
647 * @post the vertices in _polygonA and _polygonB are sorted in circular order around their
648 * respective barycenters
650 * @param poly one of the two intersection polygons
651 * @param barycenter array of three doubles with the coordinates of the barycenter
654 void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
656 LOG(3,"--- Sorting polygon ...");
658 using INTERP_KERNEL::ProjectedCentralCircularSortOrder;
659 typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
660 typedef SortOrder::CoordType CoordType;
662 // get the polygon points
663 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
665 if(polygon.size() == 0)
668 // determine type of sorting
669 CoordType type = SortOrder::XY;
670 if(poly == A && !isTriangleInclinedToFacet( OXY )) // B is on h = 0 plane -> ok
672 // NB : the following test is never true if we have eliminated the
673 // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
674 // We keep the test here anyway, to avoid interdependency.
676 // is triangle inclined to x == 0 ?
677 if(isTriangleInclinedToFacet(OZX))
679 type = SortOrder::XZ;
681 else //if(isTriangleParallelToFacet(OYZ))
683 type = SortOrder::YZ;
687 // create order object
688 SortOrder order(barycenter, type);
690 // sort vector with this object
691 // NB : do not change place of first object, with respect to which the order
693 sort((polygon.begin()), polygon.end(), order);
695 LOG(3,"Sorted polygon is ");
696 for(size_t i = 0 ; i < polygon.size() ; ++i)
698 LOG(3,vToStr(polygon[i]));
704 * Calculates the volume between the given polygon and the z = 0 plane.
706 * @pre the intersection polygones have been calculated with calculateIntersectionAndProjectionPolygons(),
707 * and they have been sorted in circular order with sortIntersectionPolygons(void)
709 * @param poly one of the two intersection polygons
710 * @param barycenter array of three doubles with the coordinates of the barycenter
711 * @return the volume between the polygon and the z = 0 plane
714 double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter)
716 LOG(2,"--- Calculating volume under polygon");
718 // get the polygon points
719 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
722 const std::size_t m = polygon.size();
724 for(std::size_t i = 0 ; i < m ; ++i)
726 const double* ptCurr = polygon[i]; // pt "i"
727 const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
729 const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
730 const double factor2 =
731 ptCurr[0]*(ptNext[1] - barycenter[1])
732 + ptNext[0]*(barycenter[1] - ptCurr[1])
733 + barycenter[0]*(ptCurr[1] - ptNext[1]);
734 vol += (factor1 * factor2) / 6.0;
737 LOG(2,"Abs. Volume is " << vol);
742 ////////////////////////////////////////////////////////////////////////////////////
743 // Detection of (very) degenerate cases /////////////
744 ////////////////////////////////////////////////////////////////////////////////////
747 * Checks if the triangle lies in the plane of a given facet
749 * @param facet one of the facets of the tetrahedron
750 * @return true if PQR lies in the plane of the facet, false if not
752 bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
755 // coordinate to check
756 const int coord = static_cast<int>(facet);
758 for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
760 if(_coords[5*c + coord] != 0.0)
770 * Checks if the triangle is parallel to the given facet
772 * @param facet one of the facets of the unit tetrahedron
773 * @return true if triangle is parallel to facet, false if not
775 bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const
777 // coordinate to check
778 const int coord = static_cast<int>(facet);
779 return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]);
783 * Checks if the triangle is not perpedicular to the given facet
785 * @param facet one of the facets of the unit tetrahedron
786 * @return zero if the triangle is perpendicular to the facet,
787 * else 1 or -1 depending on the sign of cross product of facet edges
789 int TransformedTriangle::isTriangleInclinedToFacet(const TetraFacet facet) const
791 // coordinate to check
792 const int coord = static_cast<int>(facet);
793 const int ind1 = ( coord+1 ) % 3, ind2 = ( coord+2 ) % 3;
794 const double uv_xy[4] =
797 _coords[5*Q+ind1] - _coords[5*P+ind1], _coords[5*Q+ind2] - _coords[5*P+ind2],
799 _coords[5*R+ind1] - _coords[5*P+ind1], _coords[5*R+ind2] - _coords[5*P+ind2]
802 double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
803 if(epsilonEqual(sign, 0.))
807 return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0;
811 * Determines whether the triangle is below the z-plane.
813 * @return true if the z-coordinate of the three corners of the triangle are all less than 0, false otherwise.
815 bool TransformedTriangle::isTriangleBelowTetraeder() const
817 for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
819 // check z-coords for all points
820 if(_coords[5*c + 2] >= 0.0)
829 * Prints the coordinates of the triangle to std::cout
832 void TransformedTriangle::dumpCoords() const
834 std::cout << "Coords : ";
835 for(int i = 0 ; i < 3; ++i)
837 std::cout << vToStr(&_coords[5*i]) << ",";
839 std::cout << std::endl;