/usr/include/freefoam/finiteVolume/findCellPointFaceTriangle.H is in libfreefoam-dev 0.1.0+dfsg-1build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 | /*
find the triangle in which the position is,
when the position is known to be on the face
*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
bool interpolationCellPointFace<Type>::findTriangle
(
const vector& position,
const label nFace,
label tetPointLabels[],
scalar phi[]
) const
{
bool foundTriangle = false;
vector tetPoints[3];
const labelList& facePoints = this->pMeshFaces_[nFace];
tetPoints[2] = this->pMeshFaceCentres_[nFace];
label pointi = 0;
while (pointi < facePoints.size() && !foundTriangle)
{
// The triangle is constructed from face center and one
// face edge
label nextPointLabel = (pointi + 1) % facePoints.size();
tetPointLabels[0] = facePoints[pointi];
tetPointLabels[1] = facePoints[nextPointLabel];
tetPoints[0] = this->pMeshPoints_[tetPointLabels[0]];
tetPoints[1] = this->pMeshPoints_[tetPointLabels[1]];
vector fc = (tetPoints[0] + tetPoints[1] + tetPoints[2])/3.0;
vector newPos = position + SMALL*(fc-position);
// calculate triangle edge vectors and triangle face normal
// the 'i':th edge is opposite node i
vector edge[3], normal[3];
for(label i=0; i<3; i++)
{
label ip0 = (i+1) % 3;
label ipp = (i+2) % 3;
edge[i] = tetPoints[ipp]-tetPoints[ip0];
}
vector triangleFaceNormal = edge[1] ^ edge[2];
// calculate edge normal (pointing inwards)
for(label i=0; i<3; i++)
{
normal[i] = triangleFaceNormal ^ edge[i];
normal[i] /= mag(normal[i]) + VSMALL;
}
// check if position is inside triangle
bool inside = true;
for(label i=0; i<3; i++)
{
label ip = (i+1) % 3;
inside = inside && (((newPos - tetPoints[ip]) & edge[i]) >= 0);
}
if (inside)
{
foundTriangle = true;
// calculate phi-values
for(label i=0; i<3; i++)
{
label ip = (i+1) % 3;
scalar phiMax = max(VSMALL, normal[i] & edge[ip]);
scalar phiLength = (position-tetPoints[ip]) & normal[i];
phi[i] = phiLength/phiMax;
}
}
pointi++;
}
return foundTriangle;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************ vim: set sw=4 sts=4 et: ************************ //
|