particleI.H
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-2019 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 
26 #include "polyMesh.H"
27 #include "Time.H"
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
31 void Foam::particle::stationaryTetGeometry
32 (
33  vector& centre,
34  vector& base,
35  vector& vertex1,
36  vector& vertex2
37 ) const
38 {
39  const triFace triIs(currentTetIndices().faceTriIs(mesh_));
40  const vectorField& ccs = mesh_.cellCentres();
41  const pointField& pts = mesh_.points();
42 
43  centre = ccs[celli_];
44  base = pts[triIs[0]];
45  vertex1 = pts[triIs[1]];
46  vertex2 = pts[triIs[2]];
47 }
48 
49 
50 inline Foam::barycentricTensor Foam::particle::stationaryTetTransform() const
51 {
52  vector centre, base, vertex1, vertex2;
53  stationaryTetGeometry(centre, base, vertex1, vertex2);
54 
55  return barycentricTensor(centre, base, vertex1, vertex2);
56 }
57 
58 
59 inline void Foam::particle::movingTetGeometry
60 (
61  const scalar fraction,
62  Pair<vector>& centre,
63  Pair<vector>& base,
64  Pair<vector>& vertex1,
65  Pair<vector>& vertex2
66 ) const
67 {
68  const triFace triIs(currentTetIndices().faceTriIs(mesh_));
69  const pointField& ptsOld = mesh_.oldPoints();
70  const pointField& ptsNew = mesh_.points();
71  const vector ccOld = mesh_.oldCellCentres()[celli_];
72  const vector ccNew = mesh_.cellCentres()[celli_];
73 
74  // Old and new points and cell centres are not sub-cycled. If we are sub-
75  // cycling, then we have to account for the timestep change here by
76  // modifying the fractions that we take of the old and new geometry.
77  const Pair<scalar> s = stepFractionSpan();
78  const scalar f0 = s[0] + stepFraction_*s[1], f1 = fraction*s[1];
79 
80  centre[0] = ccOld + f0*(ccNew - ccOld);
81  base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
82  vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
83  vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
84 
85  centre[1] = f1*(ccNew - ccOld);
86  base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
87  vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
88  vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
89 }
90 
91 
92 inline Foam::Pair<Foam::barycentricTensor> Foam::particle::movingTetTransform
93 (
94  const scalar fraction
95 ) const
96 {
97  Pair<vector> centre, base, vertex1, vertex2;
98  movingTetGeometry(fraction, centre, base, vertex1, vertex2);
99 
100  return
101  Pair<barycentricTensor>
102  (
103  barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]),
104  barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1])
105  );
106 }
107 
108 
109 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
110 
112 {
113  label id = particleCount_++;
114 
115  if (id == labelMax)
116  {
118  << "Particle counter has overflowed. This might cause problems"
119  << " when reconstructing particle tracks." << endl;
120  }
121  return id;
122 }
123 
124 
126 {
127  return mesh_;
128 }
129 
130 
132 {
133  return coordinates_;
134 }
135 
136 
138 {
139  return celli_;
140 }
141 
142 
144 {
145  return tetFacei_;
146 }
147 
148 
150 {
151  return tetPti_;
152 }
153 
154 
156 {
157  return facei_;
158 }
159 
160 
161 inline Foam::scalar Foam::particle::stepFraction() const
162 {
163  return stepFraction_;
164 }
165 
166 
167 inline Foam::scalar& Foam::particle::stepFraction()
168 {
169  return stepFraction_;
170 }
171 
172 
174 {
175  return origProc_;
176 }
177 
178 
180 {
181  return origProc_;
182 }
183 
184 
186 {
187  return origId_;
188 }
189 
190 
192 {
193  return origId_;
194 }
195 
196 
198 {
199  if (mesh_.time().subCycling())
200  {
201  const TimeState& tsNew = mesh_.time();
202  const TimeState& tsOld = mesh_.time().prevTimeState();
203 
204  const scalar tFrac =
205  (
206  (tsNew.value() - tsNew.deltaTValue())
207  - (tsOld.value() - tsOld.deltaTValue())
208  )/tsOld.deltaTValue();
209 
210  const scalar dtFrac = tsNew.deltaTValue()/tsOld.deltaTValue();
211 
212  return Pair<scalar>(tFrac, dtFrac);
213  }
214  else
215  {
216  return Pair<scalar>(0, 1);
217  }
218 }
219 
220 
221 inline Foam::scalar Foam::particle::currentTimeFraction() const
222 {
223  const Pair<scalar> s = stepFractionSpan();
224 
225  return s[0] + stepFraction_*s[1];
226 }
227 
228 
230 {
231  return tetIndices(celli_, tetFacei_, tetPti_);
232 }
233 
234 
236 {
237  if (mesh_.moving() && stepFraction_ != 1)
238  {
239  return movingTetTransform(0)[0];
240  }
241  else
242  {
243  return stationaryTetTransform();
244  }
245 }
246 
247 
249 {
250  return currentTetIndices().faceTri(mesh_).normal();
251 }
252 
253 
254 inline bool Foam::particle::onFace() const
255 {
256  return facei_ >= 0;
257 }
258 
259 
261 {
262  return onFace() && mesh_.isInternalFace(facei_);
263 }
264 
265 
267 {
268  return onFace() && !mesh_.isInternalFace(facei_);
269 }
270 
271 
273 {
274  return onFace() ? mesh_.boundaryMesh().whichPatch(facei_) : -1;
275 }
276 
277 
279 {
280  return currentTetTransform() & coordinates_;
281 }
282 
283 
285 {
286  stepFraction_ = 0;
287  nBehind_ = 0;
288  behind_ = 0;
289 }
290 
291 
292 void Foam::particle::patchData(vector& normal, vector& displacement) const
293 {
294  if (!onBoundaryFace())
295  {
297  << "Patch data was requested for a particle that isn't on a patch"
298  << exit(FatalError);
299  }
300 
301  if (mesh_.moving() && stepFraction_ != 1)
302  {
303  Pair<vector> centre, base, vertex1, vertex2;
304  movingTetGeometry(1, centre, base, vertex1, vertex2);
305 
306  normal = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
307 
308  // Interpolate the motion of the three face vertices to the current
309  // coordinates
310  displacement =
311  coordinates_.b()*base[1]
312  + coordinates_.c()*vertex1[1]
313  + coordinates_.d()*vertex2[1];
314  }
315  else
316  {
317  vector centre, base, vertex1, vertex2;
318  stationaryTetGeometry(centre, base, vertex1, vertex2);
319 
320  normal = triPointRef(base, vertex1, vertex2).normal();
321 
322  displacement = Zero;
323  }
324 }
325 
326 
327 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
label origProc() const
Return the originating processor ID.
Definition: particleI.H:173
const polyMesh & mesh() const
Return the mesh database.
Definition: particleI.H:125
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
bool moving() const
Is mesh moving.
Definition: polyMesh.H:512
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
The time value with time-stepping information, user-defined remapping, etc.
Definition: TimeState.H:48
bool onBoundaryFace() const
Is the particle on a boundary face?
Definition: particleI.H:266
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the face for.
Definition: tetIndicesI.H:113
const Cmpt & b() const
Definition: BarycentricI.H:66
bool subCycling() const
Return true if time currently being sub-cycled, otherwise false.
Definition: Time.H:428
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
Pair< scalar > stepFractionSpan() const
Return the step fraction change within the overall time-step.
Definition: particleI.H:197
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const TimeState & prevTimeState() const
Return previous TimeState if time is being sub-cycled.
Definition: Time.H:434
scalar f1
Definition: createFields.H:28
label tetPt() const
Return current tet face particle is in.
Definition: particleI.H:149
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
vector normal() const
Return the normal of the tri on tetFacei_ for the.
Definition: particleI.H:248
void patchData(vector &normal, vector &displacement) const
Get the normal and displacement of the current patch location.
Definition: particleI.H:292
label cell() const
Return current cell particle is in.
Definition: particleI.H:137
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
label face() const
Return current face particle is on otherwise -1.
Definition: particleI.H:155
void reset()
Set step fraction and behind data to zero in preparation for a new.
Definition: particleI.H:284
face triFace(3)
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:41
virtual const pointField & oldPoints() const
Return old points for mesh motion.
Definition: polyMesh.C:1181
const Type & value() const
Return const reference to value.
const Cmpt & d() const
Definition: BarycentricI.H:80
static const label labelMax
Definition: label.H:62
static const zero Zero
Definition: zero.H:97
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
barycentricTensor currentTetTransform() const
Return the current tet transformation tensor.
Definition: particleI.H:235
const vectorField & cellCentres() const
vector normal() const
Return unit normal.
Definition: triangleI.H:110
label origId() const
Return the particle ID on the originating processor.
Definition: particleI.H:185
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const Cmpt & c() const
Definition: BarycentricI.H:73
bool onFace() const
Is the particle on a face?
Definition: particleI.H:254
label tetFace() const
Return current tet face particle is in.
Definition: particleI.H:143
label getNewParticleID() const
Get unique particle creation id.
Definition: particleI.H:111
const Time & time() const
Return time.
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
scalar stepFraction() const
Return the fraction of time-step completed.
Definition: particleI.H:161
tetIndices currentTetIndices() const
Return the indices of the current tet that the.
Definition: particleI.H:229
virtual const pointField & oldCellCentres() const
Return old points for mesh motion.
Definition: polyMesh.C:1199
#define WarningInFunction
Report a warning using Foam::Warning.
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can represent a barycentric transfo...
const barycentric & coordinates() const
Return current particle coordinates.
Definition: particleI.H:131
static label particleCount_
Cumulative particle counter - used to provide unique ID.
Definition: particle.H:368
Field< vector > vectorField
Specialisation of Field<T> for vector.
bool onInternalFace() const
Is the particle on an internal face?
Definition: particleI.H:260
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
vector position() const
Return current particle position.
Definition: particleI.H:278
label patch() const
Return the index of patch that the particle is on.
Definition: particleI.H:272
scalar currentTimeFraction() const
Return the current fraction within the timestep. This differs.
Definition: particleI.H:221