FreeFOAM The Cross-Platform CFD Toolkit
cellPointWeight.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
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
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "cellPointWeight.H"
27 #include <OpenFOAM/polyMesh.H>
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 int Foam::cellPointWeight::debug(debug::debugSwitch("cellPointWeight", 0));
32 Foam::scalar Foam::cellPointWeight::tol(SMALL);
33 
34 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
35 
37 (
38  const polyMesh& mesh,
39  const vector& position,
40  const label cellIndex
41 )
42 {
43  if (debug)
44  {
45  Pout<< "\nFoam::cellPointWeight::findTetrahedron" << nl
46  << "position = " << position << nl
47  << "cellIndex = " << cellIndex << endl;
48  }
49 
50  // Initialise closest triangle variables
51  scalar minUVWClose = VGREAT;
52  label pointIClose = 0;
53  label faceClose = 0;
54 
55  const vector& P0 = mesh.cellCentres()[cellIndex];
56  const labelList& cellFaces = mesh.cells()[cellIndex];
57  const scalar cellVolume = mesh.cellVolumes()[cellIndex];
58 
59  // Find the tet that the point occupies
60  forAll(cellFaces, faceI)
61  {
62  // Decompose each face into triangles, making a tet when
63  // augmented by the cell centre
64  const labelList& facePoints = mesh.faces()[cellFaces[faceI]];
65 
66  label pointI = 1;
67  while ((pointI + 1) < facePoints.size())
68  {
69  // Cartesian co-ordinates of the triangle vertices
70  const vector& P1 = mesh.points()[facePoints[0]];
71  const vector& P2 = mesh.points()[facePoints[pointI]];
72  const vector& P3 = mesh.points()[facePoints[pointI + 1]];
73 
74  // Edge vectors
75  const vector e1 = P1 - P0;
76  const vector e2 = P2 - P0;
77  const vector e3 = P3 - P0;
78 
79  // Solve for interpolation weighting factors
80 
81  // Source term
82  const vector rhs = position - P0;
83 
84  // Determinant of coefficients matrix
85  // Note: if det(A) = 0 the tet is degenerate
86  const scalar detA =
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();
90 
91  if (mag(detA/cellVolume) > tol)
92  {
93  // Solve using Cramers' rule
94  const scalar u =
95  (
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()
99  )/detA;
100 
101  const scalar v =
102  (
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()
106  )/detA;
107 
108  const scalar w =
109  (
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()
113  )/detA;
114 
115  // Check if point is in tet
116  // value = 0 indicates position lies on a tet face
117  if
118  (
119  (u + tol > 0) && (v + tol > 0) && (w + tol > 0)
120  && (u + v + w < 1 + tol)
121  )
122  {
123  faceVertices_[0] = facePoints[0];
124  faceVertices_[1] = facePoints[pointI];
125  faceVertices_[2] = facePoints[pointI + 1];
126 
127  weights_[0] = u;
128  weights_[1] = v;
129  weights_[2] = w;
130  weights_[3] = 1.0 - (u + v + w);
131 
132  return;
133  }
134  else
135  {
136  scalar minU = mag(u);
137  scalar minV = mag(v);
138  scalar minW = mag(w);
139  if (minU > 1.0)
140  {
141  minU -= 1.0;
142  }
143  if (minV > 1.0)
144  {
145  minV -= 1.0;
146  }
147  if (minW > 1.0)
148  {
149  minW -= 1.0;
150  }
151  const scalar minUVW = mag(minU + minV + minW);
152 
153  if (minUVW < minUVWClose)
154  {
155  minUVWClose = minUVW;
156  pointIClose = pointI;
157  faceClose = faceI;
158  }
159  }
160  }
161 
162  pointI++;
163  }
164  }
165 
166  if (debug)
167  {
168  Pout<< "cellPointWeight::findTetrahedron" << nl
169  << " Tetrahedron search failed; using closest tet values to "
170  << "point " << nl << " cell: " << cellIndex << nl << endl;
171  }
172 
173  const labelList& facePointsClose = mesh.faces()[cellFaces[faceClose]];
174  faceVertices_[0] = facePointsClose[0];
175  faceVertices_[1] = facePointsClose[pointIClose];
176  faceVertices_[2] = facePointsClose[pointIClose + 1];
177 
178  weights_[0] = 0.25;
179  weights_[1] = 0.25;
180  weights_[2] = 0.25;
181  weights_[3] = 0.25;
182 }
183 
184 
186 (
187  const polyMesh& mesh,
188  const vector& position,
189  const label faceIndex
190 )
191 {
192  if (debug)
193  {
194  Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
195  << "position = " << position << nl
196  << "faceIndex = " << faceIndex << endl;
197  }
198 
199  // Initialise closest triangle variables
200  scalar minUVClose = VGREAT;
201  label pointIClose = 0;
202 
203  // Decompose each face into triangles, making a tet when
204  // augmented by the cell centre
205  const labelList& facePoints = mesh.faces()[faceIndex];
206 
207  const scalar faceArea2 = magSqr(mesh.faceAreas()[faceIndex]);
208 
209  label pointI = 1;
210  while ((pointI + 1) < facePoints.size())
211  {
212  // Cartesian co-ordinates of the triangle vertices
213  const vector& P1 = mesh.points()[facePoints[0]];
214  const vector& P2 = mesh.points()[facePoints[pointI]];
215  const vector& P3 = mesh.points()[facePoints[pointI + 1]];
216 
217  // Direction vectors
218  vector v1 = position - P1;
219  const vector v2 = P2 - P1;
220  const vector v3 = P3 - P1;
221 
222  // Plane normal
223  vector n = v2 ^ v3;
224  n /= mag(n);
225 
226  // Remove any offset to plane
227  v1 -= (n & v1)*v1;
228 
229  // Helper variables
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;
235 
236  // Determinant of coefficients matrix
237  // Note: if det(A) = 0 the triangle is degenerate
238  const scalar detA = d22*d33 - d23*d23;
239 
240  if (0.25*detA/faceArea2 > tol)
241  {
242  // Solve using Cramers' rule
243  const scalar u = (d12*d33 - d23*d13)/detA;
244  const scalar v = (d22*d13 - d12*d23)/detA;
245 
246  // Check if point is in triangle
247  if ((u + tol > 0) && (v + tol > 0) && (u + v < 1 + tol))
248  {
249  // Indices of the cell vertices making up the triangle
250  faceVertices_[0] = facePoints[0];
251  faceVertices_[1] = facePoints[pointI];
252  faceVertices_[2] = facePoints[pointI + 1];
253 
254  weights_[0] = u;
255  weights_[1] = v;
256  weights_[2] = 1.0 - (u + v);
257  weights_[3] = 0.0;
258 
259  return;
260  }
261  else
262  {
263  scalar minU = mag(u);
264  scalar minV = mag(v);
265  if (minU > 1.0)
266  {
267  minU -= 1.0;
268  }
269  if (minV > 1.0)
270  {
271  minV -= 1.0;
272  }
273  const scalar minUV = mag(minU + minV);
274 
275  if (minUV < minUVClose)
276  {
277  minUVClose = minUV;
278  pointIClose = pointI;
279  }
280  }
281  }
282 
283  pointI++;
284  }
285 
286  if (debug)
287  {
288  Pout<< "Foam::cellPointWeight::findTriangle"
289  << "Triangle search failed; using closest triangle to point" << nl
290  << " cell face: " << faceIndex << nl << endl;
291  }
292 
293  // Indices of the cell vertices making up the triangle
294  faceVertices_[0] = facePoints[0];
295  faceVertices_[1] = facePoints[pointIClose];
296  faceVertices_[2] = facePoints[pointIClose + 1];
297 
298  weights_[0] = 1.0/3.0;
299  weights_[1] = 1.0/3.0;
300  weights_[2] = 1.0/3.0;
301  weights_[3] = 0.0;
302 }
303 
304 
305 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
306 
308 (
309  const polyMesh& mesh,
310  const vector& position,
311  const label cellIndex,
312  const label faceIndex
313 )
314 :
315  cellIndex_(cellIndex)
316 {
317  if (faceIndex < 0)
318  {
319  // Face data not supplied
320  findTetrahedron(mesh, position, cellIndex);
321  }
322  else
323  {
324  // Face data supplied
325  findTriangle(mesh, position, faceIndex);
326  }
327 }
328 
329 
330 // ************************ vim: set sw=4 sts=4 et: ************************ //