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_;
161 nSector_ = this->coeffDict().template lookup<label>(
"nSector");
168 refDir = this->coeffDict().lookup(
"refDir");
169 refDir -= normal_[0]*(normal_[0] & refDir);
170 refDir /=
mag(refDir);
180 scalar dThetaSector = 360.0/scalar(nS);
181 label intervalPerSector =
max(1, ceil(dThetaSector/dTheta));
182 dTheta = dThetaSector/scalar(intervalPerSector);
184 label nPointPerSector = intervalPerSector + 1;
186 label nPointPerRadius = nS*(nPointPerSector - 1);
187 label nPoint = radius_.size()*nPointPerRadius;
188 label nFace = radius_.size()*nS;
193 points_.setSize(nPoint);
194 faces_.setSize(nFace);
195 area_.setSize(nFace);
197 coordSys_ = cylindricalCS(
"coordSys", origin, normal_[0], refDir,
false);
199 List<label> ptIDs(
identity(nPointPerRadius));
206 label pointOffset = radI*nPointPerRadius + 1;
208 for (
label i = 0; i < nPointPerRadius; i++)
210 label pI = i + pointOffset;
212 points_[pI] = coordSys_.globalPosition(pCyl);
217 DynamicList<label> facePts(2*nPointPerSector);
222 for (
label secI = 0; secI < nS; secI++)
229 for (
label ptI = 0; ptI < nPointPerSector; ptI++)
231 label i = ptI + secI*(nPointPerSector - 1);
232 label id = ptIDs.fcIndex(i - 1) + 1;
236 label facei = secI + radI*nS;
238 faces_[facei] = face(facePts);
239 area_[facei] = faces_[facei].mag(points_);
244 for (
label secI = 0; secI < nS; secI++)
248 label offset = (radI - 1)*nPointPerRadius + 1;
250 for (
label ptI = 0; ptI < nPointPerSector; ptI++)
252 label i = ptI + secI*(nPointPerSector - 1);
253 label id = offset + ptIDs.fcIndex(i - 1);
256 for (
label ptI = nPointPerSector-1; ptI >= 0; ptI--)
258 label i = ptI + secI*(nPointPerSector - 1);
259 label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
263 label facei = secI + radI*nS;
265 faces_[facei] = face(facePts);
266 area_[facei] = faces_[facei].mag(points_);
273 template<
class CloudType>
282 const label facePoint0 = faces_[facei][0];
284 const point& pf = points_[facePoint0];
286 const scalar d1 = normal_[facei] & (p1 - pf);
287 const scalar d2 = normal_[facei] & (p2 - pf);
296 const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
304 const face& f = faces_[facei];
305 const vector a = f.area(points_);
307 for (
label i = 0; i < f.size(); ++ i)
309 const label j = f.fcIndex(i);
310 const triPointRef t(pIntersect, points_[f[i]], points_[f[j]]);
311 if ((a & t.area()) < 0)
321 hitFaceIDs_.append(facei);
327 template<
class CloudType>
336 const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
337 const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
346 const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
350 if (r < radius_.last())
353 while (r > radius_[radI])
377 hitFaceIDs_.append(secI);
384 template<
class CloudType>
389 scalar timeNew = time.
value();
390 scalar timeElapsed = timeNew - timeOld_;
392 totalTime_ += timeElapsed;
394 const scalar
alpha = (totalTime_ - timeElapsed)/totalTime_;
395 const scalar beta = timeElapsed/totalTime_;
399 massFlowRate_[facei] =
400 alpha*massFlowRate_[facei] + beta*mass_[facei]/timeElapsed;
401 massTotal_[facei] += mass_[facei];
404 const label proci = Pstream::myProcNo();
409 this->getModelProperty(
"massTotal", faceMassTotal);
412 this->getModelProperty(
"massFlowRate", faceMassFlowRate);
415 scalar sumTotalMass = 0.0;
416 scalar sumAverageMFR = 0.0;
420 allProcMass[proci] = massTotal_[facei];
421 Pstream::gatherList(allProcMass);
422 faceMassTotal[facei] +=
sum(allProcMass);
424 scalarList allProcMassFlowRate(Pstream::nProcs());
425 allProcMassFlowRate[proci] = massFlowRate_[facei];
426 Pstream::gatherList(allProcMassFlowRate);
427 faceMassFlowRate[facei] +=
sum(allProcMassFlowRate);
429 sumTotalMass += faceMassTotal[facei];
430 sumAverageMFR += faceMassFlowRate[facei];
432 if (outputFilePtr_.valid())
437 <<
tab << faceMassTotal[facei]
438 <<
tab << faceMassFlowRate[facei]
443 Info<<
" sum(total mass) = " << sumTotalMass << nl
444 <<
" sum(average mass flow rate) = " << sumAverageMFR << nl
448 if (surfaceFormat_ !=
"none")
450 if (Pstream::master())
459 this->writeTimeDir(),
470 this->writeTimeDir(),
485 this->setModelProperty(
"massTotal", dummy);
486 this->setModelProperty(
"massFlowRate", dummy);
493 this->setModelProperty(
"massTotal", faceMassTotal);
494 this->setModelProperty(
"massFlowRate", faceMassFlowRate);
500 massTotal_[facei] = 0.0;
501 massFlowRate_[facei] = 0.0;
508 template<
class CloudType>
513 const word& modelName
518 parcelType_(this->coeffDict().lookupOrDefault(
"parcelType", -1)),
519 removeCollected_(this->coeffDict().
lookup(
"removeCollected")),
527 negateParcelsOppositeNormal_
531 surfaceFormat_(this->coeffDict().
lookup(
"surfaceFormat")),
532 resetOnWrite_(this->coeffDict().
lookup(
"resetOnWrite")),
537 log_(this->coeffDict().
lookup(
"log")),
539 timeOld_(owner.mesh().time().value()),
542 normal_ /=
mag(normal_);
545 if (mode ==
"polygon")
549 initPolygons(polygons);
554 else if (mode ==
"polygonWithNormal")
558 this->coeffDict().
lookup(
"polygons")
566 polygons[polyI] = polygonAndNormal[polyI].
first();
567 normal_[polyI] = polygonAndNormal[polyI].second();
568 normal_[polyI] /=
mag(normal_[polyI]) + rootVSmall;
571 initPolygons(polygons);
573 else if (mode ==
"concentricCircle")
578 initConcentricCircles();
583 <<
"Unknown mode " << mode <<
". Available options are " 584 <<
"polygon, polygonWithNormal and concentricCircle" 588 mass_.setSize(faces_.size(), 0.0);
589 massTotal_.setSize(faces_.size(), 0.0);
590 massFlowRate_.setSize(faces_.size(), 0.0);
592 makeLogFile(faces_, points_, area_);
596 template<
class CloudType>
604 parcelType_(pc.parcelType_),
605 removeCollected_(pc.removeCollected_),
608 faceTris_(pc.faceTris_),
609 nSector_(pc.nSector_),
611 coordSys_(pc.coordSys_),
613 negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
614 surfaceFormat_(pc.surfaceFormat_),
615 resetOnWrite_(pc.resetOnWrite_),
616 totalTime_(pc.totalTime_),
618 massTotal_(pc.massTotal_),
619 massFlowRate_(pc.massFlowRate_),
629 template<
class CloudType>
636 template<
class CloudType>
641 const point& position0,
645 if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
656 collectParcelPolygon(position0, p.position());
659 case mtConcentricCircle:
661 collectParcelConcentricCircles(position0, p.position());
670 label facei = hitFaceIDs_[i];
671 scalar m = p.nParticle()*p.mass();
673 if (negateParcelsOppositeNormal_)
681 Unormal = Uhat & normal_[facei];
684 case mtConcentricCircle:
686 Unormal = Uhat & normal_[0];
693 Uhat /=
mag(Uhat) + rootVSmall;
706 mass_[facei + 1] += m;
707 mass_[facei + 2] += m;
708 mass_[facei + 3] += m;
711 if (removeCollected_)
713 keepParticle =
false;
dimensionedScalar sign(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
mode_t mode(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file mode.
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)
stressControl lookup("compactNormalStress") >> compactNormalStress
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
A class for handling words, derived from string.
const Type & value() const
Return const reference to value.
virtual void postMove(parcelType &p, const scalar dt, const point &position0, bool &keepParticle)
Post-move hook.
const scalar twoPi(2 *pi)
IOstream::streamFormat writeFormat() const
Default write format.
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.
Vector< Cmpt > perpendicular(const Vector< Cmpt > &v)
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.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
triangle< point, const point & > triPointRef
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...
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
Templated base class for dsmc cloud.
Templated cloud function object base class.