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 #include "TransformedTriangle.hxx"
20 #include "VectorUtils.hxx"
31 namespace INTERP_KERNEL
35 * \brief Class representing a circular order of a set of points around their barycenter.
36 * It is used with the STL sort() algorithm to sort the point of the two polygons
39 class ProjectedCentralCircularSortOrder
43 /// Enumeration of different planes to project on when calculating order
44 enum CoordType { XY, XZ, YZ };
49 * @param barycenter double[3] containing the barycenter of the points to be compared
50 * @param type plane to project on when comparing. The comparison will not work if all the points are in a plane perpendicular
51 * to the plane being projected on
53 ProjectedCentralCircularSortOrder(const double* barycenter, const CoordType type)
54 : _aIdx((type == YZ) ? 1 : 0),
55 _bIdx((type == XY) ? 1 : 2),
56 _a(barycenter[_aIdx]),
62 * Comparison operator.
63 * Compares the relative position between two points in their ordering around the barycenter.
65 * @param pt1 a double[3] representing a point
66 * @param pt2 a double[3] representing a point
67 * @return true if the angle of the difference vector between pt1 and the barycenter is greater than that
68 * of the difference vector between pt2 and the barycenter.
70 bool operator()(const double* pt1, const double* pt2)
72 // calculate angles with the axis
73 const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b);
74 const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b);
80 /// index corresponding to first coordinate of plane on which points are projected
83 /// index corresponding to second coordinate of plane on which points are projected
86 /// value of first projected coordinate of the barycenter
89 /// value of second projected coordinate of the barycenter
93 // ----------------------------------------------------------------------------------
94 // TransformedTriangle PUBLIC
95 // ----------------------------------------------------------------------------------
100 * The coordinates are copied to the internal member variables
102 * @param p array of three doubles containing coordinates of P
103 * @param q array of three doubles containing coordinates of Q
104 * @param r array of three doubles containing coordinates of R
106 TransformedTriangle::TransformedTriangle(double* p, double* q, double* r)
107 : _is_double_products_calculated(false), _is_triple_products_calculated(false), _volume(0)
110 for(int i = 0 ; i < 3 ; ++i)
113 _coords[5*P + i] = p[i];
114 _coords[5*Q + i] = q[i];
115 _coords[5*R + i] = r[i];
120 _coords[5*P + 3] = 1 - p[0] - p[1] - p[2];
121 _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
122 _coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
125 _coords[5*P + 4] = 1 - p[0] - p[1];
126 _coords[5*Q + 4] = 1 - q[0] - q[1];
127 _coords[5*R + 4] = 1 - r[0] - r[1];
129 // initialise rest of data
130 preCalculateDoubleProducts();
132 preCalculateTriangleSurroundsEdge();
134 preCalculateTripleProducts();
141 * Deallocates the memory used to store the points of the polygons.
142 * This memory is allocated in calculateIntersectionPolygons().
144 TransformedTriangle::~TransformedTriangle()
146 // delete elements of polygons
147 for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
151 for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
158 * Calculates the volume of intersection between the triangle and the
161 * @return volume of intersection of this triangle with unit tetrahedron,
162 * as described in Grandy
165 double TransformedTriangle::calculateIntersectionVolume()
167 // check first that we are not below z - plane
168 if(isTriangleBelowTetraeder())
170 LOG(2, " --- Triangle is below tetraeder - V = 0.0");
174 // get the sign of the volume - equal to the sign of the z-component of the normal
175 // of the triangle, u_x * v_y - u_y * v_x, where u = q - p and v = r - p
176 // if it is zero, the triangle is perpendicular to the z - plane and so the volume is zero
177 // const double uv_xy[4] =
179 // _coords[5*Q] - _coords[5*P], _coords[5*Q + 1] - _coords[5*P + 1], // u_x, u_y
180 // _coords[5*R] - _coords[5*P], _coords[5*R + 1] - _coords[5*P + 1] // v_x, v_y
183 // double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
184 int sign = isTriangleInclinedToFacet( OXY );
188 LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0");
189 return _volume = 0.0;
194 //sign = sign > 0.0 ? 1.0 : -1.0;
196 LOG(2, "-- Calculating intersection polygons ... ");
197 calculateIntersectionPolygons();
199 double barycenter[3];
201 // calculate volume under A
203 if(_polygonA.size() > 2)
205 LOG(2, "---- Treating polygon A ... ");
206 calculatePolygonBarycenter(A, barycenter);
207 sortIntersectionPolygon(A, barycenter);
208 volA = calculateVolumeUnderPolygon(A, barycenter);
209 LOG(2, "Volume is " << sign * volA);
213 // if triangle is not in h = 0 plane, calculate volume under B
214 if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
216 LOG(2, "---- Treating polygon B ... ");
218 calculatePolygonBarycenter(B, barycenter);
219 sortIntersectionPolygon(B, barycenter);
220 volB = calculateVolumeUnderPolygon(B, barycenter);
221 LOG(2, "Volume is " << sign * volB);
224 LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
226 return _volume = sign * (volA + volB);
230 // ----------------------------------------------------------------------------------
231 // TransformedTriangle PRIVATE
232 // ----------------------------------------------------------------------------------
235 * Calculates the intersection polygons A and B, performing the intersection tests
236 * and storing the corresponding points in the vectors _polygonA and _polygonB.
238 * @post _polygonA contains the intersection polygon A and _polygonB contains the
239 * intersection polygon B.
242 void TransformedTriangle::calculateIntersectionPolygons()
244 assert(_polygonA.size() == 0);
245 assert(_polygonB.size() == 0);
246 // avoid reallocations in push_back() by pre-allocating enough memory
247 // we should never have more than 20 points
248 _polygonA.reserve(20);
249 _polygonB.reserve(20);
250 // -- surface intersections
252 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
254 if(testSurfaceEdgeIntersection(edge))
256 double* ptA = new double[3];
257 calcIntersectionPtSurfaceEdge(edge, ptA);
258 _polygonA.push_back(ptA);
259 LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
262 double* ptB = new double[3];
263 copyVector3(ptA, ptB);
264 _polygonB.push_back(ptB);
265 LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B ");
272 for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
274 if(testSurfaceRayIntersection(corner))
276 double* ptB = new double[3];
277 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
278 _polygonB.push_back(ptB);
279 LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
283 // -- segment intersections
284 for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
289 // check beforehand which double-products are zero
290 for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
292 isZero[dp] = (calcStableC(seg, dp) == 0.0);
296 for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
298 // is this test worth it?
300 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
301 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
302 !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
304 if(doTest && testSegmentFacetIntersection(seg, facet))
306 double* ptA = new double[3];
307 calcIntersectionPtSegmentFacet(seg, facet, ptA);
308 _polygonA.push_back(ptA);
309 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
312 double* ptB = new double[3];
313 copyVector3(ptA, ptB);
314 _polygonB.push_back(ptB);
315 LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
322 for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
324 const DoubleProduct edge_dp = DoubleProduct(edge);
326 if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
328 double* ptA = new double[3];
329 calcIntersectionPtSegmentEdge(seg, edge, ptA);
330 _polygonA.push_back(ptA);
331 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
334 double* ptB = new double[3];
335 copyVector3(ptA, ptB);
336 _polygonB.push_back(ptB);
342 for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
345 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
346 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
347 isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
349 if(doTest && testSegmentCornerIntersection(seg, corner))
351 double* ptA = new double[3];
352 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
353 _polygonA.push_back(ptA);
354 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
357 double* ptB = new double[3];
358 _polygonB.push_back(ptB);
359 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
360 LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
366 for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
368 if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
370 double* ptB = new double[3];
371 copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
372 _polygonB.push_back(ptB);
373 LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
377 // segment - halfstrip
378 for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
382 const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
384 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
385 !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
388 if(doTest && testSegmentHalfstripIntersection(seg, edge))
390 if(testSegmentHalfstripIntersection(seg, edge))
392 double* ptB = new double[3];
393 calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
394 _polygonB.push_back(ptB);
395 LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
401 for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
403 // { XYZ - inclusion only possible if in Tetrahedron?
405 if(testCornerInTetrahedron(corner))
407 double* ptA = new double[3];
408 copyVector3(&_coords[5*corner], ptA);
409 _polygonA.push_back(ptA);
410 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
414 if(testCornerOnXYZFacet(corner))
416 double* ptB = new double[3];
417 copyVector3(&_coords[5*corner], ptB);
418 _polygonB.push_back(ptB);
419 LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
422 // projection on XYZ - facet
423 if(testCornerAboveXYZFacet(corner))
425 double* ptB = new double[3];
426 copyVector3(&_coords[5*corner], ptB);
427 ptB[2] = 1 - ptB[0] - ptB[1];
428 assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
429 _polygonB.push_back(ptB);
430 LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
438 * Calculates the barycenters of the given intersection polygon.
440 * @pre the intersection polygons have been calculated with calculateIntersectionPolygons()
442 * @param poly one of the two intersection polygons
443 * @param barycenter array of three doubles where barycenter is stored
446 void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
448 LOG(3,"--- Calculating polygon barycenter");
450 // get the polygon points
451 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
453 // calculate barycenter
454 const int m = polygon.size();
456 for(int j = 0 ; j < 3 ; ++j)
463 for(int i = 0 ; i < m ; ++i)
465 const double* pt = polygon[i];
466 for(int j = 0 ; j < 3 ; ++j)
468 barycenter[j] += pt[j] / double(m);
472 LOG(3,"Barycenter is " << vToStr(barycenter));
476 * Sorts the given intersection polygon in circular order around its barycenter.
477 * @pre the intersection polygons have been calculated with calculateIntersectionPolygons()
478 * @post the vertices in _polygonA and _polygonB are sorted in circular order around their
479 * respective barycenters
481 * @param poly one of the two intersection polygons
482 * @param barycenter array of three doubles with the coordinates of the barycenter
485 void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
487 LOG(3,"--- Sorting polygon ...");
489 using INTERP_KERNEL::ProjectedCentralCircularSortOrder;
490 typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
491 typedef SortOrder::CoordType CoordType;
493 // get the polygon points
494 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
496 if(polygon.size() == 0)
499 // determine type of sorting
500 CoordType type = SortOrder::XY;
501 if(poly == A && !isTriangleInclinedToFacet( OXY )) // B is on h = 0 plane -> ok
503 // NB : the following test is never true if we have eliminated the
504 // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
505 // We keep the test here anyway, to avoid interdependency.
507 // is triangle inclined to x == 0 ?
508 if(isTriangleInclinedToFacet(OZX))
510 type = SortOrder::XZ;
512 else //if(isTriangleParallelToFacet(OYZ))
514 type = SortOrder::YZ;
518 // create order object
519 SortOrder order(barycenter, type);
521 // sort vector with this object
522 // NB : do not change place of first object, with respect to which the order
524 sort((polygon.begin()), polygon.end(), order);
526 LOG(3,"Sorted polygon is ");
527 for(size_t i = 0 ; i < polygon.size() ; ++i)
529 LOG(3,vToStr(polygon[i]));
535 * Calculates the volume between the given polygon and the z = 0 plane.
537 * @pre the intersection polygones have been calculated with calculateIntersectionPolygons(),
538 * and they have been sorted in circular order with sortIntersectionPolygons(void)
540 * @param poly one of the two intersection polygons
541 * @param barycenter array of three doubles with the coordinates of the barycenter
542 * @return the volume between the polygon and the z = 0 plane
545 double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter)
547 LOG(2,"--- Calculating volume under polygon");
549 // get the polygon points
550 std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
553 const int m = polygon.size();
555 for(int i = 0 ; i < m ; ++i)
557 const double* ptCurr = polygon[i]; // pt "i"
558 const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
560 const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
561 const double factor2 =
562 ptCurr[0]*(ptNext[1] - barycenter[1])
563 + ptNext[0]*(barycenter[1] - ptCurr[1])
564 + barycenter[0]*(ptCurr[1] - ptNext[1]);
565 vol += (factor1 * factor2) / 6.0;
568 LOG(2,"Abs. Volume is " << vol);
573 ////////////////////////////////////////////////////////////////////////////////////
574 // Detection of (very) degenerate cases /////////////
575 ////////////////////////////////////////////////////////////////////////////////////
578 * Checks if the triangle lies in the plane of a given facet
580 * @param facet one of the facets of the tetrahedron
581 * @return true if PQR lies in the plane of the facet, false if not
583 bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
586 // coordinate to check
587 const int coord = static_cast<int>(facet);
589 for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
591 if(_coords[5*c + coord] != 0.0)
601 * Checks if the triangle is parallel to the given facet
603 * @param facet one of the facets of the unit tetrahedron
604 * @return true if triangle is parallel to facet, false if not
606 bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const
608 // coordinate to check
609 const int coord = static_cast<int>(facet);
610 return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]);
614 * Checks if the triangle is not perpedicular to the given facet
616 * @param facet one of the facets of the unit tetrahedron
617 * @return zero if the triangle is perpendicular to the facet,
618 * else 1 or -1 depending on the sign of cross product of facet edges
620 int TransformedTriangle::isTriangleInclinedToFacet(const TetraFacet facet) const
622 // coordinate to check
623 const int coord = static_cast<int>(facet);
624 const int ind1 = ( coord+1 ) % 3, ind2 = ( coord+2 ) % 3;
625 const double uv_xy[4] =
628 _coords[5*Q+ind1] - _coords[5*P+ind1], _coords[5*Q+ind2] - _coords[5*P+ind2],
630 _coords[5*R+ind1] - _coords[5*P+ind1], _coords[5*R+ind2] - _coords[5*P+ind2]
633 double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
634 return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0;
638 * Determines whether the triangle is below the z-plane.
640 * @return true if the z-coordinate of the three corners of the triangle are all less than 0, false otherwise.
642 bool TransformedTriangle::isTriangleBelowTetraeder() const
644 for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
646 // check z-coords for all points
647 if(_coords[5*c + 2] >= 0.0)
656 * Prints the coordinates of the triangle to std::cout
659 void TransformedTriangle::dumpCoords() const
661 std::cout << "Coords : ";
662 for(int i = 0 ; i < 3; ++i)
664 std::cout << vToStr(&_coords[5*i]) << ",";
666 std::cout << std::endl;