38 template<
class CloudType>
42 const Field<point>&
points,
43 const Field<scalar>& area
51 Info<<
"Creating output file" <<
endl;
54 if (Pstream::master())
57 mkDir(this->writeTimeDir());
62 new OFstream(this->writeTimeDir()/(
type() +
".dat"))
66 <<
"# Source : " <<
type() <<
nl
67 <<
"# Bins : " << faces.size() <<
nl
68 <<
"# Total area : " <<
sum(area) <<
nl;
71 <<
"# Geometry :" <<
nl
74 <<
tab <<
"(Centre_x Centre_y Centre_z)"
90 <<
"# Output format:" <<
nl;
95 word binId =
"bin_" + id;
101 <<
tab <<
"mass[" <<
id <<
"]"
102 <<
tab <<
"massFlowRate[" <<
id <<
"]"
110 template<
class CloudType>
113 const List<Field<point>>& polygons
121 label np = polygons[polyI].size();
125 <<
"polygons must consist of at least 3 points"
132 label pointOffset = 0;
134 faces_.setSize(polygons.size());
135 area_.setSize(polygons.size());
138 const Field<point>& polyPoints = polygons[facei];
140 UIndirectList<point>(points_,
f) = polyPoints;
141 area_[facei] =
f.mag(points_);
142 faces_[facei].transfer(
f);
144 pointOffset += polyPoints.size();
149 template<
class CloudType>
152 mode_ = mtConcentricCircle;
154 vector origin(this->coeffDict().lookup(
"origin"));
156 this->coeffDict().lookup(
"radius") >> radius_;
157 nSector_ = this->coeffDict().template lookup<label>(
"nSector");
164 refDir = this->coeffDict().lookup(
"refDir");
165 refDir -= normal_[0]*(normal_[0] & refDir);
166 refDir /=
mag(refDir);
176 scalar dThetaSector = 360.0/scalar(nS);
177 label intervalPerSector =
max(1, ceil(dThetaSector/dTheta));
178 dTheta = dThetaSector/scalar(intervalPerSector);
180 label nPointPerSector = intervalPerSector + 1;
182 label nPointPerRadius = nS*(nPointPerSector - 1);
183 label nPoint = radius_.size()*nPointPerRadius;
184 label nFace = radius_.size()*nS;
189 points_.setSize(nPoint);
190 faces_.setSize(nFace);
191 area_.setSize(nFace);
193 coordSys_ = coordinateSystems::cylindrical
209 label pointOffset = radI*nPointPerRadius + 1;
211 for (
label i = 0; i < nPointPerRadius; i++)
213 label pI = i + pointOffset;
215 points_[pI] = coordSys_.globalPosition(pCyl);
220 DynamicList<label> facePts(2*nPointPerSector);
225 for (
label secI = 0; secI < nS; secI++)
232 for (
label ptI = 0; ptI < nPointPerSector; ptI++)
234 label i = ptI + secI*(nPointPerSector - 1);
235 label id = ptIDs.fcIndex(i - 1) + 1;
239 label facei = secI + radI*nS;
241 faces_[facei] = face(facePts);
242 area_[facei] = faces_[facei].mag(points_);
247 for (
label secI = 0; secI < nS; secI++)
253 for (
label ptI = 0; ptI < nPointPerSector; ptI++)
255 label i = ptI + secI*(nPointPerSector - 1);
259 for (
label ptI = nPointPerSector-1; ptI >= 0; ptI--)
261 label i = ptI + secI*(nPointPerSector - 1);
262 label id =
offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
266 label facei = secI + radI*nS;
268 faces_[facei] = face(facePts);
269 area_[facei] = faces_[facei].mag(points_);
276 template<
class CloudType>
285 const label facePoint0 = faces_[facei][0];
287 const point& pf = points_[facePoint0];
289 const scalar d1 = normal_[facei] & (p1 - pf);
290 const scalar d2 = normal_[facei] & (p2 - pf);
299 const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
307 const face&
f = faces_[facei];
308 const vector a =
f.area(points_);
310 for (
label i = 0; i <
f.size(); ++ i)
312 const label j =
f.fcIndex(i);
313 const triPointRef t(pIntersect, points_[
f[i]], points_[
f[j]]);
314 if ((a & t.area()) < 0)
324 hitFaceIDs_.append(facei);
330 template<
class CloudType>
339 const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
340 const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
349 const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
353 if (r < radius_.last())
356 while (r > radius_[radI])
380 hitFaceIDs_.append(secI);
387 template<
class CloudType>
390 const fvMesh& mesh = this->owner().mesh();
392 scalar timeNew = time.
value();
393 scalar timeElapsed = timeNew - timeOld_;
395 totalTime_ += timeElapsed;
397 const scalar
alpha = (totalTime_ - timeElapsed)/totalTime_;
398 const scalar beta = timeElapsed/totalTime_;
402 massFlowRate_[facei] =
403 alpha*massFlowRate_[facei] + beta*mass_[facei]/timeElapsed;
404 massTotal_[facei] += mass_[facei];
407 const label proci = Pstream::myProcNo();
412 this->getModelProperty(
"massTotal", faceMassTotal);
415 this->getModelProperty(
"massFlowRate", faceMassFlowRate);
418 scalar sumTotalMass = 0.0;
419 scalar sumAverageMFR = 0.0;
423 allProcMass[proci] = massTotal_[facei];
424 Pstream::gatherList(allProcMass);
425 faceMassTotal[facei] +=
sum(allProcMass);
427 scalarList allProcMassFlowRate(Pstream::nProcs());
428 allProcMassFlowRate[proci] = massFlowRate_[facei];
429 Pstream::gatherList(allProcMassFlowRate);
430 faceMassFlowRate[facei] +=
sum(allProcMassFlowRate);
432 sumTotalMass += faceMassTotal[facei];
433 sumAverageMFR += faceMassFlowRate[facei];
435 if (outputFilePtr_.valid())
440 <<
tab << faceMassTotal[facei]
441 <<
tab << faceMassFlowRate[facei]
446 Info<<
" sum(total mass) = " << sumTotalMass <<
nl
447 <<
" sum(average mass flow rate) = " << sumAverageMFR <<
nl
451 if (surfaceFormat_ !=
"none")
453 if (Pstream::master())
462 this->writeTimeDir(),
479 this->setModelProperty(
"massTotal", dummy);
480 this->setModelProperty(
"massFlowRate", dummy);
487 this->setModelProperty(
"massTotal", faceMassTotal);
488 this->setModelProperty(
"massFlowRate", faceMassFlowRate);
494 massTotal_[facei] = 0.0;
495 massFlowRate_[facei] = 0.0;
502 template<
class CloudType>
507 const word& modelName
512 parcelType_(this->coeffDict().lookupOrDefault(
"parcelType", -1)),
513 removeCollected_(this->coeffDict().lookup(
"removeCollected")),
520 negateParcelsOppositeNormal_
522 readBool(this->coeffDict().lookup(
"negateParcelsOppositeNormal"))
524 surfaceFormat_(this->coeffDict().lookup(
"surfaceFormat")),
525 resetOnWrite_(this->coeffDict().lookup(
"resetOnWrite")),
530 log_(this->coeffDict().lookup(
"log")),
532 timeOld_(owner.mesh().time().value()),
535 normal_ /=
mag(normal_);
538 if (
mode ==
"polygon")
542 initPolygons(polygons);
547 else if (
mode ==
"polygonWithNormal")
559 polygons[polyI] = polygonAndNormal[polyI].
first();
560 normal_[polyI] = polygonAndNormal[polyI].second();
561 normal_[polyI] /=
mag(normal_[polyI]) + rootVSmall;
564 initPolygons(polygons);
566 else if (
mode ==
"concentricCircle")
571 initConcentricCircles();
576 <<
"Unknown mode " << mode <<
". Available options are "
577 <<
"polygon, polygonWithNormal and concentricCircle"
585 makeLogFile(faces_, points_, area_);
589 template<
class CloudType>
597 parcelType_(pc.parcelType_),
598 removeCollected_(pc.removeCollected_),
601 nSector_(pc.nSector_),
603 coordSys_(pc.coordSys_),
605 negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
606 surfaceFormat_(pc.surfaceFormat_),
607 resetOnWrite_(pc.resetOnWrite_),
608 totalTime_(pc.totalTime_),
610 massTotal_(pc.massTotal_),
611 massFlowRate_(pc.massFlowRate_),
621 template<
class CloudType>
628 template<
class CloudType>
633 const point& position0,
637 if ((parcelType_ != -1) && (parcelType_ !=
p.typeId()))
651 p.position(this->owner().mesh())
655 case mtConcentricCircle:
657 collectParcelConcentricCircles
660 p.position(this->owner().mesh())
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;
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
#define forAll(list, i)
Loop across all elements in list.
Templated cloud function object base class.
Templated base class for dsmc cloud.
void size(const label)
Override size to be inconsistent with allocated storage.
void setSize(const label)
Reset size of List.
Function object to collect the parcel mass- and mass flow rate over a set of polygons....
virtual void postMove(parcelType &p, const scalar dt, const point &position0, bool &keepParticle)
Post-move hook.
virtual ~ParticleCollector()
Destructor.
ParticleCollector(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
void write()
Write post-processing info.
Templated 3D SphericalTensor derived from VectorSpace adding construction from 1 component,...
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
T & first()
Return the first element of the list.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
A coordinate rotation specified using global axis.
A list of keyword definitions, which are a keyword followed by any number of values (e....
const Type & value() const
Return const reference to value.
const word & name() const
Return const reference to name.
Mesh data needed to do the Finite Volume discretisation.
const Time & time() const
Return the top-level database.
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
A class for handling words, derived from string.
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
volScalarField alpha(IOobject("alpha", runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
const scalar twoPi(2 *pi)
errorManipArg< error, int > exit(error &err, const int errNo=1)
mode_t mode(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file mode.
dimensionedScalar sign(const dimensionedScalar &ds)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & endl(Ostream &os)
Add newline and flush stream.
word name(const bool)
Return a word representation of a bool.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
static const Identity< scalar > I
vector point
Point is a vector.
Vector< scalar > vector
A scalar version of the templated Vector.
Vector< Cmpt > perpendicular(const Vector< Cmpt > &v)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Field< vector > vectorField
Specialisation of Field<T> for vector.
triangle< point, const point & > triPointRef
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
void offset(label &lst, const label o)
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Unit conversion functions.