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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
29 #include "volFields.H"
30 #include "polyMesh.H"
31 #include "volPointInterpolation.H"
32 #include "linear.H"
33 #include "findCellPointFaceTet.H"
35 
36 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
37 
38 template<class Type>
40 (
42 )
43 :
45  psip_
46  (
48  (
49  psi,
50  "volPointInterpolate(" + psi.name() + ')',
51  true // use cache
52  )
53  ),
54  psis_(linearInterpolate(psi))
55 {}
56 
57 
58 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
59 
60 template<class Type>
62 (
63  const vector& position,
64  const label celli,
65  const label facei
66 ) const
67 {
68  Type ts[4];
69  vector tetPoints[4];
70  scalar phi[4], phiCandidate[4];
71  label tetLabelCandidate[2], tetPointLabels[2];
72 
73  Type t = Zero;
74 
75  // only use face information when the position is on a face
76  if (facei < 0)
77  {
78  const vector& cellCentre = this->pMesh_.cellCentres()[celli];
79  const labelList& cellFaces = this->pMesh_.cells()[celli];
80 
81  vector projection = position - cellCentre;
82  tetPoints[3] = cellCentre;
83 
84  // *********************************************************************
85  // project the cell-center through the point onto the face
86  // and get the closest face, ...
87  // *********************************************************************
88 
89  bool foundTet = false;
90  label closestFace = -1;
91  scalar minDistance = GREAT;
92 
93  for (const label nFace : cellFaces)
94  {
95  const vector normal = normalised(this->pMeshFaceAreas_[nFace]);
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
103  // be the one...
104  if (mag(multiplierDenominator) > SMALL)
105  {
106  scalar multiplier = multiplierNumerator/multiplierDenominator;
107  vector iPoint = cellCentre + multiplier*projection;
108  scalar dist = mag(position - iPoint);
109 
110  if (dist < minDistance)
111  {
112  closestFace = nFace;
113  minDistance = dist;
114  }
115  }
116  }
117 
118  // *********************************************************************
119  // find the tetrahedron containing 'position'
120  // from the cell center, face center and
121  // two other points on the face
122  // *********************************************************************
123 
124  minDistance = GREAT;
125  if (closestFace != -1)
126  {
127  label nFace = closestFace;
128  foundTet = findTet
129  (
130  position,
131  nFace,
132  tetPoints,
133  tetLabelCandidate,
134  tetPointLabels,
135  phi,
136  phiCandidate,
137  closestFace,
138  minDistance
139  );
140  }
141 
142  if (!foundTet)
143  {
144  // check if the position is 'just' outside a tet
145  if (minDistance < 1.0e-5)
146  {
147  foundTet = true;
148  for (label i=0; i<4; i++)
149  {
150  phi[i] = phiCandidate[i];
151  }
152  tetPointLabels[0] = tetLabelCandidate[0];
153  tetPointLabels[1] = tetLabelCandidate[1];
154  }
155  }
156 
157  // *********************************************************************
158  // if the search failed check all the cell-faces
159  // *********************************************************************
160 
161  if (!foundTet)
162  {
163  minDistance = GREAT;
164 
165  label facei = 0;
166  while (facei < cellFaces.size() && !foundTet)
167  {
168  label nFace = cellFaces[facei];
169  if (nFace < this->pMeshFaceAreas_.size())
170  {
171  foundTet = findTet
172  (
173  position,
174  nFace,
175  tetPoints,
176  tetLabelCandidate,
177  tetPointLabels,
178  phi,
179  phiCandidate,
180  closestFace,
181  minDistance
182  );
183  }
184  facei++;
185  }
186  }
187 
188  if (!foundTet)
189  {
190  // check if the position is 'just' outside a tet
191  // this time with a more tolerant limit
192  if (minDistance < 1.0e-3)
193  {
194  foundTet = true;
195  for (label i=0; i<4; i++)
196  {
197  phi[i] = phiCandidate[i];
198  }
199  tetPointLabels[0] = tetLabelCandidate[0];
200  tetPointLabels[1] = tetLabelCandidate[1];
201  }
202  }
203 
204  // *********************************************************************
205  // if the tet was found do the interpolation,
206  // otherwise use the closest face value
207  // *********************************************************************
208 
209  if (foundTet)
210  {
211  for (label i=0; i<2; i++)
212  {
213  ts[i] = psip_[tetPointLabels[i]];
214  }
215 
216  if (closestFace < psis_.size())
217  {
218  ts[2] = psis_[closestFace];
219  }
220  else
221  {
222  label patchi =
223  this->pMesh_.boundaryMesh().whichPatch(closestFace);
224 
225  // If the boundary patch is not empty use the face value
226  // else use the cell value
227  if (this->psi_.boundaryField()[patchi].size())
228  {
229  ts[2] = this->psi_.boundaryField()[patchi]
230  [
231  this->pMesh_.boundaryMesh()[patchi].whichFace
232  (
233  closestFace
234  )
235  ];
236  }
237  else
238  {
239  ts[2] = this->psi_[celli];
240  }
241  }
242 
243  ts[3] = this->psi_[celli];
244 
245  for (label n=0; n<4; n++)
246  {
247  phi[n] = min(1.0, phi[n]);
248  phi[n] = max(0.0, phi[n]);
249 
250  t += phi[n]*ts[n];
251  }
252  }
253  else
254  {
256  << "search failed; using closest cellFace value" << endl
257  << "cell number " << celli << tab
258  << "position " << position << endl;
259 
260  if (closestFace < psis_.size())
261  {
262  t = psis_[closestFace];
263  }
264  else
265  {
266  label patchi =
267  this->pMesh_.boundaryMesh().whichPatch(closestFace);
268 
269  // If the boundary patch is not empty use the face value
270  // else use the cell value
271  if (this->psi_.boundaryField()[patchi].size())
272  {
273  t = this->psi_.boundaryField()[patchi]
274  [
275  this->pMesh_.boundaryMesh()[patchi].whichFace
276  (
277  closestFace
278  )
279  ];
280  }
281  else
282  {
283  t = this->psi_[celli];
284  }
285  }
286  }
287  }
288  else
289  {
290  bool foundTriangle = findTriangle
291  (
292  position,
293  facei,
294  tetPointLabels,
295  phi
296  );
297 
298  if (foundTriangle)
299  {
300  // add up the point values ...
301  for (label i=0; i<2; i++)
302  {
303  Type vel = psip_[tetPointLabels[i]];
304  t += phi[i]*vel;
305  }
306 
307  // ... and the face value
308  if (facei < psis_.size())
309  {
310  t += phi[2]*psis_[facei];
311  }
312  else
313  {
314  label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
315 
316  // If the boundary patch is not empty use the face value
317  // else use the cell value
318  if (this->psi_.boundaryField()[patchi].size())
319  {
320  t += phi[2]*this->psi_.boundaryField()[patchi]
321  [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
322  }
323  else
324  {
325  t += phi[2]*this->psi_[celli];
326  }
327  }
328  }
329  else
330  {
331  // use face value only
332  if (facei < psis_.size())
333  {
334  t = psis_[facei];
335  }
336  else
337  {
338  label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
339 
340  // If the boundary patch is not empty use the face value
341  // else use the cell value
342  if (this->psi_.boundaryField()[patchi].size())
343  {
344  t = this->psi_.boundaryField()[patchi]
345  [this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
346  }
347  else
348  {
349  t = this->psi_[celli];
350  }
351  }
352  }
353  }
354 
355  return t;
356 }
357 
358 
359 // ************************************************************************* //
volFields.H
InfoInFunction
#define InfoInFunction
Report an information message using Foam::Info.
Definition: messageStream.H:350
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
polyMesh.H
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
n
label n
Definition: TABSMDCalcMethod2.H:31
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::interpolation
Abstract base class for interpolation.
Definition: mappedPatchFieldBase.H:96
Foam::tetPoints
Tet storage. Null constructable (unfortunately tetrahedron<point, point> is not)
Definition: tetPoints.H:53
Foam::normalised
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:487
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
volPointInterpolation.H
interpolationCellPointFace.H
Foam::tab
constexpr char tab
Definition: Ostream.H:403
Foam::linearInterpolate
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > linearInterpolate(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: linear.H:112
Foam::Vector< scalar >
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::interpolationCellPointFace::interpolate
Type interpolate(const vector &position, const label celli, const label facei=-1) const
Interpolate field to the given point in the given cell.
Definition: interpolationCellPointFace.C:62
findCellPointFaceTriangle.H
Foam::interpolationCellPointFace::interpolationCellPointFace
interpolationCellPointFace(const GeometricField< Type, fvPatchField, volMesh > &psi)
Construct from components.
Definition: interpolationCellPointFace.C:40
Foam::GeometricField< Type, fvPatchField, volMesh >
psi
const volScalarField & psi
Definition: createFieldRefs.H:1
Foam::fac::interpolate
static tmp< GeometricField< Type, faePatchField, edgeMesh > > interpolate(const GeometricField< Type, faPatchField, areaMesh > &tvf, const edgeScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
findCellPointFaceTet.H
find the tetrahedron in which the position is. while searching for the tet, store the tet closest to ...