1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 \*---------------------------------------------------------------------------*/
29 #include "triangleFuncs.H"
30 #include "pointField.H"
31 #include "treeBoundBox.H"
32 #include "SortableList.H"
37 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 void Foam::triangleFuncs::setIntersection
41 const point& oppositeSidePt,
42 const scalar oppositeSign,
44 const point& thisSidePt,
45 const scalar thisSign,
52 scalar denom = oppositeSign - thisSign;
56 // If almost does not cut choose one which certainly cuts.
61 pt = oppositeSidePt + oppositeSign/denom*(thisSidePt - oppositeSidePt);
66 void Foam::triangleFuncs::selectPt
85 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
87 // Intersect triangle with parallel edges aligned with axis i0.
88 // Returns true (and intersection in pInter) if any of them intersects triangle.
89 bool Foam::triangleFuncs::intersectAxesBundle
95 const pointField& origin,
96 const scalar maxLength,
100 // Based on Graphics Gems - Fast Ray Triangle intersection.
101 // Since direction is coordinate axis there is no need to do projection,
102 // we can directly check u,v components for inclusion in triangle.
104 scalar localScale = max(max(magSqr(V10), magSqr(V20)), 1.0);
106 // Get other components
107 label i1 = (i0 + 1) % 3;
108 label i2 = (i1 + 1) % 3;
116 scalar det = v2*u1 - u2*v1;
118 // Fix for V0:(-31.71428 0 -15.10714)
119 // V10:(-1.285715 8.99165e-16 -1.142858)
120 // V20:(0 0 -1.678573)
122 if (Foam::mag(det)/localScale < SMALL)
124 // Triangle parallel to dir
128 forAll(origin, originI)
130 const point& P = origin[originI];
132 scalar u0 = P[i1] - V0[i1];
133 scalar v0 = P[i2] - V0[i2];
139 if (Foam::mag(u1) < ROOTVSMALL)
142 if ((beta >= 0) && (beta <= 1))
144 alpha = (v0 - beta*v2)/v1;
145 inter = ((alpha >= 0) && ((alpha + beta) <= 1));
150 beta = (v0*u1 - u0*v1)/det;
151 if ((beta >= 0) && (beta <= 1))
153 alpha = (u0 - beta*u2)/u1;
154 inter = ((alpha >= 0) && ((alpha + beta) <= 1));
160 pInter = V0 + alpha*V10 + beta*V20;
161 scalar s = (pInter - origin[originI])[i0];
163 if ((s >= 0) && (s <= maxLength))
174 // Intersect triangle with bounding box. Return true if
175 // any of the faces of bb intersect triangle.
176 // Note: so returns false if triangle inside bb.
177 bool Foam::triangleFuncs::intersectBb
182 const treeBoundBox& cubeBb
185 const vector p10 = p1 - p0;
186 const vector p20 = p2 - p0;
188 // cubeBb points; counted as if cell with vertex0 at cubeBb.min().
189 const point& min = cubeBb.min();
190 const point& max = cubeBb.max();
192 const point& cube0 = min;
193 const point cube1(min.x(), min.y(), max.z());
194 const point cube2(max.x(), min.y(), max.z());
195 const point cube3(max.x(), min.y(), min.z());
197 const point cube4(min.x(), max.y(), min.z());
198 const point cube5(min.x(), max.y(), max.z());
199 const point cube7(max.x(), max.y(), min.z());
202 // Intersect all 12 edges of cube with triangle
206 pointField origin(4);
207 // edges in x direction
213 scalar maxSx = max.x() - min.x();
215 if (intersectAxesBundle(p0, p10, p20, 0, origin, maxSx, pInter))
220 // edges in y direction
226 scalar maxSy = max.y() - min.y();
228 if (intersectAxesBundle(p0, p10, p20, 1, origin, maxSy, pInter))
233 // edges in z direction
239 scalar maxSz = max.z() - min.z();
241 if (intersectAxesBundle(p0, p10, p20, 2, origin, maxSz, pInter))
247 // Intersect triangle edges with bounding box
248 if (cubeBb.intersects(p0, p1, pInter))
252 if (cubeBb.intersects(p1, p2, pInter))
256 if (cubeBb.intersects(p2, p0, pInter))
265 bool Foam::triangleFuncs::intersect
278 // Get triangle normal
279 vector na = va10 ^ va20;
280 scalar magArea = mag(na);
283 if (mag(na & normal) > (1 - SMALL))
289 const point va1 = va0 + va10;
290 const point va2 = va0 + va20;
292 // Find the triangle point on the other side.
293 scalar sign0 = (va0 - base) & normal;
294 scalar sign1 = (va1 - base) & normal;
295 scalar sign2 = (va2 - base) & normal;
297 label oppositeVertex = -1;
305 // All on same side of plane
310 // 2 on opposite side.
318 // 1 on opposite side.
323 // 0 on opposite side.
334 // 0 on opposite side.
339 // 1 on opposite side.
347 // 2 on opposite side.
352 // All on same side of plane
358 scalar tol = SMALL*Foam::sqrt(magArea);
360 if (oppositeVertex == 0)
362 // 0 on opposite side. Cut edges 01 and 02
363 setIntersection(va0, sign0, va1, sign1, tol, pInter0);
364 setIntersection(va0, sign0, va2, sign2, tol, pInter1);
366 else if (oppositeVertex == 1)
368 // 1 on opposite side. Cut edges 10 and 12
369 setIntersection(va1, sign1, va0, sign0, tol, pInter0);
370 setIntersection(va1, sign1, va2, sign2, tol, pInter1);
372 else // oppositeVertex == 2
374 // 2 on opposite side. Cut edges 20 and 21
375 setIntersection(va2, sign2, va0, sign0, tol, pInter0);
376 setIntersection(va2, sign2, va1, sign1, tol, pInter1);
383 bool Foam::triangleFuncs::intersect
397 // Get triangle normals
398 vector na = va10 ^ va20;
401 vector nb = vb10 ^ vb20;
404 // Calculate intersection of triangle a with plane of b
407 if (!intersect(va0, va10, va20, vb0, nb, planeB0, planeB1))
412 // ,, triangle b with plane of a
415 if (!intersect(vb0, vb10, vb20, va0, na, planeA0, planeA1))
420 // Now check if intersections overlap (w.r.t. intersection of the two
423 vector intersection(na ^ nb);
425 scalar coordB0 = planeB0 & intersection;
426 scalar coordB1 = planeB1 & intersection;
428 scalar coordA0 = planeA0 & intersection;
429 scalar coordA1 = planeA1 & intersection;
431 // Put information in indexable form for sorting.
434 SortableList<scalar> sortCoords(4);
438 sortCoords[0] = coordB0;
442 sortCoords[1] = coordB1;
446 sortCoords[2] = coordA0;
450 sortCoords[3] = coordA1;
454 const labelList& indices = sortCoords.indices();
456 if (isFromB[indices[0]] == isFromB[indices[1]])
458 // Entry 0 and 1 are of same region (both a or both b). Hence that
459 // region does not overlap.
464 // Different type. Start of overlap at indices[1], end at indices[2]
465 pInter0 = *pts[indices[1]];
466 pInter1 = *pts[indices[2]];
473 bool Foam::triangleFuncs::classify
475 const point& baseVertex,
485 // Get largest component of normal
486 scalar magX = Foam::mag(n.x());
487 scalar magY = Foam::mag(n.y());
488 scalar magZ = Foam::mag(n.z());
491 if ((magX >= magY) && (magX >= magZ))
495 else if ((magY >= magX) && (magY >= magZ))
504 // Get other components
505 label i1 = (i0 + 1) % 3;
506 label i2 = (i1 + 1) % 3;
514 scalar det = v2*u1 - u2*v1;
516 scalar u0 = pInter[i1] - baseVertex[i1];
517 scalar v0 = pInter[i2] - baseVertex[i2];
524 if (Foam::mag(u1) < ROOTVSMALL)
528 alpha = (v0 - beta*v2)/v1;
530 hit = ((alpha >= 0) && ((alpha + beta) <= 1));
534 beta = (v0*u1 - u0*v1)/det;
536 alpha = (u0 - beta*u2)/u1;
538 hit = ((alpha >= 0) && ((alpha + beta) <= 1));
542 // Now alpha, beta are the coordinates in the triangle local coordinate
549 if (mag(alpha+beta-1) < tol)
551 // On edge between vert 1-2 (=edge 1)
553 if (mag(alpha) < tol)
558 else if (mag(beta) < tol)
563 else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
569 else if (mag(alpha) < tol)
571 // On edge between vert 2-0 (=edge 2)
578 else if (mag(beta-1) < tol)
583 else if ((beta >= 0) && (beta <= 1))
589 else if (mag(beta) < tol)
591 // On edge between vert 0-1 (= edge 0)
593 if (mag(alpha) < tol)
598 else if (mag(alpha-1) < tol)
603 else if ((alpha >= 0) && (alpha <= 1))
614 // ************************************************************************* //