particleI.H
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-2018, 2020 OpenFOAM Foundation
9  Copyright (C) 2011-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "polyMesh.H"
30 #include "Time.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 void Foam::particle::stationaryTetGeometry
35 (
36  vector& centre,
37  vector& base,
38  vector& vertex1,
39  vector& vertex2
40 ) const
41 {
42  const triFace triIs(currentTetIndices().faceTriIs(mesh_));
43  const vectorField& ccs = mesh_.cellCentres();
44  const pointField& pts = mesh_.points();
45 
46  centre = ccs[celli_];
47  base = pts[triIs[0]];
48  vertex1 = pts[triIs[1]];
49  vertex2 = pts[triIs[2]];
50 }
51 
52 
53 inline Foam::barycentricTensor Foam::particle::stationaryTetTransform() const
54 {
55  vector centre, base, vertex1, vertex2;
56  stationaryTetGeometry(centre, base, vertex1, vertex2);
57 
58  return barycentricTensor(centre, base, vertex1, vertex2);
59 }
60 
61 
62 inline void Foam::particle::movingTetGeometry
63 (
64  const scalar fraction,
65  Pair<vector>& centre,
66  Pair<vector>& base,
67  Pair<vector>& vertex1,
68  Pair<vector>& vertex2
69 ) const
70 {
71  const triFace triIs(currentTetIndices().faceTriIs(mesh_));
72 
73  const pointField& ptsOld = mesh_.oldPoints();
74  const pointField& ptsNew = mesh_.points();
75 
76  // !!! <-- We would be better off using mesh_.cellCentres() here. However,
77  // we need to put a mesh_.oldCellCentres() method in for this to work. The
78  // values obtained from the mesh and those obtained from the cell do not
79  // necessarily match. See mantis #1993.
80  //const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
81  //const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
82 
83  const vector ccOld = mesh_.oldCellCentres()[celli_];
84  const vector ccNew = mesh_.cellCentres()[celli_];
85 
86  // Old and new points and cell centres are not sub-cycled. If we are sub-
87  // cycling, then we have to account for the timestep change here by
88  // modifying the fractions that we take of the old and new geometry.
89  const Pair<scalar> s = stepFractionSpan();
90  const scalar f0 = s[0] + stepFraction_*s[1], f1 = fraction*s[1];
91 
92  centre[0] = ccOld + f0*(ccNew - ccOld);
93  base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
94  vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
95  vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
96 
97  centre[1] = f1*(ccNew - ccOld);
98  base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
99  vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
100  vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
101 }
102 
103 
104 inline Foam::Pair<Foam::barycentricTensor> Foam::particle::movingTetTransform
105 (
106  const scalar fraction
107 ) const
108 {
109  Pair<vector> centre, base, vertex1, vertex2;
110  movingTetGeometry(fraction, centre, base, vertex1, vertex2);
111 
112  return
113  Pair<barycentricTensor>
114  (
115  barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]),
116  barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1])
117  );
118 }
119 
120 
121 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
122 
123 inline Foam::label Foam::particle::getNewParticleID() const
124 {
125  label id = particleCount_++;
126 
127  if (id == labelMax)
128  {
130  << "Particle counter has overflowed. This might cause problems"
131  << " when reconstructing particle tracks." << endl;
132  }
133  return id;
134 }
135 
136 
138 {
139  return mesh_;
140 }
141 
142 
144 {
145  return coordinates_;
146 }
147 
148 
149 inline Foam::label Foam::particle::cell() const
150 {
151  return celli_;
152 }
153 
154 
155 inline Foam::label& Foam::particle::cell()
156 {
157  return celli_;
158 }
159 
160 
161 inline Foam::label Foam::particle::tetFace() const
162 {
163  return tetFacei_;
164 }
165 
166 
167 inline Foam::label& Foam::particle::tetFace()
168 {
169  return tetFacei_;
170 }
171 
172 
173 inline Foam::label Foam::particle::tetPt() const
174 {
175  return tetPti_;
176 }
177 
178 
179 inline Foam::label& Foam::particle::tetPt()
180 {
181  return tetPti_;
182 }
183 
184 
185 inline Foam::label Foam::particle::face() const
186 {
187  return facei_;
188 }
189 
190 
191 inline Foam::label& Foam::particle::face()
192 {
193  return facei_;
194 }
195 
196 
197 inline Foam::scalar Foam::particle::stepFraction() const
198 {
199  return stepFraction_;
200 }
201 
202 
203 inline Foam::scalar& Foam::particle::stepFraction()
204 {
205  return stepFraction_;
206 }
207 
208 
209 inline Foam::label Foam::particle::origProc() const
210 {
211  return origProc_;
212 }
213 
214 
215 inline Foam::label& Foam::particle::origProc()
216 {
217  return origProc_;
218 }
219 
220 
221 inline Foam::label Foam::particle::origId() const
222 {
223  return origId_;
224 }
225 
226 
227 inline Foam::label& Foam::particle::origId()
228 {
229  return origId_;
230 }
231 
232 
234 {
235  if (mesh_.time().subCycling())
236  {
237  const TimeState& tsNew = mesh_.time();
238  const TimeState& tsOld = mesh_.time().prevTimeState();
239 
240  const scalar tFrac =
241  (
242  (tsNew.value() - tsNew.deltaTValue())
243  - (tsOld.value() - tsOld.deltaTValue())
244  )/tsOld.deltaTValue();
245 
246  const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue();
247 
248  return Pair<scalar>(tFrac, dtFrac);
249  }
250  else
251  {
252  return Pair<scalar>(0, 1);
253  }
254 }
255 
256 
257 inline Foam::scalar Foam::particle::currentTimeFraction() const
258 {
259  const Pair<scalar> s = stepFractionSpan();
260 
261  return s[0] + stepFraction_*s[1];
262 }
263 
264 
266 {
267  return tetIndices(celli_, tetFacei_, tetPti_);
268 }
269 
270 
272 {
273  if (mesh_.moving() && stepFraction_ != 1)
274  {
275  return movingTetTransform(0)[0];
276  }
277  else
278  {
279  return stationaryTetTransform();
280  }
281 }
282 
283 
285 {
286  return currentTetIndices().faceTri(mesh_).unitNormal();
287 }
288 
289 
290 inline bool Foam::particle::onFace() const
291 {
292  return facei_ >= 0;
293 }
294 
295 
297 {
298  return onFace() && mesh_.isInternalFace(facei_);
299 }
300 
301 
303 {
304  return onFace() && !mesh_.isInternalFace(facei_);
305 }
306 
307 
308 inline Foam::label Foam::particle::patch() const
309 {
310  return onFace() ? mesh_.boundaryMesh().whichPatch(facei_) : -1;
311 }
312 
313 
315 {
316  return currentTetTransform() & coordinates_;
317 }
318 
319 
321 {
322  stepFraction_ = 0;
323  nBehind_ = 0;
324  behind_ = 0;
325 }
326 
327 
329 {
330  if (!onBoundaryFace())
331  {
333  << "Patch data was requested for a particle that isn't on a patch"
334  << exit(FatalError);
335  }
336 
337  if ((mesh_.moving() && stepFraction_ != 1))
338  {
339  Pair<vector> centre, base, vertex1, vertex2;
340  movingTetGeometry(1, centre, base, vertex1, vertex2);
341 
342  n = triPointRef(base[0], vertex1[0], vertex2[0]).unitNormal();
343 
344  // Interpolate the motion of the three face vertices to the current
345  // coordinates
346  U =
347  coordinates_.b()*base[1]
348  + coordinates_.c()*vertex1[1]
349  + coordinates_.d()*vertex2[1];
350 
351  // The movingTetGeometry method gives the motion as a displacement
352  // across the time-step, so we divide by the time-step to get velocity
353  U /= mesh_.time().deltaTValue();
354  }
355  else
356  {
357  vector centre, base, vertex1, vertex2;
358  stationaryTetGeometry(centre, base, vertex1, vertex2);
359 
360  n = triPointRef(base, vertex1, vertex2).unitNormal();
361 
362  U = Zero;
363  }
364 }
365 
366 
367 // ************************************************************************* //
Foam::particle::reset
void reset()
Reset particle data.
Definition: particleI.H:320
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1069
Foam::labelMax
constexpr label labelMax
Definition: label.H:61
Foam::barycentricTensor
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
Definition: barycentricTensor.H:48
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::particle::onFace
bool onFace() const
Is the particle on a face?
Definition: particleI.H:290
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::particle::origId
label origId() const
Return the particle ID on the originating processor.
Definition: particleI.H:221
Foam::particle::stepFraction
scalar stepFraction() const
Return the fraction of time-step completed.
Definition: particleI.H:197
Foam::particle::tetFace
label tetFace() const
Return current tet face particle is in.
Definition: particleI.H:161
Foam::particle::getNewParticleID
label getNewParticleID() const
Get unique particle creation id.
Definition: particleI.H:123
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::particle::currentTetIndices
tetIndices currentTetIndices() const
Return the indices of the current tet that the.
Definition: particleI.H:265
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:434
polyMesh.H
Foam::particle::onInternalFace
bool onInternalFace() const
Is the particle on an internal face?
Definition: particleI.H:296
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
Foam::TimeState::deltaTValue
scalar deltaTValue() const noexcept
Return time step value.
Definition: TimeStateI.H:43
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::particle::coordinates
const barycentric & coordinates() const
Return current particle coordinates.
Definition: particleI.H:143
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::TimeState
The time value with time-stepping information, user-defined remapping, etc.
Definition: TimeState.H:51
Foam::particle::position
vector position() const
Return current particle position.
Definition: particleI.H:314
Foam::Barycentric< scalar >
Foam::particle::patch
label patch() const
Return the index of patch that the particle is on.
Definition: particleI.H:308
Foam::particle::mesh
const polyMesh & mesh() const
Return the mesh database.
Definition: particleI.H:137
Foam::particle::currentTimeFraction
scalar currentTimeFraction() const
Return the current fraction within the timestep. This differs.
Definition: particleI.H:257
Foam::particle::origProc
label origProc() const
Return the originating processor ID.
Definition: particleI.H:209
Foam::FatalError
error FatalError
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::particle::onBoundaryFace
bool onBoundaryFace() const
Is the particle on a boundary face?
Definition: particleI.H:302
Time.H
Foam::tetIndices
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:83
U
U
Definition: pEqn.H:72
Foam::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:84
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::Pair
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:54
Foam::Vector< scalar >
Foam::particle::stepFractionSpan
Pair< scalar > stepFractionSpan() const
Return the step fraction change within the overall time-step.
Definition: particleI.H:233
Foam::BarycentricTensor
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can represent a barycentric transfo...
Definition: BarycentricTensor.H:57
Foam::particle::tetPt
label tetPt() const
Return current tet face particle is in.
Definition: particleI.H:173
Foam::particle::normal
vector normal() const
The (unit) normal of the tri on tetFacei_ for the current tet.
Definition: particleI.H:284
Foam::particle::cell
label cell() const
Return current cell particle is in.
Definition: particleI.H:149
Foam::particle::face
label face() const
Return current face particle is on otherwise -1.
Definition: particleI.H:185
triFace
face triFace(3)
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::triPointRef
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
Foam::particle::patchData
void patchData(vector &n, vector &U) const
Get the normal and velocity of the current patch location.
Definition: particleI.H:328
Foam::particle::currentTetTransform
barycentricTensor currentTetTransform() const
Return the current tet transformation tensor.
Definition: particleI.H:271