39 template<
class CloudType>
43 const Field<point>&
points,
44 const Field<scalar>&
area 52 if (Pstream::master())
55 mkDir(this->writeTimeDir());
60 new OFstream(this->writeTimeDir()/(
type() +
".dat"))
64 <<
"# Source : " <<
type() <<
nl 65 <<
"# Bins : " << faces.size() <<
nl 69 <<
"# Geometry :" <<
nl 72 <<
tab <<
"(Centre_x Centre_y Centre_z)" 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;
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;
158 vector origin(this->coeffDict().lookup(
"origin"));
160 this->coeffDict().readEntry(
"radius", radius_);
161 this->coeffDict().readEntry(
"nSector", nSector_);
168 this->coeffDict().readEntry(
"refDir", refDir);
170 refDir -= normal_[0]*(normal_[0] & refDir);
179 scalar magTangent = 0.0;
181 Random& rnd = this->owner().rndGen();
182 while (magTangent < SMALL)
186 tangent = v - (v & normal_[0])*normal_[0];
187 magTangent =
mag(tangent);
190 refDir = tangent/magTangent;
194 scalar dThetaSector = 360.0/scalar(nS);
195 label intervalPerSector =
max(1, ceil(dThetaSector/dTheta));
196 dTheta = dThetaSector/scalar(intervalPerSector);
198 label nPointPerSector = intervalPerSector + 1;
200 label nPointPerRadius = nS*(nPointPerSector - 1);
201 label nPoint = radius_.size()*nPointPerRadius;
202 label nFace = radius_.size()*nS;
207 points_.setSize(nPoint);
208 faces_.setSize(nFace);
209 area_.setSize(nFace);
211 coordSys_ = coordSystem::cylindrical(origin, normal_[0], refDir);
213 List<label> ptIDs(
identity(nPointPerRadius));
220 label pointOffset = radI*nPointPerRadius + 1;
222 for (label i = 0; i < nPointPerRadius; i++)
224 label pI = i + pointOffset;
226 points_[pI] = coordSys_.globalPosition(pCyl);
231 DynamicList<label> facePts(2*nPointPerSector);
236 for (label secI = 0; secI < nS; secI++)
243 for (label ptI = 0; ptI < nPointPerSector; ptI++)
245 label i = ptI + secI*(nPointPerSector - 1);
246 label
id = ptIDs.fcIndex(i - 1) + 1;
250 label facei = secI + radI*nS;
252 faces_[facei] = face(facePts);
253 area_[facei] = faces_[facei].mag(points_);
258 for (label secI = 0; secI < nS; secI++)
262 label offset = (radI - 1)*nPointPerRadius + 1;
264 for (label ptI = 0; ptI < nPointPerSector; ptI++)
266 label i = ptI + secI*(nPointPerSector - 1);
267 label
id = offset + ptIDs.fcIndex(i - 1);
270 for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
272 label i = ptI + secI*(nPointPerSector - 1);
273 label
id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
277 label facei = secI + radI*nS;
279 faces_[facei] = face(facePts);
280 area_[facei] = faces_[facei].mag(points_);
287 template<
class CloudType>
296 const label facePoint0 = faces_[facei][0];
298 const point& pf = points_[facePoint0];
300 const scalar d1 = normal_[facei] & (p1 - pf);
301 const scalar d2 = normal_[facei] & (p2 - pf);
310 const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
318 const face&
f = faces_[facei];
319 const vector areaNorm =
f.areaNormal(points_);
321 for (label i = 0; i <
f.size(); ++i)
323 const label j =
f.fcIndex(i);
324 const triPointRef t(pIntersect, points_[
f[i]], points_[
f[j]]);
326 if ((areaNorm & t.areaNormal()) < 0)
336 hitFaceIDs_.append(facei);
342 template<
class CloudType>
351 const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
352 const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
361 const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
365 if (r < radius_.last())
368 while (r > radius_[radI])
392 hitFaceIDs_.append(secI);
399 template<
class CloudType>
402 const fvMesh&
mesh = this->owner().mesh();
403 const Time& time =
mesh.time();
404 scalar timeNew = time.value();
405 scalar timeElapsed = timeNew - timeOld_;
407 totalTime_ += timeElapsed;
409 const scalar
alpha = (totalTime_ - timeElapsed)/totalTime_;
410 const scalar
beta = timeElapsed/totalTime_;
414 massFlowRate_[facei] =
415 alpha*massFlowRate_[facei] +
beta*mass_[facei]/timeElapsed;
416 massTotal_[facei] += mass_[facei];
421 Field<scalar> faceMassTotal(mass_.size(),
Zero);
422 this->getModelProperty(
"massTotal", faceMassTotal);
424 Field<scalar> faceMassFlowRate(massFlowRate_.size(),
Zero);
425 this->getModelProperty(
"massFlowRate", faceMassFlowRate);
428 scalar sumTotalMass = 0.0;
429 scalar sumAverageMFR = 0.0;
432 faceMassTotal[facei] +=
435 faceMassFlowRate[facei] +=
438 sumTotalMass += faceMassTotal[facei];
439 sumAverageMFR += faceMassFlowRate[facei];
446 <<
tab << faceMassTotal[facei]
447 <<
tab << faceMassFlowRate[facei]
452 Log_<<
" sum(total mass) = " << sumTotalMass <<
nl 453 <<
" sum(average mass flow rate) = " << sumAverageMFR <<
nl 457 if (Pstream::master() && surfaceFormat_ !=
"none")
474 (this->writeTimeDir() /
"collector"),
479 writer->write(
"massFlowRate", faceMassFlowRate);
480 writer->write(
"massTotal", faceMassTotal);
486 Field<scalar> dummy(faceMassTotal.size(),
Zero);
487 this->setModelProperty(
"massTotal", dummy);
488 this->setModelProperty(
"massFlowRate", dummy);
495 this->setModelProperty(
"massTotal", faceMassTotal);
496 this->setModelProperty(
"massFlowRate", faceMassFlowRate);
502 massTotal_[facei] = 0.0;
503 massFlowRate_[facei] = 0.0;
510 template<
class CloudType>
515 const word& modelName
520 parcelType_(this->coeffDict().getOrDefault(
"parcelType", -1)),
521 removeCollected_(this->coeffDict().getBool(
"removeCollected")),
522 resetOnWrite_(this->coeffDict().getBool(
"resetOnWrite")),
523 logToFile_(this->coeffDict().getBool(
"log")),
531 negateParcelsOppositeNormal_
533 this->coeffDict().getBool(
"negateParcelsOppositeNormal")
535 surfaceFormat_(this->coeffDict().getWord(
"surfaceFormat")),
541 timeOld_(owner.
mesh().time().value()),
544 normal_ /=
mag(normal_);
547 if (
mode ==
"polygon")
551 initPolygons(polygons);
556 else if (
mode ==
"polygonWithNormal")
558 List<Tuple2<Field<point>,
vector>> polygonAndNormal
563 List<Field<point>> polygons(polygonAndNormal.size());
564 normal_.
setSize(polygonAndNormal.size());
568 polygons[polyI] = polygonAndNormal[polyI].first();
569 normal_[polyI] =
normalised(polygonAndNormal[polyI].second());
572 initPolygons(polygons);
574 else if (
mode ==
"concentricCircle")
579 initConcentricCircles();
584 <<
"Unknown mode " << mode <<
". Available options are " 585 <<
"polygon, polygonWithNormal and concentricCircle" 593 makeLogFile(faces_, points_, area_);
597 template<
class CloudType>
605 parcelType_(pc.parcelType_),
606 removeCollected_(pc.removeCollected_),
607 resetOnWrite_(pc.resetOnWrite_),
608 logToFile_(pc.logToFile_),
611 faceTris_(pc.faceTris_),
612 nSector_(pc.nSector_),
614 coordSys_(pc.coordSys_),
617 negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
618 surfaceFormat_(pc.surfaceFormat_),
619 totalTime_(pc.totalTime_),
621 massTotal_(pc.massTotal_),
622 massFlowRate_(pc.massFlowRate_),
631 template<
class CloudType>
636 const point& position0,
637 const typename parcelType::trackingData&
td 640 bool keepParticle =
true;
642 if ((parcelType_ != -1) && (parcelType_ !=
p.typeId()))
653 collectParcelPolygon(position0,
p.position());
656 case mtConcentricCircle:
658 collectParcelConcentricCircles(position0,
p.position());
667 label facei = hitFaceIDs_[i];
668 scalar m =
p.nParticle()*
p.mass();
670 if (negateParcelsOppositeNormal_)
678 Unormal = Uhat & normal_[facei];
681 case mtConcentricCircle:
683 Unormal = Uhat & normal_[0];
690 Uhat /=
mag(Uhat) + ROOTVSMALL;
703 mass_[facei + 1] += m;
704 mass_[facei + 2] += m;
705 mass_[facei + 3] += m;
708 if (removeCollected_)
710 keepParticle =
false;
dimensionedScalar sign(const dimensionedScalar &ds)
void size(const label n)
Older name for setAddressableSize.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
vtk::lineWriter writer(edgeCentres, edgeList::null(), fileName(aMesh.time().globalPath()/"finiteArea-edgesCentres"))
A list of keyword definitions, which are a keyword followed by a number of values (eg...
virtual bool postMove(parcelType &p, const scalar dt, const point &position0, const typename parcelType::trackingData &td)
Post-move hook.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
ParticleCollector(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Unit conversion functions.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
constexpr char tab
The tab '\t' character(0x09)
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
void write()
Write post-processing info.
Lookup type of boundary radiation properties.
#define forAll(list, i)
Loop across all elements in list.
List< face > faceList
List of faces.
void setSize(const label n)
Alias for resize()
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
#define Log_
Report write to Foam::Info if the class log switch is true.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
constexpr scalar twoPi(2 *M_PI)
const wordList area
Standard area field types (scalar, vector, tensor, etc)
A class for handling words, derived from Foam::string.
constexpr scalar pi(M_PI)
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
#define DebugInfo
Report an information message using Foam::Info.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
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.
return returnReduce(nRefine-oldNRefine, sumOp< label >())
triangle< point, const point & > triPointRef
A triangle using referred points.
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
Templated base class for dsmc cloud.
Templated cloud function object base class.
const dictionary formatOptions(propsDict.subOrEmptyDict("formatOptions", keyType::LITERAL))
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
static constexpr const zero Zero
Global zero (0)