FreeFOAM The Cross-Platform CFD Toolkit
interpolationCellPointFace.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 Description
25 
26 \*---------------------------------------------------------------------------*/
27 
29 #include <finiteVolume/volFields.H>
30 #include <OpenFOAM/polyMesh.H>
32 #include <finiteVolume/linear.H>
33 #include "findCellPointFaceTet.H"
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 
41 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
42 
43 template<class Type>
45 (
47 )
48 :
50  psip_(volPointInterpolation::New(psi.mesh()).interpolate(psi)),
51  psis_(linearInterpolate(psi))
52 {}
53 
54 
55 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
56 
57 template<class Type>
59 (
60  const vector& position,
61  const label nCell,
62  const label facei
63 ) const
64 {
65  Type ts[4];
66  vector tetPoints[4];
67  scalar phi[4], phiCandidate[4];
68  label tetLabelCandidate[2], tetPointLabels[2];
69 
70  Type t = pTraits<Type>::zero;
71 
72  // only use face information when the position is on a face
73  if (facei < 0)
74  {
75  const vector& cellCentre = this->pMesh_.cellCentres()[nCell];
76  const labelList& cellFaces = this->pMesh_.cells()[nCell];
77 
78  vector projection = position - cellCentre;
79  tetPoints[3] = cellCentre;
80 
81  // *********************************************************************
82  // project the cell-center through the point onto the face
83  // and get the closest face, ...
84  // *********************************************************************
85 
86  bool foundTet = false;
87  label closestFace = -1;
88  scalar minDistance = GREAT;
89 
90  forAll(cellFaces, facei)
91  {
92  label nFace = cellFaces[facei];
93 
94  vector normal = this->pMeshFaceAreas_[nFace];
95  normal /= mag(normal);
96 
97  const vector& faceCentreTmp = this->pMeshFaceCentres_[nFace];
98 
99  scalar multiplierNumerator = (faceCentreTmp - cellCentre) & normal;
100  scalar multiplierDenominator = projection & normal;
101 
102  // if normal and projection are not orthogonal this could be the one...
103  if (mag(multiplierDenominator) > SMALL)
104  {
105  scalar multiplier = multiplierNumerator/multiplierDenominator;
106  vector iPoint = cellCentre + multiplier*projection;
107  scalar dist = mag(position - iPoint);
108 
109  if (dist < minDistance)
110  {
111  closestFace = nFace;
112  minDistance = dist;
113  }
114  }
115  }
116 
117  // *********************************************************************
118  // find the tetrahedron containing 'position'
119  // from the cell center, face center and
120  // two other points on the face
121  // *********************************************************************
122 
123  minDistance = GREAT;
124  if (closestFace != -1)
125  {
126  label nFace = closestFace;
127  foundTet = findTet
128  (
129  position,
130  nFace,
131  tetPoints,
132  tetLabelCandidate,
133  tetPointLabels,
134  phi,
135  phiCandidate,
136  closestFace,
137  minDistance
138  );
139  }
140 
141  if (!foundTet)
142  {
143  // check if the position is 'just' outside a tet
144  if (minDistance < 1.0e-5)
145  {
146  foundTet = true;
147  for (label i=0; i<4; i++)
148  {
149  phi[i] = phiCandidate[i];
150  }
151  tetPointLabels[0] = tetLabelCandidate[0];
152  tetPointLabels[1] = tetLabelCandidate[1];
153  }
154  }
155 
156  // *********************************************************************
157  // if the search failed check all the cell-faces
158  // *********************************************************************
159 
160  if (!foundTet)
161  {
162  minDistance = GREAT;
163 
164  label facei = 0;
165  while (facei < cellFaces.size() && !foundTet)
166  {
167  label nFace = cellFaces[facei];
168  if (nFace < this->pMeshFaceAreas_.size())
169  {
170  foundTet = findTet
171  (
172  position,
173  nFace,
174  tetPoints,
175  tetLabelCandidate,
176  tetPointLabels,
177  phi,
178  phiCandidate,
179  closestFace,
180  minDistance
181  );
182  }
183  facei++;
184  }
185  }
186 
187  if (!foundTet)
188  {
189  // check if the position is 'just' outside a tet
190  // this time with a more tolerant limit
191  if (minDistance < 1.0e-3)
192  {
193  foundTet = true;
194  for (label i=0; i<4; i++)
195  {
196  phi[i] = phiCandidate[i];
197  }
198  tetPointLabels[0] = tetLabelCandidate[0];
199  tetPointLabels[1] = tetLabelCandidate[1];
200  }
201  }
202 
203  // *********************************************************************
204  // if the tet was found do the interpolation,
205  // otherwise use the closest face value
206  // *********************************************************************
207 
208  if (foundTet)
209  {
210  for (label i=0; i<2; i++)
211  {
212  ts[i] = psip_[tetPointLabels[i]];
213  }
214 
215  if (closestFace < psis_.size())
216  {
217  ts[2] = psis_[closestFace];
218  }
219  else
220  {
221  label patchi = this->pMesh_.boundaryMesh().whichPatch(closestFace);
222 
223  // If the boundary patch is not empty use the face value
224  // else use the cell value
225  if (this->psi_.boundaryField()[patchi].size())
226  {
227  ts[2] = this->psi_.boundaryField()[patchi]
228  [this->pMesh_.boundaryMesh()[patchi].whichFace(closestFace)];
229  }
230  else
231  {
232  ts[2] = this->psi_[nCell];
233  }
234  }
235 
236  ts[3] = this->psi_[nCell];
237 
238  for (label n=0; n<4; n++)
239  {
240  phi[n] = min(1.0, phi[n]);
241  phi[n] = max(0.0, phi[n]);
242 
243  t += phi[n]*ts[n];
244  }
245  }
246  else
247  {
248  Info<< "interpolationCellPointFace<Type>::interpolate"
249  << "(const vector&, const label nCell) const : "
250  << "search failed; using closest cellFace value" << endl
251  << "cell number " << nCell << tab
252  << "position " << position << endl;
253 
254  if (closestFace < psis_.size())
255  {
256  t = psis_[closestFace];
257  }
258  else
259  {
260  label patchi = this->pMesh_.boundaryMesh().whichPatch(closestFace);
261 
262  // If the boundary patch is not empty use the face value
263  // else use the cell value
264  if (this->psi_.boundaryField()[patchi].size())
265  {
266  t = this->psi_.boundaryField()[patchi]
267  [this->pMesh_.boundaryMesh()[patchi].whichFace(closestFace)];
268  }
269  else
270  {
271  t = this->psi_[nCell];
272  }
273  }
274  }
275  }
276  else
277  {
278  bool foundTriangle = findTriangle
279  (
280  position,
281  facei,
282  tetPointLabels,
283  phi
284  );
285 
286  if (foundTriangle)
287  {
288  // add up the point values ...
289  for (label i=0; i<2; i++)
290  {
291  Type vel = psip_[tetPointLabels[i]];
292  t += phi[i]*vel;
293  }
294 
295  // ... and the face value
296  if (facei < psis_.size())
297  {
298  t += phi[2]*psis_[facei];
299  }
300  else
301  {
302  label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
303 
304  // If the boundary patch is not empty use the face value
305  // else use the cell value
306  if (this->psi_.boundaryField()[patchi].size())
307  {
308  t += phi[2]*this->psi_.boundaryField()[patchi]
309  [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
310  }
311  else
312  {
313  t += phi[2]*this->psi_[nCell];
314  }
315  }
316  }
317  else
318  {
319  // use face value only
320  if (facei < psis_.size())
321  {
322  t = psis_[facei];
323  }
324  else
325  {
326  label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
327 
328  // If the boundary patch is not empty use the face value
329  // else use the cell value
330  if (this->psi_.boundaryField()[patchi].size())
331  {
332  t = this->psi_.boundaryField()[patchi]
333  [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
334  }
335  else
336  {
337  t = this->psi_[nCell];
338  }
339  }
340  }
341  }
342 
343  return t;
344 }
345 
346 
347 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
348 
349 } // End namespace Foam
350 
351 // ************************ vim: set sw=4 sts=4 et: ************************ //