1 // Copyright (C) 2007-2024 CEA, EDF
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 #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 // A ***much*** faster alternative to atan2 to get a pseudo-angle suitable for sorting:
74 // https://stackoverflow.com/questions/16542042/fastest-way-to-sort-vectors-by-angle-without-actually-computing-that-angle
75 const double dy1 = pt1[_aIdx] - _a, dx1 = pt1[_bIdx] - _b,
76 dy2 = pt2[_aIdx] - _a, dx2 = pt2[_bIdx] - _b;
77 const double deno1 = std::fabs(dx1)+fabs(dy1),
78 deno2 = std::fabs(dx2)+fabs(dy2);
80 const double ang1 = deno1 == 0.0 ? 0.0 : std::copysign(1. - dx1/deno1, dy1),
81 ang2 = deno2 == 0.0 ? 0.0 : std::copysign(1. - dx2/deno2, dy2);
87 /// index corresponding to first coordinate of plane on which points are projected
90 /// index corresponding to second coordinate of plane on which points are projected
93 /// value of first projected coordinate of the barycenter
96 /// value of second projected coordinate of the barycenter
100 // ----------------------------------------------------------------------------------
101 // TransformedTriangle PUBLIC
102 // ----------------------------------------------------------------------------------
107 * The coordinates are copied to the internal member variables
109 * @param p array of three doubles containing coordinates of P
110 * @param q array of three doubles containing coordinates of Q
111 * @param r array of three doubles containing coordinates of R
113 TransformedTriangle::TransformedTriangle(double* p, double* q, double* r)
114 : _is_double_products_calculated(false), _is_triple_products_calculated(false), _volume(0)
117 for(int i = 0 ; i < 3 ; ++i)
120 _coords[5*P + i] = p[i];
121 _coords[5*Q + i] = q[i];
122 _coords[5*R + i] = r[i];
127 _coords[5*P + 3] = 1 - p[0] - p[1] - p[2];
128 _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
129 _coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
131 // Handle degenerated case where one of the seg of PQR is (almost) inside XYZ plane,
132 // and hence by extension when the whole PQR triangle is in the XYZ plane
133 handleDegenerateCases();
136 _coords[5*P + 4] = 1 - p[0] - p[1];
137 _coords[5*Q + 4] = 1 - q[0] - q[1];
138 _coords[5*R + 4] = 1 - r[0] - r[1];
140 // initialise rest of data
141 preCalculateDoubleProducts();
143 preCalculateTriangleSurroundsEdge();
145 preCalculateTripleProducts();
152 * Deallocates the memory used to store the points of the polygons.
153 * This memory is allocated in calculateIntersectionAndProjectionPolygons().
155 TransformedTriangle::~TransformedTriangle()
157 // delete elements of polygons
158 for(auto& it: _polygonA)
160 for(auto& it: _polygonB)
165 * Calculates the volume of intersection between the triangle and the
168 * @return volume of intersection of this triangle with unit tetrahedron,
169 * as described in Grandy
172 double TransformedTriangle::calculateIntersectionVolume()
174 // check first that we are not below z - plane
175 if(isTriangleBelowTetraeder())
177 LOG(2, " --- Triangle is below tetraeder - V = 0.0");
181 // get the sign of the volume - equal to the sign of the z-component of the normal
182 // of the triangle, u_x * v_y - u_y * v_x, where u = q - p and v = r - p
183 // if it is zero, the triangle is perpendicular to the z - plane and so the volume is zero
184 // const double uv_xy[4] =
186 // _coords[5*Q] - _coords[5*P], _coords[5*Q + 1] - _coords[5*P + 1], // u_x, u_y
187 // _coords[5*R] - _coords[5*P], _coords[5*R + 1] - _coords[5*P + 1] // v_x, v_y
190 // double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
191 int sign = isTriangleInclinedToFacet( OXY );
195 LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0");
196 return _volume = 0.0;
201 //sign = sign > 0.0 ? 1.0 : -1.0;
203 LOG(2, "-- Calculating intersection polygons ... ");
204 calculateIntersectionAndProjectionPolygons();
206 double barycenter[3];
208 // calculate volume under A
210 if(_polygonA.size() > 2)
212 LOG(2, "---- Treating polygon A ... ");
214 LOG(3, " --- Final points in polygon A");
215 for(const auto& pt: _polygonA)
218 calculatePolygonBarycenter(A, barycenter);
219 sortIntersectionPolygon(A, barycenter);
220 volA = calculateVolumeUnderPolygon(A, barycenter);
221 LOG(2, "Volume is " << sign * volA);
225 // if triangle is not in h = 0 plane, calculate volume under B
226 if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ)) // second condition avoids double counting in case triangle fully included in h=0 facet
228 LOG(2, "---- Treating polygon B ... ");
230 LOG(3, " --- Final points in polygon B");
231 for(const auto& pt: _polygonB)
234 calculatePolygonBarycenter(B, barycenter);
235 sortIntersectionPolygon(B, barycenter);
236 volB = calculateVolumeUnderPolygon(B, barycenter);
237 LOG(2, "Volume is " << sign * volB);
241 LOG(2, "############ Triangle :")
243 LOG(2, "vol A = " << volA);
244 LOG(2, "vol B = " << volB);
245 LOG(2, "TOTAL = " << sign*(volA+volB));
248 return _volume = sign * (volA + volB);
253 * Calculates the volume of intersection between the triangle and the
256 * @return volume of intersection of this triangle with unit tetrahedron,
257 * as described in Grandy
260 double TransformedTriangle::calculateIntersectionSurface(TetraAffineTransform* tat)
262 // check first that we are not below z - plane
263 if(isTriangleBelowTetraeder())
265 LOG(2, " --- Triangle is below tetraeder - V = 0.0");
269 LOG(2, "-- Calculating intersection polygon ... ");
270 calculateIntersectionPolygon();
273 if(_polygonA.size() > 2) {
274 double barycenter[3];
275 calculatePolygonBarycenter(A, barycenter);
276 sortIntersectionPolygon(A, barycenter);
277 const std::size_t nbPoints = _polygonA.size();
278 for(std::size_t i = 0 ; i < nbPoints ; ++i)
279 tat->reverseApply(_polygonA[i], _polygonA[i]);
280 _volume = calculateSurfacePolygon();
286 // ----------------------------------------------------------------------------------
287 // TransformedTriangle PROTECTED
288 // ----------------------------------------------------------------------------------
291 * Calculates the intersection polygons A and B, performing the intersection tests
292 * and storing the corresponding points in the vectors _polygonA and _polygonB.
294 * @post _polygonA contains the intersection polygon A and _polygonB contains the
295 * intersection polygon B.
298 void TransformedTriangle::calculateIntersectionAndProjectionPolygons()
301 std::cout << " @@@@@@@@ COORDS @@@@@@ " << std::endl;
305 assert(_polygonA.size() == 0);
306 assert(_polygonB.size() == 0);
307 // avoid reallocations in push_back() by pre-allocating enough memory
308 // we should never have more than 20 points
309 _polygonA.reserve(20);
310 _polygonB.reserve(20);
311 // -- surface intersections
313 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
315 if(testSurfaceEdgeIntersection(edge))
317 double* ptA = new double[3];
318 calcIntersectionPtSurfaceEdge(edge, ptA);
319 _polygonA.push_back(ptA);
320 LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptA) << " added to A ");
323 double* ptB = new double[3];
324 copyVector3(ptA, ptB);
325 _polygonB.push_back(ptB);
326 LOG(3,"Surface-edge (edge " << strTE(edge) << "): " << vToStr(ptB) << " added to B ");
333 for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
335 if(testSurfaceRayIntersection(corner))
337 double* ptB = new double[3];
338 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
339 _polygonB.push_back(ptB);
340 LOG(3,"Surface-ray (corner " << strTC(corner) << "): " << vToStr(ptB) << " added to B");
344 // -- segment intersections
345 for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
350 // check beforehand which double-products are zero.
351 for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
352 isZero[dp] = (calcStableC(seg, dp) == 0.0);
355 for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
357 // is this test worth it?
359 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
360 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
361 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
363 if(doTest && testSegmentFacetIntersection(seg, facet))
365 double* ptA = new double[3];
366 calcIntersectionPtSegmentFacet(seg, facet, ptA);
367 _polygonA.push_back(ptA);
368 LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
371 double* ptB = new double[3];
372 copyVector3(ptA, ptB);
373 _polygonB.push_back(ptB);
374 LOG(3,"Segment-facet (facet " << strTF(facet) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
381 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
383 const DoubleProduct edge_dp = DoubleProduct(edge);
385 if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
387 double* ptA = new double[3];
388 calcIntersectionPtSegmentEdge(seg, edge, ptA);
389 _polygonA.push_back(ptA);
390 LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
393 double* ptB = new double[3];
394 copyVector3(ptA, ptB);
395 _polygonB.push_back(ptB);
396 LOG(3,"Segment-edge (edge " << strTE(edge) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to B");
402 for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
405 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
406 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
407 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
409 if(doTest && testSegmentCornerIntersection(seg, corner))
411 double* ptA = new double[3];
412 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
413 _polygonA.push_back(ptA);
414 LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptA) << " added to A");
417 double* ptB = new double[3];
418 _polygonB.push_back(ptB);
419 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
420 LOG(3,"Segment-corner (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
426 for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
428 if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
430 double* ptB = new double[3];
431 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
432 _polygonB.push_back(ptB);
433 LOG(3,"Segment-ray (corner " << strTC(corner) << ", seg " << strTriS(seg) << "): " << vToStr(ptB) << " added to B");
437 // segment - halfstrip
438 for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
442 const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
444 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
445 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
448 if(doTest && testSegmentHalfstripIntersection(seg, edge))
450 if(testSegmentHalfstripIntersection(seg, edge))
452 double* ptB = new double[3];
453 calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
454 _polygonB.push_back(ptB);
455 LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
461 for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
463 // { XYZ - inclusion only possible if in Tetrahedron?
465 if(testCornerInTetrahedron(corner))
467 double* ptA = new double[3];
468 copyVector3(&_coords[5*corner], ptA);
469 _polygonA.push_back(ptA);
470 LOG(3,"Inclusion tetrahedron (corner " << strTriC(corner) << "): " << vToStr(ptA) << " added to A");
474 if(testCornerOnXYZFacet(corner))
476 double* ptB = new double[3];
477 copyVector3(&_coords[5*corner], ptB);
478 _polygonB.push_back(ptB);
479 LOG(3,"Inclusion XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B");
482 // projection on XYZ - facet
483 if(testCornerAboveXYZFacet(corner))
485 double* ptB = new double[3];
486 copyVector3(&_coords[5*corner], ptB);
487 ptB[2] = 1 - ptB[0] - ptB[1]; // lower z to project on XYZ
488 assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
489 _polygonB.push_back(ptB);
490 LOG(3,"Projection XYZ-plane (corner " << strTriC(corner) << "): " << vToStr(ptB) << " added to B");
498 * Calculates the intersection polygon A, performing the intersection tests
499 * and storing the corresponding point in the vector _polygonA.
501 * @post _polygonA contains the intersection polygon A.
504 void TransformedTriangle::calculateIntersectionPolygon()
506 assert(_polygonA.size() == 0);
507 // avoid reallocations in push_back() by pre-allocating enough memory
508 // we should never have more than 20 points
509 _polygonA.reserve(20);
510 // -- surface intersections
512 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
514 if(testSurfaceEdgeIntersection(edge))
516 double* ptA = new double[3];
517 calcIntersectionPtSurfaceEdge(edge, ptA);
518 _polygonA.push_back(ptA);
519 LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
523 // -- segment intersections
524 for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
529 // check beforehand which double-products are zero
530 // Test for "== 0.0" here is OK since doubleProduct has been fixed for rounding to zero already.
531 for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
532 isZero[dp] = (calcStableC(seg, dp) == 0.0);
535 for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
537 // is this test worth it?
539 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
540 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
541 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
543 if(doTest && testSegmentFacetIntersection(seg, facet))
545 double* ptA = new double[3];
546 calcIntersectionPtSegmentFacet(seg, facet, ptA);
547 _polygonA.push_back(ptA);
548 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
553 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
555 const DoubleProduct edge_dp = DoubleProduct(edge);
557 if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
559 double* ptA = new double[3];
560 calcIntersectionPtSegmentEdge(seg, edge, ptA);
561 _polygonA.push_back(ptA);
562 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
567 for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
570 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
571 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
572 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
574 if(doTest && testSegmentCornerIntersection(seg, corner))
576 double* ptA = new double[3];
577 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
578 _polygonA.push_back(ptA);
579 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
586 for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
588 // { XYZ - inclusion only possible if in Tetrahedron?
590 if(testCornerInTetrahedron(corner))
592 double* ptA = new double[3];
593 copyVector3(&_coords[5*corner], ptA);
594 _polygonA.push_back(ptA);
595 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
604 * Returns the surface of polygon A.
606 * @return the surface of polygon A.
608 double TransformedTriangle::calculateSurfacePolygon()
610 const std::size_t nbPoints = _polygonA.size();
612 double sum[3] = {0., 0., 0.};
614 for(std::size_t i = 0 ; i < nbPoints ; ++i)
616 const double *const ptCurr = _polygonA[i]; // pt "i"
617 const double *const ptNext = _polygonA[(i + 1) % nbPoints]; // pt "i+1" (pt nbPoints == pt 0)
619 cross(ptCurr, ptNext, pdt);
623 const double surface = norm(sum) * 0.5;
624 LOG(2,"Surface is " << surface);
629 * Calculates the barycenters of the given intersection polygon.
631 * @pre the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons()
633 * @param poly one of the two intersection polygons
634 * @param barycenter array of three doubles where barycenter is stored
637 void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
639 LOG(3,"--- Calculating polygon barycenter");
641 // get the polygon points
642 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
644 // calculate barycenter
645 const std::size_t m = polygon.size();
647 for(int j = 0 ; j < 3 ; ++j)
650 for(std::size_t i = 0 ; i < m ; ++i)
652 const double* pt = polygon[i];
653 for(int j = 0 ; j < 3 ; ++j)
654 barycenter[j] += pt[j];
659 for(int j = 0 ; j < 3 ; ++j)
660 barycenter[j] = barycenter[j] / double(m);
662 LOG(3,"Barycenter is " << vToStr(barycenter));
666 * Sorts the given intersection polygon in circular order around its barycenter.
667 * @pre the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons()
668 * @post the vertices in _polygonA and _polygonB are sorted in circular order around their
669 * respective barycenters
671 * @param poly one of the two intersection polygons
672 * @param barycenter array of three doubles with the coordinates of the barycenter
675 void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
677 LOG(3,"--- Sorting polygon ...");
679 using INTERP_KERNEL::ProjectedCentralCircularSortOrder;
680 typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
681 typedef SortOrder::CoordType CoordType;
683 // get the polygon points
684 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
686 if(polygon.size() == 0)
689 // determine type of sorting
690 CoordType type = SortOrder::XY;
691 if(poly == A && !isTriangleInclinedToFacet( OXY )) // B is on h = 0 plane -> ok
693 // NB : the following test is never true if we have eliminated the
694 // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
695 // We keep the test here anyway, to avoid interdependency.
697 // is triangle inclined to x == 0 ?
698 type = isTriangleInclinedToFacet(OZX) ? SortOrder::XZ : SortOrder::YZ;
701 // create order object
702 SortOrder order(barycenter, type);
704 // sort vector with this object
705 // NB : do not change place of first object, with respect to which the order
707 sort((polygon.begin()), polygon.end(), order);
709 LOG(3,"Sorted polygon is ");
710 for(size_t i = 0 ; i < polygon.size() ; ++i)
712 LOG(3,vToStr(polygon[i]));
718 * Calculates the volume between the given polygon and the z = 0 plane.
720 * @pre the intersection polygones have been calculated with calculateIntersectionAndProjectionPolygons(),
721 * and they have been sorted in circular order with sortIntersectionPolygons(void)
723 * @param poly one of the two intersection polygons
724 * @param barycenter array of three doubles with the coordinates of the barycenter
725 * @return the volume between the polygon and the z = 0 plane
728 double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter)
730 LOG(2,"--- Calculating volume under polygon");
732 // get the polygon points
733 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
736 const std::size_t m = polygon.size();
738 for(std::size_t i = 0 ; i < m ; ++i)
740 const double* ptCurr = polygon[i]; // pt "i"
741 const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
743 const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
744 const double factor2 =
745 ptCurr[0]*(ptNext[1] - barycenter[1])
746 + ptNext[0]*(barycenter[1] - ptCurr[1])
747 + barycenter[0]*(ptCurr[1] - ptNext[1]);
748 vol += (factor1 * factor2) / 6.0;
751 LOG(2,"Abs. Volume is " << vol);
756 ////////////////////////////////////////////////////////////////////////////////////
757 // Detection of (very) degenerate cases /////////////
758 ////////////////////////////////////////////////////////////////////////////////////
761 * Checks if the triangle lies in the plane of a given facet
763 * @param facet one of the facets of the tetrahedron
764 * @return true if PQR lies in the plane of the facet, false if not
766 bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
768 // coordinate to check
769 const int coord = static_cast<int>(facet);
771 for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
772 if(_coords[5*c + coord] != 0.0)
779 * Checks if the triangle is parallel to the given facet
781 * @param facet one of the facets of the unit tetrahedron
782 * @return true if triangle is parallel to facet, false if not
784 bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const
786 // coordinate to check
787 const int coord = static_cast<int>(facet);
788 return (epsilonEqual(_coords[5*P + coord], _coords[5*Q + coord])) && (epsilonEqual(_coords[5*P + coord], _coords[5*R + coord]));
792 * Checks if the triangle is not perpedicular to the given facet
794 * @param facet one of the facets of the unit tetrahedron
795 * @return zero if the triangle is perpendicular to the facet,
796 * else 1 or -1 depending on the sign of cross product of facet edges
798 int TransformedTriangle::isTriangleInclinedToFacet(const TetraFacet facet) const
800 // coordinate to check
801 const int coord = static_cast<int>(facet);
802 const int ind1 = ( coord+1 ) % 3, ind2 = ( coord+2 ) % 3;
803 const double uv_xy[4] =
806 _coords[5*Q+ind1] - _coords[5*P+ind1], _coords[5*Q+ind2] - _coords[5*P+ind2],
808 _coords[5*R+ind1] - _coords[5*P+ind1], _coords[5*R+ind2] - _coords[5*P+ind2]
811 double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
812 if(epsilonEqual(sign, 0.))
814 return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0;
818 * Determines whether the triangle is below the z-plane.
820 * @return true if the z-coordinate of the three corners of the triangle are all less than 0, false otherwise.
822 bool TransformedTriangle::isTriangleBelowTetraeder() const
824 for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
825 // check z-coords for all points
826 if(_coords[5*c + 2] >= 0.0)
833 * Prints the coordinates of the triangle to std::cout
836 void TransformedTriangle::dumpCoords() const
838 std::cout << "Coords : ";
839 for(int i = 0 ; i < 3; ++i)
840 std::cout << vToStr(&_coords[5*i]) << ",";
842 std::cout << std::endl;