45 Pout<<
"\nFoam::cellPointWeight::findTetrahedron" <<
nl
46 <<
"position = " << position <<
nl
47 <<
"cellIndex = " << cellIndex <<
endl;
51 scalar minUVWClose = VGREAT;
52 label pointIClose = 0;
57 const scalar cellVolume = mesh.
cellVolumes()[cellIndex];
67 while ((pointI + 1) < facePoints.
size())
72 const vector& P3 = mesh.
points()[facePoints[pointI + 1]];
82 const vector rhs = position - P0;
87 e1.
x()*e2.
y()*e3.
z() + e2.
x()*e3.
y()*e1.
z()
88 + e3.
x()*e1.
y()*e2.
z() - e1.
x()*e3.
y()*e2.
z()
89 - e2.
x()*e1.
y()*e3.
z() - e3.
x()*e2.
y()*e1.
z();
91 if (
mag(detA/cellVolume) > tol)
96 rhs.
x()*e2.
y()*e3.
z() + e2.
x()*e3.
y()*rhs.
z()
97 + e3.
x()*rhs.
y()*e2.
z() - rhs.
x()*e3.
y()*e2.
z()
98 - e2.
x()*rhs.
y()*e3.
z() - e3.
x()*e2.
y()*rhs.
z()
103 e1.
x()*rhs.
y()*e3.
z() + rhs.
x()*e3.
y()*e1.
z()
104 + e3.
x()*e1.
y()*rhs.
z() - e1.
x()*e3.
y()*rhs.
z()
105 - rhs.
x()*e1.
y()*e3.
z() - e3.
x()*rhs.
y()*e1.
z()
110 e1.
x()*e2.
y()*rhs.
z() + e2.
x()*rhs.
y()*e1.
z()
111 + rhs.
x()*e1.
y()*e2.
z() - e1.
x()*rhs.
y()*e2.
z()
112 - e2.
x()*e1.
y()*rhs.
z() - rhs.
x()*e2.
y()*e1.
z()
119 (u + tol > 0) && (v + tol > 0) && (w + tol > 0)
120 && (u + v + w < 1 + tol)
123 faceVertices_[0] = facePoints[0];
124 faceVertices_[1] = facePoints[pointI];
125 faceVertices_[2] = facePoints[pointI + 1];
130 weights_[3] = 1.0 - (u + v + w);
136 scalar minU =
mag(u);
137 scalar minV =
mag(v);
138 scalar minW =
mag(w);
151 const scalar minUVW =
mag(minU + minV + minW);
153 if (minUVW < minUVWClose)
155 minUVWClose = minUVW;
156 pointIClose = pointI;
168 Pout<<
"cellPointWeight::findTetrahedron" <<
nl
169 <<
" Tetrahedron search failed; using closest tet values to "
170 <<
"point " <<
nl <<
" cell: " << cellIndex <<
nl <<
endl;
173 const labelList& facePointsClose = mesh.
faces()[cellFaces[faceClose]];
174 faceVertices_[0] = facePointsClose[0];
175 faceVertices_[1] = facePointsClose[pointIClose];
176 faceVertices_[2] = facePointsClose[pointIClose + 1];
189 const label faceIndex
194 Pout<<
"\nbool Foam::cellPointWeight::findTriangle" <<
nl
195 <<
"position = " << position <<
nl
196 <<
"faceIndex = " << faceIndex <<
endl;
200 scalar minUVClose = VGREAT;
201 label pointIClose = 0;
210 while ((pointI + 1) < facePoints.
size())
215 const vector& P3 = mesh.
points()[facePoints[pointI + 1]];
220 const vector v3 = P3 - P1;
230 const scalar d12 = v1 &
v2;
231 const scalar d13 = v1 & v3;
232 const scalar d22 = v2 &
v2;
233 const scalar d23 = v2 & v3;
234 const scalar d33 = v3 & v3;
238 const scalar detA = d22*d33 - d23*d23;
240 if (0.25*detA/faceArea2 > tol)
243 const scalar u = (d12*d33 - d23*d13)/detA;
244 const scalar v = (d22*d13 - d12*d23)/detA;
247 if ((u + tol > 0) && (v + tol > 0) && (u + v < 1 + tol))
250 faceVertices_[0] = facePoints[0];
251 faceVertices_[1] = facePoints[pointI];
252 faceVertices_[2] = facePoints[pointI + 1];
256 weights_[2] = 1.0 - (u + v);
263 scalar minU =
mag(u);
264 scalar minV =
mag(v);
273 const scalar minUV =
mag(minU + minV);
275 if (minUV < minUVClose)
278 pointIClose = pointI;
288 Pout<<
"Foam::cellPointWeight::findTriangle"
289 <<
"Triangle search failed; using closest triangle to point" <<
nl
290 <<
" cell face: " << faceIndex <<
nl <<
endl;
294 faceVertices_[0] = facePoints[0];
295 faceVertices_[1] = facePoints[pointIClose];
296 faceVertices_[2] = facePoints[pointIClose + 1];
298 weights_[0] = 1.0/3.0;
299 weights_[1] = 1.0/3.0;
300 weights_[2] = 1.0/3.0;
311 const label cellIndex,
312 const label faceIndex
315 cellIndex_(cellIndex)
320 findTetrahedron(mesh, position, cellIndex);
325 findTriangle(mesh, position, faceIndex);