wallBoundedParticleTemplates.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 Copyright (C) 2017-2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
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 "wallBoundedParticle.H"
30
31// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32
33template<class TrackCloudType>
35(
36 TrackCloudType& cloud,
37 trackingData& td,
38 const scalar trackFraction
39)
40{
41 typename TrackCloudType::particleType& p =
42 static_cast<typename TrackCloudType::particleType&>(*this);
43 typename TrackCloudType::particleType::trackingData& ttd =
44 static_cast<typename TrackCloudType::particleType::trackingData&>(td);
45
46 if (!mesh().isInternalFace(face()))
47 {
48 label origFacei = face();
49 label patchi = patch();
50
51 // Did patch interaction model switch patches?
52 // Note: recalculate meshEdgeStart_, diagEdge_!
53 if (face() != origFacei)
54 {
55 patchi = patch();
56 }
57
58 const polyPatch& patch = mesh().boundaryMesh()[patchi];
59
60 if (isA<processorPolyPatch>(patch))
61 {
62 p.hitProcessorPatch(cloud, ttd);
63 }
64 else if (isA<wallPolyPatch>(patch))
65 {
66 p.hitWallPatch(cloud, ttd);
67 }
68 else
69 {
70 td.keepParticle = false;
71 }
72 }
73}
74
75
76// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
77
78template<class TrackCloudType>
80(
81 TrackCloudType& cloud,
82 trackingData& td,
83 const vector& endPosition
84)
85{
86 // Track particle to a given position and returns 1.0 if the
87 // trajectory is completed without hitting a face otherwise
88 // stops at the face and returns the fraction of the trajectory
89 // completed.
90 // on entry 'stepFraction()' should be set to the fraction of the
91 // time-step at which the tracking starts.
92
93 // Are we on a track face? If not we do a topological walk.
94
95 // Particle:
96 // - cell_ always set
97 // - tetFace_, tetPt_ always set (these identify tet particle is in)
98 // - optionally meshEdgeStart_ or diagEdge_ set (edge particle is on)
99
100 //checkInside();
101 //checkOnTriangle(localPosition_);
102 //if (meshEdgeStart_ != -1 || diagEdge_ != -1)
103 //{
104 // checkOnEdge();
105 //}
106
107 scalar trackFraction = 0.0;
108
109 if (!td.isWallPatch_[tetFace()])
110 {
111 // Don't track across face. Just walk in cell. Particle is on
112 // mesh edge (as indicated by meshEdgeStart_).
113 const edge meshEdge(currentEdge());
114
115 // If internal face check whether to go to neighbour cell or just
116 // check to the other internal tet on the edge.
117 if (mesh().isInternalFace(tetFace()))
118 {
119 label nbrCelli =
120 (
121 this->cell() == mesh().faceOwner()[face()]
122 ? mesh().faceNeighbour()[face()]
123 : mesh().faceOwner()[face()]
124 );
125 // Check angle to nbrCell tet. Is it in the direction of the
126 // endposition? i.e. since volume of nbr tet is positive the
127 // tracking direction should be into the tet.
128 tetIndices nbrTi(nbrCelli, tetFace(), tetPt());
129
130 const bool posVol = (nbrTi.tet(mesh()).mag() > 0);
131 const vector path(endPosition - localPosition_);
132
133 if (posVol == ((nbrTi.faceTri(mesh()).areaNormal() & path) < 0))
134 {
135 // Change into nbrCell. No need to change tetFace, tetPt.
136 //Pout<< " crossed from cell:" << celli_
137 // << " into " << nbrCelli << endl;
138 this->cell() = nbrCelli;
139 patchInteraction(cloud, td, trackFraction);
140 }
141 else
142 {
143 // Walk to other face on edge. Changes tetFace, tetPt but not
144 // cell.
145 crossEdgeConnectedFace(meshEdge);
146 patchInteraction(cloud, td, trackFraction);
147 }
148 }
149 else
150 {
151 // Walk to other face on edge. This might give loop since
152 // particle should have been removed?
153 crossEdgeConnectedFace(meshEdge);
154 patchInteraction(cloud, td, trackFraction);
155 }
156 }
157 else
158 {
159 // We're inside a tet on the wall. Check if the current tet is
160 // the one to cross. If not we cross into the neighbouring triangle.
161
162 if (mesh().isInternalFace(tetFace()))
163 {
165 << "Can only track on boundary faces."
166 << " Face:" << tetFace()
167 << " at:" << mesh().faceCentres()[tetFace()]
168 << abort(FatalError);
169 }
170
171 const triFace tri(currentTetIndices().faceTriIs(mesh(), false));
172 const vector n = tri.unitNormal(mesh().points());
173
174 point projectedEndPosition = endPosition;
175
176 const bool posVol = (currentTetIndices().tet(mesh()).mag() > 0);
177
178 if (!posVol)
179 {
180 // Negative tet volume. Track back by setting the end point
181 projectedEndPosition =
182 localPosition_ - (endPosition - localPosition_);
183
184 // Make sure to use a large enough vector to cross the negative
185 // face. Bit overkill.
186 const vector d(endPosition - localPosition_);
187 const scalar magD(mag(d));
188 if (magD > ROOTVSMALL)
189 {
190 // Get overall mesh bounding box
191 treeBoundBox meshBb(mesh().bounds());
192 // Extend to make 3D
193 meshBb.inflate(ROOTSMALL);
194
195 // Create vector guaranteed to cross mesh bounds
196 projectedEndPosition = localPosition_ - meshBb.mag()*d/magD;
197
198 // Clip to mesh bounds
199 point intPt;
200 direction intPtBits;
201 bool ok = meshBb.intersects
202 (
203 projectedEndPosition,
204 localPosition_ - projectedEndPosition,
205 projectedEndPosition,
206 localPosition_,
207 intPt,
208 intPtBits
209 );
210 if (ok)
211 {
212 // Should always be the case
213 projectedEndPosition = intPt;
214 }
215 }
216 }
217
218 // Remove normal component
219 {
220 const point& basePt = mesh().points()[tri[0]];
221 projectedEndPosition -= ((projectedEndPosition - basePt)&n)*n;
222 }
223
224
225 bool doTrack = false;
226 if (meshEdgeStart_ == -1 && diagEdge_ == -1)
227 {
228 // We're starting and not yet on an edge.
229 doTrack = true;
230 }
231 else
232 {
233 // See if the current triangle has got a point on the
234 // correct side of the edge.
235 doTrack = isTriAlongTrack(n, projectedEndPosition);
236 }
237
238
239 if (doTrack)
240 {
241 // Track across triangle. Return triangle edge crossed.
242 label triEdgei = -1;
243 trackFraction = trackFaceTri(n, projectedEndPosition, triEdgei);
244
245 if (triEdgei == -1)
246 {
247 // Reached endpoint
248 //checkInside();
249 diagEdge_ = -1;
250 meshEdgeStart_ = -1;
251 return trackFraction;
252 }
253
254 const tetIndices ti(currentTetIndices());
255 const triFace trif(ti.triIs(mesh(), false));
256 // Triangle (faceTriIs) gets constructed from
257 // f[faceBasePtI_],
258 // f[facePtAI_],
259 // f[facePtBI_]
260 //
261 // So edge indices are:
262 // 0 : edge between faceBasePtI_ and facePtAI_
263 // 1 : edge between facePtAI_ and facePtBI_ (is always a real edge)
264 // 2 : edge between facePtBI_ and faceBasePtI_
265
266 const Foam::face& f = mesh().faces()[ti.face()];
267 const label fp0 = trif[0];
268
269 if (triEdgei == 0)
270 {
271 if (trif[1] == f.fcIndex(fp0))
272 {
273 //Pout<< "Real edge." << endl;
274 diagEdge_ = -1;
275 meshEdgeStart_ = fp0;
276 //checkOnEdge();
277 crossEdgeConnectedFace(currentEdge());
278 patchInteraction(cloud, td, trackFraction);
279 }
280 else if (trif[1] == f.rcIndex(fp0))
281 {
282 //Note: should not happen since boundary face so owner
283 //Pout<< "Real edge." << endl;
285 << abort(FatalError);
286
287 diagEdge_ = -1;
288 meshEdgeStart_ = f.rcIndex(fp0);
289 //checkOnEdge();
290 crossEdgeConnectedFace(currentEdge());
291 patchInteraction(cloud, td, trackFraction);
292 }
293 else
294 {
295 // Get index of triangle on other side of edge.
296 diagEdge_ = trif[1] - fp0;
297 if (diagEdge_ < 0)
298 {
299 diagEdge_ += f.size();
300 }
301 meshEdgeStart_ = -1;
302 //checkOnEdge();
303 crossDiagonalEdge();
304 }
305 }
306 else if (triEdgei == 1)
307 {
308 //Pout<< "Real edge." << endl;
309 diagEdge_ = -1;
310 meshEdgeStart_ = trif[1];
311 //checkOnEdge();
312 crossEdgeConnectedFace(currentEdge());
313 patchInteraction(cloud, td, trackFraction);
314 }
315 else // if (triEdgei == 2)
316 {
317 if (trif[2] == f.rcIndex(fp0))
318 {
319 //Pout<< "Real edge." << endl;
320 diagEdge_ = -1;
321 meshEdgeStart_ = trif[2];
322 //checkOnEdge();
323 crossEdgeConnectedFace(currentEdge());
324 patchInteraction(cloud, td, trackFraction);
325 }
326 else if (trif[2] == f.fcIndex(fp0))
327 {
328 //Note: should not happen since boundary face so owner
329 //Pout<< "Real edge." << endl;
331
332 diagEdge_ = -1;
333 meshEdgeStart_ = fp0;
334 //checkOnEdge();
335 crossEdgeConnectedFace(currentEdge());
336 patchInteraction(cloud, td, trackFraction);
337 }
338 else
339 {
340 //Pout<< "Triangle edge." << endl;
341 // Get index of triangle on other side of edge.
342 diagEdge_ = trif[2] - fp0;
343 if (diagEdge_ < 0)
344 {
345 diagEdge_ += f.size();
346 }
347 meshEdgeStart_ = -1;
348 //checkOnEdge();
349 crossDiagonalEdge();
350 }
351 }
352 }
353 else
354 {
355 // Current tet is not the right one. Check the neighbour tet.
356
357 if (meshEdgeStart_ != -1)
358 {
359 // Particle is on mesh edge so change into other face on cell
360 crossEdgeConnectedFace(currentEdge());
361 //checkOnEdge();
362 patchInteraction(cloud, td, trackFraction);
363 }
364 else
365 {
366 // Particle is on diagonal edge so change into the other
367 // triangle.
368 crossDiagonalEdge();
369 //checkOnEdge();
370 }
371 }
372 }
373
374 //checkInside();
375
376 return trackFraction;
377}
378
379
380template<class TrackCloudType>
382(
383 TrackCloudType& cloud,
384 trackingData& td
385)
386{
387 // Switch particle
388 td.switchProcessor = true;
389
390 // Adapt edgeStart_ for other side.
391 // E.g. if edgeStart_ is 1 then the edge is between vertex 1 and 2 so
392 // on the other side between 2 and 3 so edgeStart_ should be
393 // f.size()-edgeStart_-1.
394
395 const Foam::face& f = mesh().faces()[face()];
396
397 if (meshEdgeStart_ != -1)
398 {
399 meshEdgeStart_ = f.size() - meshEdgeStart_-1;
400 }
401 else
402 {
403 // diagEdge_ is relative to faceBasePt
404 diagEdge_ = f.size() - diagEdge_;
405 }
406}
407
408
409template<class TrackCloudType>
411(
412 TrackCloudType& cloud,
413 trackingData& td
414)
415{}
416
417
418template<class TrackCloudType>
420{
421 if (!c.size())
422 {
423 return;
424 }
425
427
428 IOField<point> localPosition
429 (
430 c.fieldIOobject("position", IOobject::MUST_READ)
431 );
432 c.checkFieldIOobject(c, localPosition);
433
434 IOField<label> meshEdgeStart
435 (
436 c.fieldIOobject("meshEdgeStart", IOobject::MUST_READ)
437 );
438 c.checkFieldIOobject(c, meshEdgeStart);
439
440 IOField<label> diagEdge
441 (
442 c.fieldIOobject("diagEdge", IOobject::MUST_READ)
443 );
444 c.checkFieldIOobject(c, diagEdge);
445
446 label i = 0;
447 for (wallBoundedParticle& p : c)
448 {
449 p.localPosition_ = localPosition[i];
450 p.meshEdgeStart_ = meshEdgeStart[i];
451 p.diagEdge_ = diagEdge[i];
452
453 ++i;
454 }
455}
456
457
458template<class TrackCloudType>
459void Foam::wallBoundedParticle::writeFields(const TrackCloudType& c)
460{
462
463 label np = c.size();
464
465 IOField<point> localPosition
466 (
467 c.fieldIOobject("position", IOobject::NO_READ),
468 np
469 );
470 IOField<label> meshEdgeStart
471 (
472 c.fieldIOobject("meshEdgeStart", IOobject::NO_READ),
473 np
474 );
475 IOField<label> diagEdge
476 (
477 c.fieldIOobject("diagEdge", IOobject::NO_READ),
478 np
479 );
480
481 label i = 0;
482 for (const wallBoundedParticle& p : c)
483 {
484 localPosition[i] = p.localPosition_;
485 meshEdgeStart[i] = p.meshEdgeStart_;
486 diagEdge[i] = p.diagEdge_;
487
488 ++i;
489 }
490
491 localPosition.write();
492 meshEdgeStart.write();
493 diagEdge.write();
494}
495
496
497// ************************************************************************* //
label n
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
const PatchInteractionModel< KinematicCloud< CloudType > > & patchInteraction() const
Return const-access to the patch interaction model.
label rcIndex(const label i) const noexcept
Definition: UListI.H:67
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
A cloud is a registry collection of lagrangian particles.
Definition: cloud.H:60
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Reads fields from the time directories and adds them to the mesh database for further post-processing...
Definition: readFields.H:158
bool switchProcessor
Flag to switch processor.
Definition: particle.H:102
bool keepParticle
Flag to indicate whether to keep particle (false = delete)
Definition: particle.H:105
label face() const noexcept
Return current face particle is on otherwise -1.
Definition: particleI.H:185
const polyMesh & mesh() const noexcept
Return the mesh database.
Definition: particleI.H:137
label patch() const
Return the index of patch that the particle is on.
Definition: particleI.H:308
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const vectorField & faceCentres() const
virtual bool write(const bool valid=true) const
Write using setting from DB.
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:84
label face() const noexcept
Return the face index.
Definition: tetIndices.H:139
triFace triIs(const polyMesh &mesh, const bool warn=true) const
Definition: tetIndicesI.H:95
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet.
Definition: tetIndicesI.H:134
triPointRef faceTri(const polyMesh &mesh) const
Definition: tetIndicesI.H:149
scalar mag() const
Return volume.
Definition: tetrahedronI.H:172
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:72
vector unitNormal(const UList< point > &points) const
The unit normal.
Definition: triFaceI.H:209
Class used to pass tracking data to the trackToFace function.
Particle class that tracks on triangles of boundary faces. Use trackToEdge similar to trackToFace on ...
void hitProcessorPatch(TrackCloudType &cloud, trackingData &td)
scalar trackToEdge(TrackCloudType &cloud, trackingData &td, const vector &endPosition)
Equivalent of trackToFace.
void hitWallPatch(TrackCloudType &cloud, trackingData &td)
Overridable function to handle the particle hitting a wallPatch.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
volScalarField & p
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
uint8_t direction
Definition: direction.H:56
error FatalError
labelList f(nPoints)
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(rndGen, 1e-3))