36 template<
class CloudType>
40 const Field<point>& points,
41 const Field<scalar>& area
49 Info<<
"Creating output file" <<
endl;
52 if (Pstream::master())
55 mkDir(this->writeTimeDir());
60 new OFstream(this->writeTimeDir()/(
type() +
".dat"))
64 <<
"# Source : " <<
type() <<
nl 65 <<
"# Bins : " << faces.size() <<
nl 66 <<
"# Total area : " <<
sum(area) <<
nl;
69 <<
"# Geometry :" << nl
72 <<
tab <<
"(Centre_x Centre_y Centre_z)" 81 <<
tab << faces[i].centre(points)
88 <<
"# Output format:" <<
nl;
93 word binId =
"bin_" + id;
99 <<
tab <<
"mass[" <<
id <<
"]" 100 <<
tab <<
"massFlowRate[" <<
id <<
"]" 108 template<
class CloudType>
111 const List<Field<point>>& polygons
119 label np = polygons[polyI].size();
123 <<
"polygons must consist of at least 3 points" 130 label pointOffset = 0;
131 points_.setSize(nPoints);
132 faces_.setSize(polygons.size());
133 faceTris_.setSize(polygons.size());
134 area_.setSize(polygons.size());
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_);
142 DynamicList<face> tris;
143 f.triangles(points_, tris);
144 faceTris_[facei].transfer(tris);
146 faces_[facei].transfer(f);
148 pointOffset += polyPoints.size();
153 template<
class CloudType>
156 mode_ = mtConcentricCircle;
160 this->coeffDict().lookup(
"radius") >> radius_;
168 refDir = this->coeffDict().lookup(
"refDir");
169 refDir -= normal_[0]*(normal_[0] & refDir);
170 refDir /=
mag(refDir);
178 scalar magTangent = 0.0;
181 while (magTangent < small)
185 tangent = v - (v & normal_[0])*normal_[0];
186 magTangent =
mag(tangent);
189 refDir = tangent/magTangent;
193 scalar dThetaSector = 360.0/scalar(nS);
194 label intervalPerSector =
max(1, ceil(dThetaSector/dTheta));
195 dTheta = dThetaSector/scalar(intervalPerSector);
197 label nPointPerSector = intervalPerSector + 1;
199 label nPointPerRadius = nS*(nPointPerSector - 1);
200 label nPoint = radius_.size()*nPointPerRadius;
201 label nFace = radius_.size()*nS;
206 points_.setSize(nPoint);
207 faces_.setSize(nFace);
208 area_.setSize(nFace);
210 coordSys_ = cylindricalCS(
"coordSys", origin, normal_[0], refDir,
false);
212 List<label> ptIDs(
identity(nPointPerRadius));
219 label pointOffset = radI*nPointPerRadius + 1;
221 for (
label i = 0; i < nPointPerRadius; i++)
223 label pI = i + pointOffset;
225 points_[pI] = coordSys_.globalPosition(pCyl);
230 DynamicList<label> facePts(2*nPointPerSector);
235 for (
label secI = 0; secI < nS; secI++)
242 for (
label ptI = 0; ptI < nPointPerSector; ptI++)
244 label i = ptI + secI*(nPointPerSector - 1);
245 label id = ptIDs.fcIndex(i - 1) + 1;
249 label facei = secI + radI*nS;
251 faces_[facei] = face(facePts);
252 area_[facei] = faces_[facei].mag(points_);
257 for (
label secI = 0; secI < nS; secI++)
261 label offset = (radI - 1)*nPointPerRadius + 1;
263 for (
label ptI = 0; ptI < nPointPerSector; ptI++)
265 label i = ptI + secI*(nPointPerSector - 1);
266 label id = offset + ptIDs.fcIndex(i - 1);
269 for (
label ptI = nPointPerSector-1; ptI >= 0; ptI--)
271 label i = ptI + secI*(nPointPerSector - 1);
272 label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
276 label facei = secI + radI*nS;
278 faces_[facei] = face(facePts);
279 area_[facei] = faces_[facei].mag(points_);
286 template<
class CloudType>
295 const label facePoint0 = faces_[facei][0];
297 const point& pf = points_[facePoint0];
299 const scalar d1 = normal_[facei] & (p1 - pf);
300 const scalar d2 = normal_[facei] & (p2 - pf);
309 const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
317 const face& f = faces_[facei];
318 const vector a = f.area(points_);
320 for (
label i = 0; i < f.size(); ++ i)
322 const label j = f.fcIndex(i);
323 const triPointRef t(pIntersect, points_[f[i]], points_[f[j]]);
324 if ((a & t.area()) < 0)
334 hitFaceIDs_.append(facei);
340 template<
class CloudType>
349 const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
350 const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
359 const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
363 if (r < radius_.last())
366 while (r > radius_[radI])
390 hitFaceIDs_.append(secI);
397 template<
class CloudType>
402 scalar timeNew = time.
value();
403 scalar timeElapsed = timeNew - timeOld_;
405 totalTime_ += timeElapsed;
407 const scalar
alpha = (totalTime_ - timeElapsed)/totalTime_;
408 const scalar
beta = timeElapsed/totalTime_;
412 massFlowRate_[facei] =
413 alpha*massFlowRate_[facei] + beta*mass_[facei]/timeElapsed;
414 massTotal_[facei] += mass_[facei];
417 const label proci = Pstream::myProcNo();
422 this->getModelProperty(
"massTotal", faceMassTotal);
425 this->getModelProperty(
"massFlowRate", faceMassFlowRate);
428 scalar sumTotalMass = 0.0;
429 scalar sumAverageMFR = 0.0;
433 allProcMass[proci] = massTotal_[facei];
434 Pstream::gatherList(allProcMass);
435 faceMassTotal[facei] +=
sum(allProcMass);
437 scalarList allProcMassFlowRate(Pstream::nProcs());
438 allProcMassFlowRate[proci] = massFlowRate_[facei];
439 Pstream::gatherList(allProcMassFlowRate);
440 faceMassFlowRate[facei] +=
sum(allProcMassFlowRate);
442 sumTotalMass += faceMassTotal[facei];
443 sumAverageMFR += faceMassFlowRate[facei];
445 if (outputFilePtr_.valid())
450 <<
tab << faceMassTotal[facei]
451 <<
tab << faceMassFlowRate[facei]
456 Info<<
" sum(total mass) = " << sumTotalMass << nl
457 <<
" sum(average mass flow rate) = " << sumAverageMFR << nl
461 if (surfaceFormat_ !=
"none")
463 if (Pstream::master())
469 this->writeTimeDir(),
480 this->writeTimeDir(),
495 this->setModelProperty(
"massTotal", dummy);
496 this->setModelProperty(
"massFlowRate", dummy);
503 this->setModelProperty(
"massTotal", faceMassTotal);
504 this->setModelProperty(
"massFlowRate", faceMassFlowRate);
510 massTotal_[facei] = 0.0;
511 massFlowRate_[facei] = 0.0;
518 template<
class CloudType>
523 const word& modelName
528 parcelType_(this->coeffDict().lookupOrDefault(
"parcelType", -1)),
529 removeCollected_(this->coeffDict().
lookup(
"removeCollected")),
537 negateParcelsOppositeNormal_
541 surfaceFormat_(this->coeffDict().
lookup(
"surfaceFormat")),
542 resetOnWrite_(this->coeffDict().
lookup(
"resetOnWrite")),
547 log_(this->coeffDict().
lookup(
"log")),
549 timeOld_(owner.mesh().time().value()),
552 normal_ /=
mag(normal_);
555 if (mode ==
"polygon")
559 initPolygons(polygons);
564 else if (mode ==
"polygonWithNormal")
568 this->coeffDict().
lookup(
"polygons")
576 polygons[polyI] = polygonAndNormal[polyI].
first();
577 normal_[polyI] = polygonAndNormal[polyI].second();
578 normal_[polyI] /=
mag(normal_[polyI]) + rootVSmall;
581 initPolygons(polygons);
583 else if (mode ==
"concentricCircle")
588 initConcentricCircles();
593 <<
"Unknown mode " << mode <<
". Available options are " 594 <<
"polygon, polygonWithNormal and concentricCircle" 598 mass_.setSize(faces_.size(), 0.0);
599 massTotal_.setSize(faces_.size(), 0.0);
600 massFlowRate_.setSize(faces_.size(), 0.0);
602 makeLogFile(faces_, points_, area_);
606 template<
class CloudType>
614 parcelType_(pc.parcelType_),
615 removeCollected_(pc.removeCollected_),
618 faceTris_(pc.faceTris_),
619 nSector_(pc.nSector_),
621 coordSys_(pc.coordSys_),
623 negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
624 surfaceFormat_(pc.surfaceFormat_),
625 resetOnWrite_(pc.resetOnWrite_),
626 totalTime_(pc.totalTime_),
628 massTotal_(pc.massTotal_),
629 massFlowRate_(pc.massFlowRate_),
639 template<
class CloudType>
646 template<
class CloudType>
651 const point& position0,
655 if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
666 collectParcelPolygon(position0, p.position());
669 case mtConcentricCircle:
671 collectParcelConcentricCircles(position0, p.position());
680 label facei = hitFaceIDs_[i];
681 scalar m = p.nParticle()*p.mass();
683 if (negateParcelsOppositeNormal_)
691 Unormal = Uhat & normal_[facei];
694 case mtConcentricCircle:
696 Unormal = Uhat & normal_[0];
703 Uhat /=
mag(Uhat) + rootVSmall;
716 mass_[facei + 1] += m;
717 mass_[facei + 2] += m;
718 mass_[facei + 3] += m;
721 if (removeCollected_)
723 keepParticle =
false;
dimensionedScalar sign(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A list of keyword definitions, which are a keyword followed by any number of values (e...
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.
Ostream & endl(Ostream &os)
Add newline and flush stream.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
volVectorField vectorField(fieldObject, mesh)
Base class for graphics format writing. Entry points are.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Vector< scalar > vector
A scalar version of the templated Vector.
void write()
Write post-processing info.
const Time & time() const
Return the top-level database.
T & first()
Return the first element of the list.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
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
A class for handling words, derived from string.
mode_t mode(const fileName &, const bool followLink=true)
Return the file mode.
const Type & value() const
Return const reference to value.
label readLabel(Istream &is)
virtual void postMove(parcelType &p, const scalar dt, const point &position0, bool &keepParticle)
Post-move hook.
const scalar twoPi(2 *pi)
virtual ~ParticleCollector()
Destructor.
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
word name(const complex &)
Return a string representation of a complex.
void setSize(const label)
Reset size of List.
vector point
Point is a vector.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
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 'collected', it can be flagged to be removed from the domain using the removeCollected entry.
Mesh data needed to do the Finite Volume discretisation.
triangle< point, const point & > triPointRef
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...
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Templated base class for dsmc cloud.
Templated cloud function object base class.