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-2023 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 "particle.H"
27 #include "polyMesh.H"
28 #include "Time.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 void Foam::particle::stationaryTetGeometry
33 (
34  const polyMesh& mesh,
35  vector& centre,
36  vector& base,
37  vector& vertex1,
38  vector& vertex2
39 ) const
40 {
41  const triFace triIs(currentTetIndices(mesh).faceTriIs(mesh));
42  const vectorField& ccs = mesh.cellCentres();
43  const pointField& pts = mesh.points();
44 
45  centre = ccs[celli_];
46  base = pts[triIs[0]];
47  vertex1 = pts[triIs[1]];
48  vertex2 = pts[triIs[2]];
49 }
50 
51 
52 inline Foam::barycentricTensor Foam::particle::stationaryTetTransform
53 (
54  const polyMesh& mesh
55 ) const
56 {
57  vector centre, base, vertex1, vertex2;
58  stationaryTetGeometry(mesh, centre, base, vertex1, vertex2);
59 
60  return barycentricTensor(centre, base, vertex1, vertex2);
61 }
62 
63 
64 inline void Foam::particle::movingTetGeometry
65 (
66  const polyMesh& mesh,
67  const scalar fraction,
68  Pair<vector>& centre,
69  Pair<vector>& base,
70  Pair<vector>& vertex1,
71  Pair<vector>& vertex2
72 ) const
73 {
74  const triFace triIs(currentTetIndices(mesh).faceTriIs(mesh));
75  const pointField& ptsOld = mesh.oldPoints();
76  const pointField& ptsNew = mesh.points();
77  const vector ccOld = mesh.oldCellCentres()[celli_];
78  const vector ccNew = mesh.cellCentres()[celli_];
79 
80  const scalar f0 = stepFraction_, f1 = fraction;
81 
82  centre[0] = ccOld + f0*(ccNew - ccOld);
83  base[0] = ptsOld[triIs[0]] + f0*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
84  vertex1[0] = ptsOld[triIs[1]] + f0*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
85  vertex2[0] = ptsOld[triIs[2]] + f0*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
86 
87  centre[1] = f1*(ccNew - ccOld);
88  base[1] = f1*(ptsNew[triIs[0]] - ptsOld[triIs[0]]);
89  vertex1[1] = f1*(ptsNew[triIs[1]] - ptsOld[triIs[1]]);
90  vertex2[1] = f1*(ptsNew[triIs[2]] - ptsOld[triIs[2]]);
91 }
92 
93 
94 inline Foam::Pair<Foam::barycentricTensor> Foam::particle::movingTetTransform
95 (
96  const polyMesh& mesh,
97  const scalar fraction
98 ) const
99 {
100  Pair<vector> centre, base, vertex1, vertex2;
101  movingTetGeometry(mesh, fraction, centre, base, vertex1, vertex2);
102 
103  return
104  Pair<barycentricTensor>
105  (
106  barycentricTensor(centre[0], base[0], vertex1[0], vertex2[0]),
107  barycentricTensor(centre[1], base[1], vertex1[1], vertex2[1])
108  );
109 }
110 
111 
112 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
113 
115 {
116  label id = particleCount_++;
117 
118  if (id == labelMax)
119  {
121  << "Particle counter has overflowed. This might cause problems"
122  << " when reconstructing particle tracks." << endl;
123  }
124  return id;
125 }
126 
127 
129 {
130  return coordinates_;
131 }
132 
133 
135 {
136  return celli_;
137 }
138 
139 
141 {
142  return tetFacei_;
143 }
144 
145 
147 {
148  return tetPti_;
149 }
150 
151 
153 {
154  return facei_;
155 }
156 
157 
159 {
160  return facei_;
161 }
162 
163 
164 inline Foam::scalar Foam::particle::stepFraction() const
165 {
166  return stepFraction_;
167 }
168 
169 
170 inline Foam::scalar& Foam::particle::stepFraction()
171 {
172  return stepFraction_;
173 }
174 
175 
177 {
178  return origProc_;
179 }
180 
181 
183 {
184  return origProc_;
185 }
186 
187 
189 {
190  return origId_;
191 }
192 
193 
195 {
196  return origId_;
197 }
198 
199 
201 (
202  const polyMesh& mesh
203 ) const
204 {
205  return tetIndices(celli_, tetFacei_, tetPti_);
206 }
207 
208 
210 (
211  const polyMesh& mesh
212 ) const
213 {
214  if (mesh.moving() && stepFraction_ != 1)
215  {
216  return movingTetTransform(mesh, 0)[0];
217  }
218  else
219  {
220  return stationaryTetTransform(mesh);
221  }
222 }
223 
224 
226 {
227  return currentTetIndices(mesh).faceTri(mesh).normal();
228 }
229 
230 
231 inline bool Foam::particle::onFace() const
232 {
233  return facei_ >= 0;
234 }
235 
236 
237 inline bool Foam::particle::onInternalFace(const polyMesh& mesh) const
238 {
239  return onFace() && mesh.isInternalFace(facei_);
240 }
241 
242 
243 inline bool Foam::particle::onBoundaryFace(const polyMesh& mesh) const
244 {
245  return onFace() && !mesh.isInternalFace(facei_);
246 }
247 
248 
249 inline Foam::label Foam::particle::patch(const polyMesh& mesh) const
250 {
251  return onFace() ? mesh.boundaryMesh().whichPatch(facei_) : -1;
252 }
253 
254 
256 {
257  return currentTetTransform(mesh) & coordinates_;
258 }
259 
260 
261 inline void Foam::particle::reset(const scalar stepFraction)
262 {
263  stepFraction_ = stepFraction;
264  stepFractionBehind_ = 0;
265  nTracksBehind_ = 0;
266 }
267 
268 
270 (
271  const polyMesh& mesh,
272  vector& normal,
273  vector& displacement
274 ) const
275 {
276  if (!onBoundaryFace(mesh))
277  {
279  << "Patch data was requested for a particle that isn't on a patch"
280  << exit(FatalError);
281  }
282 
283  if (mesh.moving() && stepFraction_ != 1)
284  {
285  Pair<vector> centre, base, vertex1, vertex2;
286  movingTetGeometry(mesh, 1, centre, base, vertex1, vertex2);
287 
288  normal = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
289 
290  // Interpolate the motion of the three face vertices to the current
291  // coordinates
292  displacement =
293  coordinates_.b()*base[1]
294  + coordinates_.c()*vertex1[1]
295  + coordinates_.d()*vertex2[1];
296  }
297  else
298  {
299  vector centre, base, vertex1, vertex2;
300  stationaryTetGeometry(mesh, centre, base, vertex1, vertex2);
301 
302  normal = triPointRef(base, vertex1, vertex2).normal();
303 
304  displacement = Zero;
305  }
306 }
307 
308 
309 // ************************************************************************* //
Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can represent a barycentric transfo...
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:65
label face() const
Return current face particle is on otherwise -1.
Definition: particleI.H:152
vector normal(const polyMesh &mesh) const
Return the normal of the tri on tetFacei_ for the.
Definition: particleI.H:225
label tetPt() const
Return current tet face particle is in.
Definition: particleI.H:146
label origProc() const
Return the originating processor ID.
Definition: particleI.H:176
const barycentric & coordinates() const
Return current particle coordinates.
Definition: particleI.H:128
label cell() const
Return current cell particle is in.
Definition: particleI.H:134
void reset(const scalar stepFraction)
Set the step fraction and clear the behind data in preparation.
Definition: particleI.H:261
label patch(const polyMesh &mesh) const
Return the index of patch that the particle is on.
Definition: particleI.H:249
scalar stepFraction() const
Return the fraction of time-step completed.
Definition: particleI.H:164
bool onFace() const
Is the particle on a face?
Definition: particleI.H:231
void patchData(const polyMesh &mesh, vector &normal, vector &displacement) const
Get the normal and displacement of the current patch location.
Definition: particleI.H:270
bool onBoundaryFace(const polyMesh &mesh) const
Is the particle on a boundary face?
Definition: particleI.H:243
tetIndices currentTetIndices(const polyMesh &mesh) const
Return the indices of the current tet that the.
Definition: particleI.H:201
bool onInternalFace(const polyMesh &mesh) const
Is the particle on an internal face?
Definition: particleI.H:237
vector position(const polyMesh &mesh) const
Return current particle position.
Definition: particleI.H:255
label getNewParticleIndex() const
Get unique particle creation id.
Definition: particleI.H:114
label origId() const
Return the particle ID on the originating processor.
Definition: particleI.H:188
barycentricTensor currentTetTransform(const polyMesh &mesh) const
Return the current tet transformation tensor.
Definition: particleI.H:210
label tetFace() const
Return current tet face particle is in.
Definition: particleI.H:140
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:404
bool moving() const
Is mesh moving.
Definition: polyMesh.H:476
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:82
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
#define WarningInFunction
Report a warning using Foam::Warning.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Field< vector > vectorField
Specialisation of Field<T> for vector.
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
BarycentricTensor< scalar > barycentricTensor
A scalar version of the templated BarycentricTensor.
error FatalError
static const label labelMax
Definition: label.H:62
face triFace(3)