39template<
class CloudType>
42 const faceList& faces,
43 const Field<point>&
points,
44 const Field<scalar>& area
52 Info<<
"Creating output file" <<
endl;
55 if (Pstream::master())
58 mkDir(this->writeTimeDir());
63 new OFstream(this->writeTimeDir()/(
type() +
".dat"))
67 <<
"# Source : " <<
type() <<
nl
68 <<
"# Bins : " << faces.size() <<
nl
69 <<
"# Total area : " <<
sum(area) <<
nl;
72 <<
"# Geometry :" <<
nl
75 <<
tab <<
"(Centre_x Centre_y Centre_z)"
91 <<
"# Output format:" <<
nl;
96 word binId =
"bin_" + id;
102 <<
tab <<
"mass[" <<
id <<
"]"
103 <<
tab <<
"massFlowRate[" <<
id <<
"]"
111template<
class CloudType>
114 const List<Field<point>>& polygons
122 label np = polygons[polyI].size();
126 <<
"polygons must consist of at least 3 points"
127 <<
exit(FatalIOError);
133 label pointOffset = 0;
135 faces_.setSize(polygons.size());
136 faceTris_.setSize(polygons.size());
137 area_.setSize(polygons.size());
140 const Field<point>& polyPoints = polygons[facei];
141 face
f(
identity(polyPoints.size(), pointOffset));
142 UIndirectList<point>(points_,
f) = polyPoints;
143 area_[facei] =
f.mag(points_);
145 DynamicList<face> tris;
146 f.triangles(points_, tris);
147 faceTris_[facei].transfer(tris);
149 faces_[facei].transfer(
f);
151 pointOffset += polyPoints.size();
156template<
class CloudType>
159 mode_ = mtConcentricCircle;
161 vector origin(this->coeffDict().lookup(
"origin"));
163 this->coeffDict().readEntry(
"radius", radius_);
164 this->coeffDict().readEntry(
"nSector", nSector_);
171 this->coeffDict().readEntry(
"refDir", refDir);
173 refDir -= normal_[0]*(normal_[0] & refDir);
182 scalar magTangent = 0.0;
184 Random& rnd = this->owner().rndGen();
185 while (magTangent < SMALL)
189 tangent = v - (v & normal_[0])*normal_[0];
190 magTangent =
mag(tangent);
193 refDir = tangent/magTangent;
197 scalar dThetaSector = 360.0/scalar(nS);
198 label intervalPerSector =
max(1, ceil(dThetaSector/dTheta));
199 dTheta = dThetaSector/scalar(intervalPerSector);
201 label nPointPerSector = intervalPerSector + 1;
203 label nPointPerRadius = nS*(nPointPerSector - 1);
204 label nPoint = radius_.size()*nPointPerRadius;
205 label nFace = radius_.size()*nS;
210 points_.setSize(nPoint);
211 faces_.setSize(nFace);
212 area_.setSize(nFace);
214 coordSys_ = coordSystem::cylindrical(origin, normal_[0], refDir);
216 List<label> ptIDs(
identity(nPointPerRadius));
223 label pointOffset = radI*nPointPerRadius + 1;
225 for (label i = 0; i < nPointPerRadius; i++)
227 label pI = i + pointOffset;
229 points_[pI] = coordSys_.globalPosition(pCyl);
234 DynamicList<label> facePts(2*nPointPerSector);
239 for (label secI = 0; secI < nS; secI++)
246 for (label ptI = 0; ptI < nPointPerSector; ptI++)
248 label i = ptI + secI*(nPointPerSector - 1);
249 label
id = ptIDs.fcIndex(i - 1) + 1;
253 label facei = secI + radI*nS;
255 faces_[facei] = face(facePts);
256 area_[facei] = faces_[facei].mag(points_);
261 for (label secI = 0; secI < nS; secI++)
265 label offset = (radI - 1)*nPointPerRadius + 1;
267 for (label ptI = 0; ptI < nPointPerSector; ptI++)
269 label i = ptI + secI*(nPointPerSector - 1);
270 label
id = offset + ptIDs.fcIndex(i - 1);
273 for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
275 label i = ptI + secI*(nPointPerSector - 1);
276 label
id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
280 label facei = secI + radI*nS;
282 faces_[facei] = face(facePts);
283 area_[facei] = faces_[facei].mag(points_);
290template<
class CloudType>
299 const label facePoint0 = faces_[facei][0];
301 const point& pf = points_[facePoint0];
303 const scalar d1 = normal_[facei] & (p1 - pf);
304 const scalar d2 = normal_[facei] & (p2 - pf);
313 const point pIntersect = p1 + (d1/(d1 - d2))*(p2 - p1);
321 const face&
f = faces_[facei];
322 const vector areaNorm =
f.areaNormal(points_);
324 for (label i = 0; i <
f.size(); ++i)
326 const label j =
f.fcIndex(i);
327 const triPointRef t(pIntersect, points_[
f[i]], points_[
f[j]]);
329 if ((areaNorm & t.areaNormal()) < 0)
339 hitFaceIDs_.append(facei);
345template<
class CloudType>
354 const scalar d1 = normal_[0] & (p1 - coordSys_.origin());
355 const scalar d2 = normal_[0] & (p2 - coordSys_.origin());
364 const point pCyl = coordSys_.localPosition(p1 + (d1/(d1 - d2))*(p2 - p1));
368 if (r < radius_.last())
371 while (r > radius_[radI])
382 scalar theta = pCyl[1] + constant::mathematical::pi;
388 scalar(nSector_)*theta/constant::mathematical::twoPi
395 hitFaceIDs_.append(secI);
402template<
class CloudType>
407 scalar timeNew = time.
value();
408 scalar timeElapsed = timeNew - timeOld_;
410 totalTime_ += timeElapsed;
412 const scalar
alpha = (totalTime_ - timeElapsed)/totalTime_;
413 const scalar
beta = timeElapsed/totalTime_;
417 massFlowRate_[facei] =
418 alpha*massFlowRate_[facei] +
beta*mass_[facei]/timeElapsed;
419 massTotal_[facei] += mass_[facei];
422 Info<< type() <<
" output:" <<
nl;
425 this->getModelProperty(
"massTotal", faceMassTotal);
428 this->getModelProperty(
"massFlowRate", faceMassFlowRate);
431 scalar sumTotalMass = 0.0;
432 scalar sumAverageMFR = 0.0;
435 faceMassTotal[facei] +=
438 faceMassFlowRate[facei] +=
441 sumTotalMass += faceMassTotal[facei];
442 sumAverageMFR += faceMassFlowRate[facei];
449 <<
tab << faceMassTotal[facei]
450 <<
tab << faceMassFlowRate[facei]
455 Info<<
" sum(total mass) = " << sumTotalMass <<
nl
456 <<
" sum(average mass flow rate) = " << sumAverageMFR <<
nl
460 if (surfaceFormat_ !=
"none" && Pstream::master())
462 auto writer = surfaceWriter::New
465 this->coeffDict().subOrEmptyDict(
"formatOptions")
466 .subOrEmptyDict(surfaceFormat_)
478 (this->writeTimeDir() /
"collector"),
483 writer->write(
"massFlowRate", faceMassFlowRate);
484 writer->write(
"massTotal", faceMassTotal);
491 this->setModelProperty(
"massTotal", dummy);
492 this->setModelProperty(
"massFlowRate", dummy);
499 this->setModelProperty(
"massTotal", faceMassTotal);
500 this->setModelProperty(
"massFlowRate", faceMassFlowRate);
506 massTotal_[facei] = 0.0;
507 massFlowRate_[facei] = 0.0;
514template<
class CloudType>
519 const word& modelName
524 parcelType_(this->coeffDict().getOrDefault(
"parcelType", -1)),
525 removeCollected_(this->coeffDict().getBool(
"removeCollected")),
526 resetOnWrite_(this->coeffDict().getBool(
"resetOnWrite")),
527 log_(this->coeffDict().getBool(
"log")),
535 negateParcelsOppositeNormal_
537 this->coeffDict().getBool(
"negateParcelsOppositeNormal")
539 surfaceFormat_(this->coeffDict().
lookup(
"surfaceFormat")),
545 timeOld_(owner.
mesh().time().value()),
548 normal_ /=
mag(normal_);
551 if (
mode ==
"polygon")
555 initPolygons(polygons);
560 else if (
mode ==
"polygonWithNormal")
572 polygons[polyI] = polygonAndNormal[polyI].
first();
573 normal_[polyI] =
normalised(polygonAndNormal[polyI].second());
576 initPolygons(polygons);
578 else if (
mode ==
"concentricCircle")
583 initConcentricCircles();
588 <<
"Unknown mode " << mode <<
". Available options are "
589 <<
"polygon, polygonWithNormal and concentricCircle"
597 makeLogFile(faces_, points_, area_);
601template<
class CloudType>
609 parcelType_(pc.parcelType_),
610 removeCollected_(pc.removeCollected_),
611 resetOnWrite_(pc.resetOnWrite_),
615 faceTris_(pc.faceTris_),
616 nSector_(pc.nSector_),
618 coordSys_(pc.coordSys_),
621 negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
622 surfaceFormat_(pc.surfaceFormat_),
623 totalTime_(pc.totalTime_),
625 massTotal_(pc.massTotal_),
626 massFlowRate_(pc.massFlowRate_),
635template<
class CloudType>
640 const point& position0,
644 if ((parcelType_ != -1) && (parcelType_ !=
p.typeId()))
655 collectParcelPolygon(position0,
p.position());
658 case mtConcentricCircle:
660 collectParcelConcentricCircles(position0,
p.position());
669 label facei = hitFaceIDs_[i];
670 scalar m =
p.nParticle()*
p.mass();
672 if (negateParcelsOppositeNormal_)
680 Unormal = Uhat & normal_[facei];
683 case mtConcentricCircle:
685 Unormal = Uhat & normal_[0];
692 Uhat /=
mag(Uhat) + ROOTVSMALL;
705 mass_[facei + 1] += m;
706 mass_[facei + 2] += m;
707 mass_[facei + 3] += m;
710 if (removeCollected_)
712 keepParticle =
false;
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
Templated cloud function object base class.
Templated base class for dsmc cloud.
Generic templated field type.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void setSize(const label n)
Alias for resize()
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.
void write()
Write post-processing info.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
static word timeName(const scalar t, const int precision=precision_)
T & first()
Return the first element of the list.
void size(const label n)
Older name for setAddressableSize.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
const Type & value() const
Return const reference to value.
Mesh data needed to do the Finite Volume discretisation.
Lookup type of boundary radiation properties.
const dictionary & coeffDict() const
Return const access to the coefficients dictionary.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
const wordList area
Standard area field types (scalar, vector, tensor, etc)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
dimensionedScalar sign(const dimensionedScalar &ds)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
vector point
Point is a vector.
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
triangle< point, const point & > triPointRef
A triangle using referred points.
Field< vector > vectorField
Specialisation of Field<T> for vector.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
constexpr char tab
The tab '\t' character(0x09)
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
#define forAll(list, i)
Loop across all elements in list.
Unit conversion functions.