particleTemplates.C
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-2018 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 "IOPosition.H"
27 
28 #include "cyclicPolyPatch.H"
29 #include "cyclicAMIPolyPatch.H"
30 #include "cyclicACMIPolyPatch.H"
32 #include "processorPolyPatch.H"
33 #include "symmetryPlanePolyPatch.H"
34 #include "symmetryPolyPatch.H"
35 #include "wallPolyPatch.H"
36 #include "wedgePolyPatch.H"
37 #include "meshTools.H"
38 
39 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
40 
41 template<class TrackCloudType>
42 void Foam::particle::readFields(TrackCloudType& c)
43 {
44  bool valid = c.size();
45 
46  IOobject procIO(c.fieldIOobject("origProcId", IOobject::MUST_READ));
47 
48  bool haveFile = procIO.typeHeaderOk<IOField<label>>(true);
49 
50  IOField<label> origProcId(procIO, valid && haveFile);
51  c.checkFieldIOobject(c, origProcId);
53  (
54  c.fieldIOobject("origId", IOobject::MUST_READ),
55  valid && haveFile
56  );
57  c.checkFieldIOobject(c, origId);
58 
59  label i = 0;
60  forAllIter(typename TrackCloudType, c, iter)
61  {
62  particle& p = iter();
63 
64  p.origProc_ = origProcId[i];
65  p.origId_ = origId[i];
66  i++;
67  }
68 }
69 
70 
71 template<class TrackCloudType>
72 void Foam::particle::writeFields(const TrackCloudType& c)
73 {
74  label np = c.size();
75 
77  ioP.write(np > 0);
78 
80  (
81  c.fieldIOobject("origProcId", IOobject::NO_READ),
82  np
83  );
85  (
86  c.fieldIOobject("origId", IOobject::NO_READ),
87  np
88  );
89 
90  label i = 0;
91  forAllConstIter(typename TrackCloudType, c, iter)
92  {
93  origProc[i] = iter().origProc_;
94  origId[i] = iter().origId_;
95  i++;
96  }
97 
98  origProc.write(np > 0);
99  origId.write(np > 0);
100 }
101 
102 
103 template<class TrackCloudType>
105 (
106  const vector& direction,
107  TrackCloudType& cloud,
108  trackingData& td
109 )
110 {
111  typename TrackCloudType::particleType& p =
112  static_cast<typename TrackCloudType::particleType&>(*this);
113  typename TrackCloudType::particleType::trackingData& ttd =
114  static_cast<typename TrackCloudType::particleType::trackingData&>(td);
115 
116  if (!onFace())
117  {
118  return;
119  }
120  else if (onInternalFace())
121  {
122  changeCell();
123  }
124  else if (onBoundaryFace())
125  {
126  changeToMasterPatch();
127 
128  if (!p.hitPatch(cloud, ttd))
129  {
130  const polyPatch& patch = mesh_.boundaryMesh()[p.patch()];
131 
132  if (isA<wedgePolyPatch>(patch))
133  {
134  p.hitWedgePatch(cloud, ttd);
135  }
136  else if (isA<symmetryPlanePolyPatch>(patch))
137  {
138  p.hitSymmetryPlanePatch(cloud, ttd);
139  }
140  else if (isA<symmetryPolyPatch>(patch))
141  {
142  p.hitSymmetryPatch(cloud, ttd);
143  }
144  else if (isA<cyclicPolyPatch>(patch))
145  {
146  p.hitCyclicPatch(cloud, ttd);
147  }
148  else if (isA<cyclicACMIPolyPatch>(patch))
149  {
150  p.hitCyclicACMIPatch(cloud, ttd, direction);
151  }
152  else if (isA<cyclicAMIPolyPatch>(patch))
153  {
154  p.hitCyclicAMIPatch(cloud, ttd, direction);
155  }
156  else if (isA<cyclicRepeatAMIPolyPatch>(patch))
157  {
158  p.hitCyclicRepeatAMIPatch(cloud, ttd, direction);
159  }
160  else if (isA<processorPolyPatch>(patch))
161  {
162  p.hitProcessorPatch(cloud, ttd);
163  }
164  else if (isA<wallPolyPatch>(patch))
165  {
166  p.hitWallPatch(cloud, ttd);
167  }
168  else
169  {
170  td.keepParticle = false;
171  }
172  }
173  }
174 }
175 
176 
177 template<class TrackCloudType>
179 (
180  const vector& direction,
181  const scalar fraction,
182  TrackCloudType& cloud,
183  trackingData& td
184 )
185 {
186  trackToFace(direction, fraction);
187 
188  hitFace(direction, cloud, td);
189 }
190 
191 
192 template<class TrackCloudType>
193 bool Foam::particle::hitPatch(TrackCloudType&, trackingData&)
194 {
195  return false;
196 }
197 
198 
199 template<class TrackCloudType>
200 void Foam::particle::hitWedgePatch(TrackCloudType& cloud, trackingData& td)
201 {
203  << "Hitting a wedge patch should not be possible."
204  << abort(FatalError);
205 
206  hitSymmetryPatch(cloud, td);
207 }
208 
209 
210 template<class TrackCloudType>
212 (
213  TrackCloudType& cloud,
214  trackingData& td
215 )
216 {
217  hitSymmetryPatch(cloud, td);
218 }
219 
220 
221 template<class TrackCloudType>
223 {
224  const vector nf = normal();
225  transformProperties(I - 2.0*nf*nf);
226 }
227 
228 
229 template<class TrackCloudType>
231 {
232  const cyclicPolyPatch& cpp =
233  static_cast<const cyclicPolyPatch&>(mesh_.boundaryMesh()[patch()]);
234  const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
235  const label receiveFacei = receiveCpp.whichFace(facei_);
236 
237  // Set the topology
238  facei_ = tetFacei_ = cpp.transformGlobalFace(facei_);
239  celli_ = mesh_.faceOwner()[facei_];
240  // See note in correctAfterParallelTransfer for tetPti addressing ...
241  tetPti_ = mesh_.faces()[tetFacei_].size() - 1 - tetPti_;
242 
243  // Reflect to account for the change of triangle orientation in the new cell
244  reflect();
245 
246  // Transform the properties
247  if (!receiveCpp.parallel())
248  {
249  const tensor& T =
250  (
251  receiveCpp.forwardT().size() == 1
252  ? receiveCpp.forwardT()[0]
253  : receiveCpp.forwardT()[receiveFacei]
254  );
256  }
257  else if (receiveCpp.separated())
258  {
259  const vector& s =
260  (
261  (receiveCpp.separation().size() == 1)
262  ? receiveCpp.separation()[0]
263  : receiveCpp.separation()[receiveFacei]
264  );
266  }
267 }
268 
269 
270 template<class TrackCloudType>
272 (
273  TrackCloudType&,
274  trackingData& td,
275  const vector& direction
276 )
277 {
278  vector pos = position();
279 
280  const cyclicAMIPolyPatch& cpp =
281  static_cast<const cyclicAMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
282  const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
283  const label sendFacei = cpp.whichFace(facei_);
284  const labelPair receiveIs = cpp.pointAMIAndFace(sendFacei, direction, pos);
285  const label receiveAMIi = receiveIs.first();
286  const label receiveFacei = receiveIs.second();
287 
288  if (receiveFacei < 0)
289  {
290  // If the patch face of the particle is not known assume that the
291  // particle is lost and mark it to be deleted.
292  td.keepParticle = false;
294  << "Particle lost across " << cyclicAMIPolyPatch::typeName
295  << " patches " << cpp.name() << " and " << receiveCpp.name()
296  << " at position " << pos << endl;
297  return;
298  }
299 
300  // Set the topology
301  facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
302 
303  // Locate the particle on the receiving side
304  vector directionT = direction;
305  cpp.reverseTransformDirection(directionT, sendFacei);
306  locate
307  (
308  pos,
309  &directionT,
310  mesh_.faceOwner()[facei_],
311  false,
312  "Particle crossed between " + cyclicAMIPolyPatch::typeName +
313  " patches " + cpp.name() + " and " + receiveCpp.name() +
314  " to a location outside of the mesh."
315  );
316 
317  // The particle must remain associated with a face for the tracking to
318  // register as incomplete
319  facei_ = tetFacei_;
320 
321  // Transform the properties
322  if (!receiveCpp.parallel())
323  {
324  const tensor& T =
325  (
326  receiveCpp.forwardT().size() == 1
327  ? receiveCpp.forwardT()[0]
328  : receiveCpp.forwardT()[receiveFacei]
329  );
331  }
332  else if (receiveCpp.separated())
333  {
334  const vector& s =
335  (
336  (receiveCpp.separation().size() == 1)
337  ? receiveCpp.separation()[0]
338  : receiveCpp.separation()[receiveFacei]
339  );
341  }
342  const vectorTensorTransform& T =
343  cpp.owner()
344  ? cpp.AMITransforms()[receiveAMIi]
345  : cpp.neighbPatch().AMITransforms()[receiveAMIi];
346  if (T.hasR())
347  {
348  transformProperties(T.R());
349  }
350  else if (T.t() != vector::zero)
351  {
352  transformProperties(T.t());
353  }
354 }
355 
356 
357 template<class TrackCloudType>
359 (
360  TrackCloudType& cloud,
361  trackingData& td,
362  const vector& direction
363 )
364 {
365  const cyclicACMIPolyPatch& cpp =
366  static_cast<const cyclicACMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
367 
368  const label localFacei = cpp.whichFace(facei_);
369 
370  // If the mask is within the patch tolerance at either end, then we can
371  // assume an interaction with the appropriate part of the ACMI pair.
372  const scalar mask = cpp.mask()[localFacei];
373  bool couple = mask >= 1 - cpp.tolerance();
374  bool nonOverlap = mask <= cpp.tolerance();
375 
376  // If the mask is an intermediate value, then we search for a location on
377  // the other side of the AMI. If we can't find a location, then we assume
378  // that we have hit the non-overlap patch.
379  if (!couple && !nonOverlap)
380  {
381  vector pos = position();
382  couple = cpp.pointAMIAndFace(localFacei, direction, pos).first() >= 0;
383  nonOverlap = !couple;
384  }
385 
386  if (couple)
387  {
388  hitCyclicAMIPatch(cloud, td, direction);
389  }
390  else
391  {
392  // Move to the face associated with the non-overlap patch and redo the
393  // face interaction.
394  tetFacei_ = facei_ = cpp.nonOverlapPatch().start() + localFacei;
395  hitFace(direction, cloud, td);
396  }
397 }
398 
399 
400 template<class TrackCloudType>
402 (
403  TrackCloudType& cloud,
404  trackingData& td,
405  const vector& direction
406 )
407 {
408  hitCyclicAMIPatch(cloud, td, direction);
409 }
410 
411 
412 template<class TrackCloudType>
414 {}
415 
416 
417 template<class TrackCloudType>
419 {}
420 
421 
422 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
void hitCyclicACMIPatch(TrackCloudType &, trackingData &, const vector &)
Overridable function to handle the particle hitting a.
label origProc() const
Return the originating processor ID.
Definition: particleI.H:178
Cyclic patch for Arbitrarily Coupled Mesh Interface (ACMI)
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
const word & name() const
Return name.
bool typeHeaderOk(const bool checkType=true)
Read header (uses typeFilePath to find file) and check header.
virtual bool separated() const
Are the planes separated.
bool onBoundaryFace() const
Is the particle on a boundary face?
Definition: particleI.H:271
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Vector-tensor class used to perform translations and rotations in 3D space.
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:453
uint8_t direction
Definition: direction.H:45
scalar trackToFace(const vector &displacement, const scalar fraction)
As particle::track, but stops when a face is hit.
Definition: particle.C:625
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
labelPair pointAMIAndFace(const label facei, const vector &n, point &p) const
Return the transform and face indices on neighbour patch which.
void hitWallPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a wallPatch.
void hitSymmetryPlanePatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a.
virtual void reverseTransformDirection(vector &d, const label facei) const
Transform a patch-based direction from this side to nbr side.
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
Definition: particle.C:972
void hitSymmetryPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a symmetryPatch.
virtual const vectorField & separation() const
If the planes are separated the separation vector.
virtual const tensorField & forwardT() const
Return face transformation tensor.
Base particle class.
Definition: particle.H:84
const polyPatch & nonOverlapPatch() const
Return a const reference to the non-overlapping patch.
const List< vectorTensorTransform > & AMITransforms() const
Return a reference to the AMI transforms.
vector normal() const
Return the normal of the tri on tetFacei_ for the.
Definition: particleI.H:253
virtual bool parallel() const
Are the cyclic planes parallel.
dimensionedScalar pos(const dimensionedScalar &ds)
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))
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
void hitFace(const vector &direction, TrackCloudType &cloud, trackingData &td)
Hit the current face. If the current face is internal than this.
static const Identity< scalar > I
Definition: Identity.H:93
void hitProcessorPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a processorPatch.
Cyclic plane patch.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1041
void hitCyclicPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a cyclicPatch.
Cyclic patch for Arbitrary Mesh Interface (AMI)
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1028
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void hitCyclicAMIPatch(TrackCloudType &, trackingData &, const vector &)
Overridable function to handle the particle hitting a cyclicAMIPatch.
label origId() const
Return the particle ID on the originating processor.
Definition: particleI.H:190
virtual bool write(const bool valid=true) const
Write using setting from DB.
Definition: IOPosition.C:51
bool onFace() const
Is the particle on a face?
Definition: particleI.H:259
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
const Type & second() const
Return second.
Definition: Pair.H:99
void hitWedgePatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a wedgePatch.
const scalarField & mask() const
Mask field where 1 = overlap, 0 = no-overlap.
#define WarningInFunction
Report a warning using Foam::Warning.
bool keepParticle
Flag to indicate whether to keep particle (false = delete)
Definition: particle.H:119
virtual const cyclicAMIPolyPatch & neighbPatch() const
Return a reference to the neighbour patch.
static void readFields(TrackCloudType &c)
Read the fields associated with the owner cloud.
virtual bool owner() const
Does this side own the patch?
const dimensionedScalar c
Speed of light in a vacuum.
static void writeFields(const TrackCloudType &c)
Write the fields associated with the owner cloud.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
bool hitPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a patch.
bool onInternalFace() const
Is the particle on an internal face?
Definition: particleI.H:265
label transformGlobalFace(const label facei) const
void trackToAndHitFace(const vector &direction, const scalar fraction, TrackCloudType &cloud, trackingData &td)
Convenience function. Cobines trackToFace and hitFace.
volScalarField & p
Helper IO class to read and write particle positions.
Definition: Cloud.H:55
const cyclicPolyPatch & neighbPatch() const
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
const Type & first() const
Return first.
Definition: Pair.H:87
A primitive field of type <T> with automated input and output.
Definition: IOField.H:50
vector position() const
Return current particle position.
Definition: particleI.H:283
void hitCyclicRepeatAMIPatch(TrackCloudType &, trackingData &, const vector &)
Overridable function to handle the particle hitting an.
label patch() const
Return the index of patch that the particle is on.
Definition: particleI.H:277
static scalar tolerance()
Overlap tolerance.
label whichFace(const label l) const
Return label of face in patch from global face label.
Definition: polyPatch.H:377