ParticleCollector.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) 2012-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 "ParticleCollector.H"
27 #include "Pstream.H"
28 #include "surfaceWriter.H"
29 #include "unitConversion.H"
30 #include "Random.H"
31 #include "triangle.H"
32 #include "cloud.H"
33 
34 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
35 
36 template<class CloudType>
38 (
39  const faceList& faces,
40  const Field<point>& points,
41  const Field<scalar>& area
42 )
43 {
44  // Create the output file if not already created
45  if (log_)
46  {
47  if (debug)
48  {
49  Info<< "Creating output file" << endl;
50  }
51 
52  if (Pstream::master())
53  {
54  // Create directory if does not exist
55  mkDir(this->writeTimeDir());
56 
57  // Open new file at start up
58  outputFilePtr_.reset
59  (
60  new OFstream(this->writeTimeDir()/(type() + ".dat"))
61  );
62 
63  outputFilePtr_()
64  << "# Source : " << type() << nl
65  << "# Bins : " << faces.size() << nl
66  << "# Total area : " << sum(area) << nl;
67 
68  outputFilePtr_()
69  << "# Geometry :" << nl
70  << '#'
71  << tab << "Bin"
72  << tab << "(Centre_x Centre_y Centre_z)"
73  << tab << "Area"
74  << nl;
75 
76  forAll(faces, i)
77  {
78  outputFilePtr_()
79  << '#'
80  << tab << i
81  << tab << faces[i].centre(points)
82  << tab << area[i]
83  << nl;
84  }
85 
86  outputFilePtr_()
87  << '#' << nl
88  << "# Output format:" << nl;
89 
90  forAll(faces, i)
91  {
92  word id = Foam::name(i);
93  word binId = "bin_" + id;
94 
95  outputFilePtr_()
96  << '#'
97  << tab << "Time"
98  << tab << binId
99  << tab << "mass[" << id << "]"
100  << tab << "massFlowRate[" << id << "]"
101  << endl;
102  }
103  }
104  }
105 }
106 
107 
108 template<class CloudType>
110 (
111  const List<Field<point>>& polygons
112 )
113 {
114  mode_ = mtPolygon;
115 
116  label nPoints = 0;
117  forAll(polygons, polyI)
118  {
119  label np = polygons[polyI].size();
120  if (np < 3)
121  {
122  FatalIOErrorInFunction(this->coeffDict())
123  << "polygons must consist of at least 3 points"
124  << exit(FatalIOError);
125  }
126 
127  nPoints += np;
128  }
129 
130  label pointOffset = 0;
131  points_.setSize(nPoints);
132  faces_.setSize(polygons.size());
133  faceTris_.setSize(polygons.size());
134  area_.setSize(polygons.size());
135  forAll(faces_, facei)
136  {
137  const Field<point>& polyPoints = polygons[facei];
138  face f(identity(polyPoints.size()) + pointOffset);
139  UIndirectList<point>(points_, f) = polyPoints;
140  area_[facei] = f.mag(points_);
141 
142  DynamicList<face> tris;
143  f.triangles(points_, tris);
144  faceTris_[facei].transfer(tris);
145 
146  faces_[facei].transfer(f);
147 
148  pointOffset += polyPoints.size();
149  }
150 }
151 
152 
153 template<class CloudType>
155 {
156  mode_ = mtConcentricCircle;
157 
158  vector origin(this->coeffDict().lookup("origin"));
159 
160  this->coeffDict().lookup("radius") >> radius_;
161  nSector_ = readLabel(this->coeffDict().lookup("nSector"));
162 
163  label nS = nSector_;
164 
165  vector refDir;
166  if (nSector_ > 1)
167  {
168  refDir = this->coeffDict().lookup("refDir");
169  refDir -= normal_[0]*(normal_[0] & refDir);
170  refDir /= mag(refDir);
171  }
172  else
173  {
174  // Set 4 quadrants for single sector cases
175  nS = 4;
176 
177  vector tangent = Zero;
178  scalar magTangent = 0.0;
179 
180  Random rnd(1234);
181  while (magTangent < small)
182  {
183  vector v = rnd.sample01<vector>();
184 
185  tangent = v - (v & normal_[0])*normal_[0];
186  magTangent = mag(tangent);
187  }
188 
189  refDir = tangent/magTangent;
190  }
191 
192  scalar dTheta = 5.0;
193  scalar dThetaSector = 360.0/scalar(nS);
194  label intervalPerSector = max(1, ceil(dThetaSector/dTheta));
195  dTheta = dThetaSector/scalar(intervalPerSector);
196 
197  label nPointPerSector = intervalPerSector + 1;
198 
199  label nPointPerRadius = nS*(nPointPerSector - 1);
200  label nPoint = radius_.size()*nPointPerRadius;
201  label nFace = radius_.size()*nS;
202 
203  // Add origin
204  nPoint++;
205 
206  points_.setSize(nPoint);
207  faces_.setSize(nFace);
208  area_.setSize(nFace);
209 
210  coordSys_ = cylindricalCS("coordSys", origin, normal_[0], refDir, false);
211 
212  List<label> ptIDs(identity(nPointPerRadius));
213 
214  points_[0] = origin;
215 
216  // Points
217  forAll(radius_, radI)
218  {
219  label pointOffset = radI*nPointPerRadius + 1;
220 
221  for (label i = 0; i < nPointPerRadius; i++)
222  {
223  label pI = i + pointOffset;
224  point pCyl(radius_[radI], degToRad(i*dTheta), 0.0);
225  points_[pI] = coordSys_.globalPosition(pCyl);
226  }
227  }
228 
229  // Faces
230  DynamicList<label> facePts(2*nPointPerSector);
231  forAll(radius_, radI)
232  {
233  if (radI == 0)
234  {
235  for (label secI = 0; secI < nS; secI++)
236  {
237  facePts.clear();
238 
239  // Append origin point
240  facePts.append(0);
241 
242  for (label ptI = 0; ptI < nPointPerSector; ptI++)
243  {
244  label i = ptI + secI*(nPointPerSector - 1);
245  label id = ptIDs.fcIndex(i - 1) + 1;
246  facePts.append(id);
247  }
248 
249  label facei = secI + radI*nS;
250 
251  faces_[facei] = face(facePts);
252  area_[facei] = faces_[facei].mag(points_);
253  }
254  }
255  else
256  {
257  for (label secI = 0; secI < nS; secI++)
258  {
259  facePts.clear();
260 
261  label offset = (radI - 1)*nPointPerRadius + 1;
262 
263  for (label ptI = 0; ptI < nPointPerSector; ptI++)
264  {
265  label i = ptI + secI*(nPointPerSector - 1);
266  label id = offset + ptIDs.fcIndex(i - 1);
267  facePts.append(id);
268  }
269  for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
270  {
271  label i = ptI + secI*(nPointPerSector - 1);
272  label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
273  facePts.append(id);
274  }
275 
276  label facei = secI + radI*nS;
277 
278  faces_[facei] = face(facePts);
279  area_[facei] = faces_[facei].mag(points_);
280  }
281  }
282  }
283 }
284 
285 
286 template<class CloudType>
288 (
289  const point& p1,
290  const point& p2
291 ) const
292 {
293  forAll(faces_, facei)
294  {
295  const label facePoint0 = faces_[facei][0];
296 
297  const point& pf = points_[facePoint0];
298 
299  const scalar d1 = normal_[facei] & (p1 - pf);
300  const scalar d2 = normal_[facei] & (p2 - pf);
301 
302  if (sign(d1) == sign(d2))
303  {
304  // Did not cross polygon plane
305  continue;
306  }
307 
308  // Intersection point
309  const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
310 
311  // Identify if point is within the bounds of the face. Create triangles
312  // between the intersection point and each edge of the face. If all the
313  // triangle normals point in the same direction as the face normal, then
314  // the particle is within the face. Note that testing for pointHits on
315  // the face's decomposed triangles does not work due to ambiguity along
316  // the diagonals.
317  const face& f = faces_[facei];
318  const vector a = f.area(points_);
319  bool inside = true;
320  for (label i = 0; i < f.size(); ++ i)
321  {
322  const label j = f.fcIndex(i);
323  const triPointRef t(pIntersect, points_[f[i]], points_[f[j]]);
324  if ((a & t.area()) < 0)
325  {
326  inside = false;
327  break;
328  }
329  }
330 
331  // Add to the list of hits
332  if (inside)
333  {
334  hitFaceIDs_.append(facei);
335  }
336  }
337 }
338 
339 
340 template<class CloudType>
342 (
343  const point& p1,
344  const point& p2
345 ) const
346 {
347  label secI = -1;
348 
349  const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
350  const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
351 
352  if (sign(d1) == sign(d2))
353  {
354  // Did not cross plane
355  return;
356  }
357 
358  // Intersection point in cylindrical co-ordinate system
359  const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
360 
361  scalar r = pCyl[0];
362 
363  if (r < radius_.last())
364  {
365  label radI = 0;
366  while (r > radius_[radI])
367  {
368  radI++;
369  }
370 
371  if (nSector_ == 1)
372  {
373  secI = 4*radI;
374  }
375  else
376  {
377  scalar theta = pCyl[1] + constant::mathematical::pi;
378 
379  secI =
380  nSector_*radI
381  + floor
382  (
383  scalar(nSector_)*theta/constant::mathematical::twoPi
384  );
385  }
386  }
387 
388  if (secI != -1)
389  {
390  hitFaceIDs_.append(secI);
391  }
392 }
393 
394 
395 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
396 
397 template<class CloudType>
399 {
400  const fvMesh& mesh = this->owner().mesh();
401  const Time& time = mesh.time();
402  scalar timeNew = time.value();
403  scalar timeElapsed = timeNew - timeOld_;
404 
405  totalTime_ += timeElapsed;
406 
407  const scalar alpha = (totalTime_ - timeElapsed)/totalTime_;
408  const scalar beta = timeElapsed/totalTime_;
409 
410  forAll(faces_, facei)
411  {
412  massFlowRate_[facei] =
413  alpha*massFlowRate_[facei] + beta*mass_[facei]/timeElapsed;
414  massTotal_[facei] += mass_[facei];
415  }
416 
417  const label proci = Pstream::myProcNo();
418 
419  Info<< type() << " output:" << nl;
420 
421  Field<scalar> faceMassTotal(mass_.size(), 0.0);
422  this->getModelProperty("massTotal", faceMassTotal);
423 
424  Field<scalar> faceMassFlowRate(massFlowRate_.size(), 0.0);
425  this->getModelProperty("massFlowRate", faceMassFlowRate);
426 
427 
428  scalar sumTotalMass = 0.0;
429  scalar sumAverageMFR = 0.0;
430  forAll(faces_, facei)
431  {
432  scalarList allProcMass(Pstream::nProcs());
433  allProcMass[proci] = massTotal_[facei];
434  Pstream::gatherList(allProcMass);
435  faceMassTotal[facei] += sum(allProcMass);
436 
437  scalarList allProcMassFlowRate(Pstream::nProcs());
438  allProcMassFlowRate[proci] = massFlowRate_[facei];
439  Pstream::gatherList(allProcMassFlowRate);
440  faceMassFlowRate[facei] += sum(allProcMassFlowRate);
441 
442  sumTotalMass += faceMassTotal[facei];
443  sumAverageMFR += faceMassFlowRate[facei];
444 
445  if (outputFilePtr_.valid())
446  {
447  outputFilePtr_()
448  << time.timeName()
449  << tab << facei
450  << tab << faceMassTotal[facei]
451  << tab << faceMassFlowRate[facei]
452  << endl;
453  }
454  }
455 
456  Info<< " sum(total mass) = " << sumTotalMass << nl
457  << " sum(average mass flow rate) = " << sumAverageMFR << nl
458  << endl;
459 
460 
461  if (surfaceFormat_ != "none")
462  {
463  if (Pstream::master())
464  {
466 
467  writer->write
468  (
469  this->writeTimeDir(),
470  "collector",
471  points_,
472  faces_,
473  "massTotal",
474  faceMassTotal,
475  false
476  );
477 
478  writer->write
479  (
480  this->writeTimeDir(),
481  "collector",
482  points_,
483  faces_,
484  "massFlowRate",
485  faceMassFlowRate,
486  false
487  );
488  }
489  }
490 
491 
492  if (resetOnWrite_)
493  {
494  Field<scalar> dummy(faceMassTotal.size(), 0.0);
495  this->setModelProperty("massTotal", dummy);
496  this->setModelProperty("massFlowRate", dummy);
497 
498  timeOld_ = timeNew;
499  totalTime_ = 0.0;
500  }
501  else
502  {
503  this->setModelProperty("massTotal", faceMassTotal);
504  this->setModelProperty("massFlowRate", faceMassFlowRate);
505  }
506 
507  forAll(faces_, facei)
508  {
509  mass_[facei] = 0.0;
510  massTotal_[facei] = 0.0;
511  massFlowRate_[facei] = 0.0;
512  }
513 }
514 
515 
516 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
517 
518 template<class CloudType>
520 (
521  const dictionary& dict,
522  CloudType& owner,
523  const word& modelName
524 )
525 :
526  CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
527  mode_(mtUnknown),
528  parcelType_(this->coeffDict().lookupOrDefault("parcelType", -1)),
529  removeCollected_(this->coeffDict().lookup("removeCollected")),
530  points_(),
531  faces_(),
532  faceTris_(),
533  nSector_(0),
534  radius_(),
535  coordSys_(false),
536  normal_(),
537  negateParcelsOppositeNormal_
538  (
539  readBool(this->coeffDict().lookup("negateParcelsOppositeNormal"))
540  ),
541  surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
542  resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
543  totalTime_(0.0),
544  mass_(),
545  massTotal_(),
546  massFlowRate_(),
547  log_(this->coeffDict().lookup("log")),
548  outputFilePtr_(),
549  timeOld_(owner.mesh().time().value()),
550  hitFaceIDs_()
551 {
552  normal_ /= mag(normal_);
553 
554  word mode(this->coeffDict().lookup("mode"));
555  if (mode == "polygon")
556  {
557  List<Field<point>> polygons(this->coeffDict().lookup("polygons"));
558 
559  initPolygons(polygons);
560 
561  vector n0(this->coeffDict().lookup("normal"));
562  normal_ = vectorField(faces_.size(), n0);
563  }
564  else if (mode == "polygonWithNormal")
565  {
566  List<Tuple2<Field<point>, vector>> polygonAndNormal
567  (
568  this->coeffDict().lookup("polygons")
569  );
570 
571  List<Field<point>> polygons(polygonAndNormal.size());
572  normal_.setSize(polygonAndNormal.size());
573 
574  forAll(polygons, polyI)
575  {
576  polygons[polyI] = polygonAndNormal[polyI].first();
577  normal_[polyI] = polygonAndNormal[polyI].second();
578  normal_[polyI] /= mag(normal_[polyI]) + rootVSmall;
579  }
580 
581  initPolygons(polygons);
582  }
583  else if (mode == "concentricCircle")
584  {
585  vector n0(this->coeffDict().lookup("normal"));
586  normal_ = vectorField(1, n0);
587 
588  initConcentricCircles();
589  }
590  else
591  {
592  FatalIOErrorInFunction(this->coeffDict())
593  << "Unknown mode " << mode << ". Available options are "
594  << "polygon, polygonWithNormal and concentricCircle"
595  << exit(FatalIOError);
596  }
597 
598  mass_.setSize(faces_.size(), 0.0);
599  massTotal_.setSize(faces_.size(), 0.0);
600  massFlowRate_.setSize(faces_.size(), 0.0);
601 
602  makeLogFile(faces_, points_, area_);
603 }
604 
605 
606 template<class CloudType>
608 (
610 )
611 :
613  mode_(pc.mode_),
614  parcelType_(pc.parcelType_),
615  removeCollected_(pc.removeCollected_),
616  points_(pc.points_),
617  faces_(pc.faces_),
618  faceTris_(pc.faceTris_),
619  nSector_(pc.nSector_),
620  radius_(pc.radius_),
621  coordSys_(pc.coordSys_),
622  normal_(pc.normal_),
623  negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
624  surfaceFormat_(pc.surfaceFormat_),
625  resetOnWrite_(pc.resetOnWrite_),
626  totalTime_(pc.totalTime_),
627  mass_(pc.mass_),
628  massTotal_(pc.massTotal_),
629  massFlowRate_(pc.massFlowRate_),
630  log_(pc.log_),
631  outputFilePtr_(),
632  timeOld_(0.0),
633  hitFaceIDs_()
634 {}
635 
636 
637 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
638 
639 template<class CloudType>
641 {}
642 
643 
644 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
645 
646 template<class CloudType>
648 (
649  parcelType& p,
650  const scalar dt,
651  const point& position0,
652  bool& keepParticle
653 )
654 {
655  if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
656  {
657  return;
658  }
659 
660  hitFaceIDs_.clear();
661 
662  switch (mode_)
663  {
664  case mtPolygon:
665  {
666  collectParcelPolygon(position0, p.position());
667  break;
668  }
669  case mtConcentricCircle:
670  {
671  collectParcelConcentricCircles(position0, p.position());
672  break;
673  }
674  default:
675  {}
676  }
677 
678  forAll(hitFaceIDs_, i)
679  {
680  label facei = hitFaceIDs_[i];
681  scalar m = p.nParticle()*p.mass();
682 
683  if (negateParcelsOppositeNormal_)
684  {
685  scalar Unormal = 0;
686  vector Uhat = p.U();
687  switch (mode_)
688  {
689  case mtPolygon:
690  {
691  Unormal = Uhat & normal_[facei];
692  break;
693  }
694  case mtConcentricCircle:
695  {
696  Unormal = Uhat & normal_[0];
697  break;
698  }
699  default:
700  {}
701  }
702 
703  Uhat /= mag(Uhat) + rootVSmall;
704 
705  if (Unormal < 0)
706  {
707  m = -m;
708  }
709  }
710 
711  // Add mass contribution
712  mass_[facei] += m;
713 
714  if (nSector_ == 1)
715  {
716  mass_[facei + 1] += m;
717  mass_[facei + 2] += m;
718  mass_[facei + 3] += m;
719  }
720 
721  if (removeCollected_)
722  {
723  keepParticle = false;
724  }
725  }
726 }
727 
728 
729 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
type
Types of root.
Definition: Roots.H:52
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const char tab
Definition: Ostream.H:264
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
ParticleCollector(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
volVectorField vectorField(fieldObject, mesh)
Base class for graphics format writing. Entry points are.
Definition: writer.H:78
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
bool readBool(Istream &)
Definition: boolIO.C:60
void write()
Write post-processing info.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
T & first()
Return the first element of the list.
Definition: UListI.H:114
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
A class for handling words, derived from string.
Definition: word.H:59
mode_t mode(const fileName &, const bool followLink=true)
Return the file mode.
Definition: POSIX.C:459
const Type & value() const
Return const reference to value.
static const zero Zero
Definition: zero.H:97
label readLabel(Istream &is)
Definition: label.H:64
virtual void postMove(parcelType &p, const scalar dt, const point &position0, bool &keepParticle)
Post-move hook.
const scalar twoPi(2 *pi)
static const char nl
Definition: Ostream.H:265
virtual ~ParticleCollector()
Destructor.
labelList f(nPoints)
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
Definition: POSIX.C:289
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void setSize(const label)
Reset size of List.
Definition: List.C:281
vector point
Point is a vector.
Definition: point.H:41
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:331
Function object to collect the parcel mass- and mass flow rate over a set of polygons. The polygons can either be specified by sets of user- supplied points, or in a concentric circles arrangement. If a parcel is &#39;collected&#39;, it can be flagged to be removed from the domain using the removeCollected entry.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
messageStream Info
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
Templated cloud function object base class.
IOerror FatalIOError