mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
b4a482751b
- combine header definitions, more pass-through methods - face/triFace: support += operator (vertex offset)
417 lines
9.9 KiB
C
417 lines
9.9 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2011 OpenFOAM Foundation
|
|
Copyright (C) 2017 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "triangleFuncs.H"
|
|
#include "triangle.H"
|
|
#include "pointField.H"
|
|
#include "treeBoundBox.H"
|
|
#include "SortableList.H"
|
|
#include "boolList.H"
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
void Foam::triangleFuncs::setIntersection
|
|
(
|
|
const point& oppositeSidePt,
|
|
const scalar oppositeSign,
|
|
|
|
const point& thisSidePt,
|
|
const scalar thisSign,
|
|
|
|
const scalar tol,
|
|
|
|
point& pt
|
|
)
|
|
{
|
|
const scalar denom = oppositeSign - thisSign;
|
|
|
|
if (mag(denom) < tol)
|
|
{
|
|
// If almost does not cut choose one which certainly cuts.
|
|
pt = oppositeSidePt;
|
|
}
|
|
else
|
|
{
|
|
pt = oppositeSidePt + oppositeSign/denom*(thisSidePt - oppositeSidePt);
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
bool Foam::triangleFuncs::intersectAxesBundle
|
|
(
|
|
const point& V0,
|
|
const point& V10,
|
|
const point& V20,
|
|
const label i0,
|
|
const pointField& origin,
|
|
const scalar maxLength,
|
|
point& pInter
|
|
)
|
|
{
|
|
// Based on Graphics Gems - Fast Ray Triangle intersection.
|
|
// Since direction is coordinate axis there is no need to do projection,
|
|
// we can directly check u,v components for inclusion in triangle.
|
|
|
|
// Get other components
|
|
label i1 = (i0 + 1) % 3;
|
|
label i2 = (i1 + 1) % 3;
|
|
|
|
scalar u1 = V10[i1];
|
|
scalar v1 = V10[i2];
|
|
|
|
scalar u2 = V20[i1];
|
|
scalar v2 = V20[i2];
|
|
|
|
scalar localScale = mag(u1)+mag(v1)+mag(u2)+mag(v2);
|
|
|
|
scalar det = v2*u1 - u2*v1;
|
|
|
|
// Fix for V0:(-31.71428 0 -15.10714)
|
|
// V10:(-1.285715 8.99165e-16 -1.142858)
|
|
// V20:(0 0 -1.678573)
|
|
// i0:0
|
|
if (localScale < VSMALL || Foam::mag(det)/localScale < SMALL)
|
|
{
|
|
// Triangle parallel to dir
|
|
return false;
|
|
}
|
|
|
|
forAll(origin, originI)
|
|
{
|
|
const point& P = origin[originI];
|
|
|
|
scalar u0 = P[i1] - V0[i1];
|
|
scalar v0 = P[i2] - V0[i2];
|
|
|
|
scalar alpha = 0;
|
|
scalar beta = 0;
|
|
bool inter = false;
|
|
|
|
if (Foam::mag(u1) < ROOTVSMALL)
|
|
{
|
|
beta = u0/u2;
|
|
if ((beta >= 0) && (beta <= 1))
|
|
{
|
|
alpha = (v0 - beta*v2)/v1;
|
|
inter = ((alpha >= 0) && ((alpha + beta) <= 1));
|
|
}
|
|
}
|
|
else
|
|
{
|
|
beta = (v0*u1 - u0*v1)/det;
|
|
if ((beta >= 0) && (beta <= 1))
|
|
{
|
|
alpha = (u0 - beta*u2)/u1;
|
|
inter = ((alpha >= 0) && ((alpha + beta) <= 1));
|
|
}
|
|
}
|
|
|
|
if (inter)
|
|
{
|
|
pInter = V0 + alpha*V10 + beta*V20;
|
|
scalar s = (pInter - origin[originI])[i0];
|
|
|
|
if ((s >= 0) && (s <= maxLength))
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
|
|
bool Foam::triangleFuncs::intersectBb
|
|
(
|
|
const point& p0,
|
|
const point& p1,
|
|
const point& p2,
|
|
const treeBoundBox& cubeBb
|
|
)
|
|
{
|
|
// Slow (edge by edge) bounding box intersection. TBD: replace with call
|
|
// to above intersectAxesBundle. However this function is not fully
|
|
// correct and misses intersection between some triangles.
|
|
{
|
|
const triPointRef tri(p0, p1, p2);
|
|
|
|
const edgeList& es = treeBoundBox::edges;
|
|
const pointField points(cubeBb.points());
|
|
|
|
forAll(es, i)
|
|
{
|
|
const edge& e = es[i];
|
|
const point& start = points[e[0]];
|
|
const point& end = points[e[1]];
|
|
|
|
pointHit inter = tri.intersection
|
|
(
|
|
start,
|
|
end-start,
|
|
intersection::HALF_RAY
|
|
);
|
|
|
|
if (inter.hit() && inter.distance() <= 1)
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Intersect triangle edges with bounding box
|
|
point pInter;
|
|
if (cubeBb.intersects(p0, p1, pInter))
|
|
{
|
|
return true;
|
|
}
|
|
if (cubeBb.intersects(p1, p2, pInter))
|
|
{
|
|
return true;
|
|
}
|
|
if (cubeBb.intersects(p2, p0, pInter))
|
|
{
|
|
return true;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
bool Foam::triangleFuncs::intersect
|
|
(
|
|
const point& va0,
|
|
const point& va10,
|
|
const point& va20,
|
|
|
|
const point& base,
|
|
const point& normal,
|
|
|
|
point& pInter0,
|
|
point& pInter1
|
|
)
|
|
{
|
|
// Get triangle normal
|
|
vector na = va10 ^ va20;
|
|
scalar magArea = mag(na);
|
|
na/magArea;
|
|
|
|
if (mag(na & normal) > (1 - SMALL))
|
|
{
|
|
// Parallel
|
|
return false;
|
|
}
|
|
|
|
const point va1 = va0 + va10;
|
|
const point va2 = va0 + va20;
|
|
|
|
// Find the triangle point on the other side.
|
|
scalar sign0 = (va0 - base) & normal;
|
|
scalar sign1 = (va1 - base) & normal;
|
|
scalar sign2 = (va2 - base) & normal;
|
|
|
|
label oppositeVertex = -1;
|
|
|
|
if (sign0 < 0)
|
|
{
|
|
if (sign1 < 0)
|
|
{
|
|
if (sign2 < 0)
|
|
{
|
|
// All on same side of plane
|
|
return false;
|
|
}
|
|
else // sign2 >= 0
|
|
{
|
|
// 2 on opposite side.
|
|
oppositeVertex = 2;
|
|
}
|
|
}
|
|
else // sign1 >= 0
|
|
{
|
|
if (sign2 < 0)
|
|
{
|
|
// 1 on opposite side.
|
|
oppositeVertex = 1;
|
|
}
|
|
else
|
|
{
|
|
// 0 on opposite side.
|
|
oppositeVertex = 0;
|
|
}
|
|
}
|
|
}
|
|
else // sign0 >= 0
|
|
{
|
|
if (sign1 < 0)
|
|
{
|
|
if (sign2 < 0)
|
|
{
|
|
// 0 on opposite side.
|
|
oppositeVertex = 0;
|
|
}
|
|
else // sign2 >= 0
|
|
{
|
|
// 1 on opposite side.
|
|
oppositeVertex = 1;
|
|
}
|
|
}
|
|
else // sign1 >= 0
|
|
{
|
|
if (sign2 < 0)
|
|
{
|
|
// 2 on opposite side.
|
|
oppositeVertex = 2;
|
|
}
|
|
else // sign2 >= 0
|
|
{
|
|
// All on same side of plane
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
|
|
scalar tol = SMALL*Foam::sqrt(magArea);
|
|
|
|
if (oppositeVertex == 0)
|
|
{
|
|
// 0 on opposite side. Cut edges 01 and 02
|
|
setIntersection(va0, sign0, va1, sign1, tol, pInter0);
|
|
setIntersection(va0, sign0, va2, sign2, tol, pInter1);
|
|
}
|
|
else if (oppositeVertex == 1)
|
|
{
|
|
// 1 on opposite side. Cut edges 10 and 12
|
|
setIntersection(va1, sign1, va0, sign0, tol, pInter0);
|
|
setIntersection(va1, sign1, va2, sign2, tol, pInter1);
|
|
}
|
|
else // oppositeVertex == 2
|
|
{
|
|
// 2 on opposite side. Cut edges 20 and 21
|
|
setIntersection(va2, sign2, va0, sign0, tol, pInter0);
|
|
setIntersection(va2, sign2, va1, sign1, tol, pInter1);
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
bool Foam::triangleFuncs::intersect
|
|
(
|
|
const point& va0,
|
|
const point& va10,
|
|
const point& va20,
|
|
|
|
const point& vb0,
|
|
const point& vb10,
|
|
const point& vb20,
|
|
|
|
point& pInter0,
|
|
point& pInter1
|
|
)
|
|
{
|
|
// Get triangle normals
|
|
vector na = va10 ^ va20;
|
|
na/mag(na);
|
|
|
|
vector nb = vb10 ^ vb20;
|
|
nb/mag(nb);
|
|
|
|
// Calculate intersection of triangle a with plane of b
|
|
point planeB0;
|
|
point planeB1;
|
|
if (!intersect(va0, va10, va20, vb0, nb, planeB0, planeB1))
|
|
{
|
|
return false;
|
|
}
|
|
|
|
// ,, triangle b with plane of a
|
|
point planeA0;
|
|
point planeA1;
|
|
if (!intersect(vb0, vb10, vb20, va0, na, planeA0, planeA1))
|
|
{
|
|
return false;
|
|
}
|
|
|
|
// Now check if intersections overlap (w.r.t. intersection of the two
|
|
// planes)
|
|
|
|
vector intersection(na ^ nb);
|
|
|
|
scalar coordB0 = planeB0 & intersection;
|
|
scalar coordB1 = planeB1 & intersection;
|
|
|
|
scalar coordA0 = planeA0 & intersection;
|
|
scalar coordA1 = planeA1 & intersection;
|
|
|
|
// Put information in indexable form for sorting.
|
|
List<point*> pts(4);
|
|
boolList isFromB(4);
|
|
SortableList<scalar> sortCoords(4);
|
|
|
|
pts[0] = &planeB0;
|
|
isFromB[0] = true;
|
|
sortCoords[0] = coordB0;
|
|
|
|
pts[1] = &planeB1;
|
|
isFromB[1] = true;
|
|
sortCoords[1] = coordB1;
|
|
|
|
pts[2] = &planeA0;
|
|
isFromB[2] = false;
|
|
sortCoords[2] = coordA0;
|
|
|
|
pts[3] = &planeA1;
|
|
isFromB[3] = false;
|
|
sortCoords[3] = coordA1;
|
|
|
|
sortCoords.sort();
|
|
|
|
const labelList& indices = sortCoords.indices();
|
|
|
|
if (isFromB[indices[0]] == isFromB[indices[1]])
|
|
{
|
|
// Entry 0 and 1 are of same region (both a or both b). Hence that
|
|
// region does not overlap.
|
|
return false;
|
|
}
|
|
else
|
|
{
|
|
// Different type. Start of overlap at indices[1], end at indices[2]
|
|
pInter0 = *pts[indices[1]];
|
|
pInter1 = *pts[indices[2]];
|
|
|
|
return true;
|
|
}
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|