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