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) 2011-2013 OpenFOAM Foundation
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 
27 #include "volFields.H"
28 #include "polyMesh.H"
29 #include "volPointInterpolation.H"
30 #include "linear.H"
31 #include "findCellPointFaceTet.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38 
39 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
40 
41 template<class Type>
43 (
45 )
46 :
48  psip_
49  (
51  (
52  psi,
53  "volPointInterpolate(" + psi.name() + ')',
54  true // use cache
55  )
56  ),
57  psis_(linearInterpolate(psi))
58 {}
59 
60 
61 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
62 
63 template<class Type>
65 (
66  const vector& position,
67  const label cellI,
68  const label faceI
69 ) const
70 {
71  Type ts[4];
72  vector tetPoints[4];
73  scalar phi[4], phiCandidate[4];
74  label tetLabelCandidate[2], tetPointLabels[2];
75 
76  Type t = pTraits<Type>::zero;
77 
78  // only use face information when the position is on a face
79  if (faceI < 0)
80  {
81  const vector& cellCentre = this->pMesh_.cellCentres()[cellI];
82  const labelList& cellFaces = this->pMesh_.cells()[cellI];
83 
84  vector projection = position - cellCentre;
85  tetPoints[3] = cellCentre;
86 
87  // *********************************************************************
88  // project the cell-center through the point onto the face
89  // and get the closest face, ...
90  // *********************************************************************
91 
92  bool foundTet = false;
93  label closestFace = -1;
94  scalar minDistance = GREAT;
95 
96  forAll(cellFaces, faceI)
97  {
98  label nFace = cellFaces[faceI];
99 
100  vector normal = this->pMeshFaceAreas_[nFace];
101  normal /= mag(normal);
102 
103  const vector& faceCentreTmp = this->pMeshFaceCentres_[nFace];
104 
105  scalar multiplierNumerator = (faceCentreTmp - cellCentre) & normal;
106  scalar multiplierDenominator = projection & normal;
107 
108  // if normal and projection are not orthogonal this could
109  // be the one...
110  if (mag(multiplierDenominator) > SMALL)
111  {
112  scalar multiplier = multiplierNumerator/multiplierDenominator;
113  vector iPoint = cellCentre + multiplier*projection;
114  scalar dist = mag(position - iPoint);
115 
116  if (dist < minDistance)
117  {
118  closestFace = nFace;
119  minDistance = dist;
120  }
121  }
122  }
123 
124  // *********************************************************************
125  // find the tetrahedron containing 'position'
126  // from the cell center, face center and
127  // two other points on the face
128  // *********************************************************************
129 
130  minDistance = GREAT;
131  if (closestFace != -1)
132  {
133  label nFace = closestFace;
134  foundTet = findTet
135  (
136  position,
137  nFace,
138  tetPoints,
139  tetLabelCandidate,
140  tetPointLabels,
141  phi,
142  phiCandidate,
143  closestFace,
144  minDistance
145  );
146  }
147 
148  if (!foundTet)
149  {
150  // check if the position is 'just' outside a tet
151  if (minDistance < 1.0e-5)
152  {
153  foundTet = true;
154  for (label i=0; i<4; i++)
155  {
156  phi[i] = phiCandidate[i];
157  }
158  tetPointLabels[0] = tetLabelCandidate[0];
159  tetPointLabels[1] = tetLabelCandidate[1];
160  }
161  }
162 
163  // *********************************************************************
164  // if the search failed check all the cell-faces
165  // *********************************************************************
166 
167  if (!foundTet)
168  {
169  minDistance = GREAT;
170 
171  label faceI = 0;
172  while (faceI < cellFaces.size() && !foundTet)
173  {
174  label nFace = cellFaces[faceI];
175  if (nFace < this->pMeshFaceAreas_.size())
176  {
177  foundTet = findTet
178  (
179  position,
180  nFace,
181  tetPoints,
182  tetLabelCandidate,
183  tetPointLabels,
184  phi,
185  phiCandidate,
186  closestFace,
187  minDistance
188  );
189  }
190  faceI++;
191  }
192  }
193 
194  if (!foundTet)
195  {
196  // check if the position is 'just' outside a tet
197  // this time with a more tolerant limit
198  if (minDistance < 1.0e-3)
199  {
200  foundTet = true;
201  for (label i=0; i<4; i++)
202  {
203  phi[i] = phiCandidate[i];
204  }
205  tetPointLabels[0] = tetLabelCandidate[0];
206  tetPointLabels[1] = tetLabelCandidate[1];
207  }
208  }
209 
210  // *********************************************************************
211  // if the tet was found do the interpolation,
212  // otherwise use the closest face value
213  // *********************************************************************
214 
215  if (foundTet)
216  {
217  for (label i=0; i<2; i++)
218  {
219  ts[i] = psip_[tetPointLabels[i]];
220  }
221 
222  if (closestFace < psis_.size())
223  {
224  ts[2] = psis_[closestFace];
225  }
226  else
227  {
228  label patchI =
229  this->pMesh_.boundaryMesh().whichPatch(closestFace);
230 
231  // If the boundary patch is not empty use the face value
232  // else use the cell value
233  if (this->psi_.boundaryField()[patchI].size())
234  {
235  ts[2] = this->psi_.boundaryField()[patchI]
236  [
237  this->pMesh_.boundaryMesh()[patchI].whichFace
238  (
239  closestFace
240  )
241  ];
242  }
243  else
244  {
245  ts[2] = this->psi_[cellI];
246  }
247  }
248 
249  ts[3] = this->psi_[cellI];
250 
251  for (label n=0; n<4; n++)
252  {
253  phi[n] = min(1.0, phi[n]);
254  phi[n] = max(0.0, phi[n]);
255 
256  t += phi[n]*ts[n];
257  }
258  }
259  else
260  {
261  Info<< "interpolationCellPointFace<Type>::interpolate"
262  << "(const vector&, const label cellI) const : "
263  << "search failed; using closest cellFace value" << endl
264  << "cell number " << cellI << tab
265  << "position " << position << endl;
266 
267  if (closestFace < psis_.size())
268  {
269  t = psis_[closestFace];
270  }
271  else
272  {
273  label patchI =
274  this->pMesh_.boundaryMesh().whichPatch(closestFace);
275 
276  // If the boundary patch is not empty use the face value
277  // else use the cell value
278  if (this->psi_.boundaryField()[patchI].size())
279  {
280  t = this->psi_.boundaryField()[patchI]
281  [
282  this->pMesh_.boundaryMesh()[patchI].whichFace
283  (
284  closestFace
285  )
286  ];
287  }
288  else
289  {
290  t = this->psi_[cellI];
291  }
292  }
293  }
294  }
295  else
296  {
297  bool foundTriangle = findTriangle
298  (
299  position,
300  faceI,
301  tetPointLabels,
302  phi
303  );
304 
305  if (foundTriangle)
306  {
307  // add up the point values ...
308  for (label i=0; i<2; i++)
309  {
310  Type vel = psip_[tetPointLabels[i]];
311  t += phi[i]*vel;
312  }
313 
314  // ... and the face value
315  if (faceI < psis_.size())
316  {
317  t += phi[2]*psis_[faceI];
318  }
319  else
320  {
321  label patchI = this->pMesh_.boundaryMesh().whichPatch(faceI);
322 
323  // If the boundary patch is not empty use the face value
324  // else use the cell value
325  if (this->psi_.boundaryField()[patchI].size())
326  {
327  t += phi[2]*this->psi_.boundaryField()[patchI]
328  [this->pMesh_.boundaryMesh()[patchI].whichFace(faceI)];
329  }
330  else
331  {
332  t += phi[2]*this->psi_[cellI];
333  }
334  }
335  }
336  else
337  {
338  // use face value only
339  if (faceI < psis_.size())
340  {
341  t = psis_[faceI];
342  }
343  else
344  {
345  label patchI = this->pMesh_.boundaryMesh().whichPatch(faceI);
346 
347  // If the boundary patch is not empty use the face value
348  // else use the cell value
349  if (this->psi_.boundaryField()[patchI].size())
350  {
351  t = this->psi_.boundaryField()[patchI]
352  [this->pMesh_.boundaryMesh()[patchI].whichFace(faceI)];
353  }
354  else
355  {
356  t = this->psi_[cellI];
357  }
358  }
359  }
360  }
361 
362  return t;
363 }
364 
365 
366 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
367 
368 } // End namespace Foam
369 
370 // ************************************************************************* //
dimensioned< scalar > mag(const dimensioned< Type > &)
static const volPointInterpolation & New(const fvMesh &mesh)
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
Definition: curveTools.C:75
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Type interpolate(const vector &position, const label cellI, const label faceI=-1) const
Interpolate field to the given point in the given cell.
Abstract base class for interpolation.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
messageStream Info
const Mesh & mesh() const
Return mesh.
find the tetrahedron in which the position is. while searching for the tet, store the tet closest to ...
Namespace for OpenFOAM.
label n
const double e
Elementary charge.
Definition: doubleFloat.H:78
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
const volScalarField & psi
Definition: createFields.H:24
const word & name() const
Return name.
Definition: IOobject.H:260
surfaceScalarField & phi
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > linearInterpolate(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: linear.H:107
Generic GeometricField class.
static const char tab
Definition: Ostream.H:259
Traits class for primitives.
Definition: pTraits.H:50
A normal distribution model.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Tet storage. Null constructable (unfortunately tetrahedron<point, point> is not)
Definition: tetPoints.H:50
interpolationCellPointFace(const GeometricField< Type, fvPatchField, volMesh > &psi)
Construct from components.