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