39 template<
class CloudType>
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
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 <<
"]"
111 template<
class CloudType>
114 const List<Field<point>>& polygons
122 label np = polygons[polyI].size();
126 <<
"polygons must consist of at least 3 points"
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();
156 template<
class CloudType>
159 mode_ = mtConcentricCircle;
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_);
290 template<
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);
345 template<
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])
395 hitFaceIDs_.append(secI);
402 template<
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 const label proci = Pstream::myProcNo();
427 this->getModelProperty(
"massTotal", faceMassTotal);
430 this->getModelProperty(
"massFlowRate", faceMassFlowRate);
433 scalar sumTotalMass = 0.0;
434 scalar sumAverageMFR = 0.0;
438 allProcMass[proci] = massTotal_[facei];
439 Pstream::gatherList(allProcMass);
440 faceMassTotal[facei] +=
sum(allProcMass);
442 scalarList allProcMassFlowRate(Pstream::nProcs());
443 allProcMassFlowRate[proci] = massFlowRate_[facei];
444 Pstream::gatherList(allProcMassFlowRate);
445 faceMassFlowRate[facei] +=
sum(allProcMassFlowRate);
447 sumTotalMass += faceMassTotal[facei];
448 sumAverageMFR += faceMassFlowRate[facei];
450 if (outputFilePtr_.valid())
455 <<
tab << faceMassTotal[facei]
456 <<
tab << faceMassFlowRate[facei]
461 Info<<
" sum(total mass) = " << sumTotalMass <<
nl
462 <<
" sum(average mass flow rate) = " << sumAverageMFR <<
nl
466 if (surfaceFormat_ !=
"none" && Pstream::master())
471 this->coeffDict().subOrEmptyDict(
"formatOptions")
472 .subOrEmptyDict(surfaceFormat_)
484 (this->writeTimeDir() /
"collector"),
497 this->setModelProperty(
"massTotal", dummy);
498 this->setModelProperty(
"massFlowRate", dummy);
505 this->setModelProperty(
"massTotal", faceMassTotal);
506 this->setModelProperty(
"massFlowRate", faceMassFlowRate);
512 massTotal_[facei] = 0.0;
513 massFlowRate_[facei] = 0.0;
520 template<
class CloudType>
525 const word& modelName
530 parcelType_(this->coeffDict().getOrDefault(
"parcelType", -1)),
531 removeCollected_(this->coeffDict().getBool(
"removeCollected")),
532 resetOnWrite_(this->coeffDict().getBool(
"resetOnWrite")),
533 log_(this->coeffDict().getBool(
"log")),
541 negateParcelsOppositeNormal_
543 this->coeffDict().getBool(
"negateParcelsOppositeNormal")
545 surfaceFormat_(this->coeffDict().
lookup(
"surfaceFormat")),
554 normal_ /=
mag(normal_);
557 if (
mode ==
"polygon")
561 initPolygons(polygons);
566 else if (
mode ==
"polygonWithNormal")
570 this->coeffDict().
lookup(
"polygons")
574 normal_.
setSize(polygonAndNormal.size());
578 polygons[polyI] = polygonAndNormal[polyI].first();
579 normal_[polyI] =
normalised(polygonAndNormal[polyI].second());
582 initPolygons(polygons);
584 else if (
mode ==
"concentricCircle")
589 initConcentricCircles();
594 <<
"Unknown mode " <<
mode <<
". Available options are "
595 <<
"polygon, polygonWithNormal and concentricCircle"
599 mass_.setSize(faces_.size(), 0.0);
600 massTotal_.setSize(faces_.size(), 0.0);
601 massFlowRate_.setSize(faces_.size(), 0.0);
603 makeLogFile(faces_, points_, area_);
607 template<
class CloudType>
615 parcelType_(pc.parcelType_),
616 removeCollected_(pc.removeCollected_),
617 resetOnWrite_(pc.resetOnWrite_),
621 faceTris_(pc.faceTris_),
622 nSector_(pc.nSector_),
624 coordSys_(pc.coordSys_),
627 negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
628 surfaceFormat_(pc.surfaceFormat_),
629 totalTime_(pc.totalTime_),
631 massTotal_(pc.massTotal_),
632 massFlowRate_(pc.massFlowRate_),
641 template<
class CloudType>
646 const point& position0,
650 if ((parcelType_ != -1) && (parcelType_ !=
p.typeId()))
661 collectParcelPolygon(position0,
p.position());
664 case mtConcentricCircle:
666 collectParcelConcentricCircles(position0,
p.position());
675 label facei = hitFaceIDs_[i];
676 scalar m =
p.nParticle()*
p.mass();
678 if (negateParcelsOppositeNormal_)
686 Unormal = Uhat & normal_[facei];
689 case mtConcentricCircle:
691 Unormal = Uhat & normal_[0];
698 Uhat /=
mag(Uhat) + ROOTVSMALL;
711 mass_[facei + 1] += m;
712 mass_[facei + 2] += m;
713 mass_[facei + 3] += m;
716 if (removeCollected_)
718 keepParticle =
false;