37 template<
class CloudType>
43 const fvMesh&
mesh = this->owner().mesh();
45 const label nCells = cellZoneCells.size();
46 Random& rnd = this->owner().rndGen();
48 DynamicList<vector> positions(nCells);
49 DynamicList<label> injectorCells(nCells);
50 DynamicList<label> injectorTetFaces(nCells);
51 DynamicList<label> injectorTetPts(nCells);
53 scalar newParticlesTotal = 0.0;
54 label addParticlesTotal = 0;
58 const label celli = cellZoneCells[i];
61 const scalar newParticles = V[celli]*numberDensity_;
62 newParticlesTotal += newParticles;
63 label addParticles = floor(newParticles);
64 addParticlesTotal += addParticles;
66 const scalar
diff = newParticlesTotal - addParticlesTotal;
69 label corr = floor(
diff);
71 addParticlesTotal += corr;
75 const List<tetIndices> cellTetIs =
76 polyMeshTetDecomposition::cellTetIndices(
mesh, celli);
80 for (label tetI = 1; tetI < cellTetIs.size() - 1; tetI++)
83 cTetVFrac[tetI-1] + cellTetIs[tetI].tet(
mesh).mag()/V[celli];
85 cTetVFrac.
last() = 1.0;
88 for (label pI = 0; pI < addParticles; pI++)
90 const scalar volFrac = rnd.sample01<scalar>();
94 if (cTetVFrac[vfI] > volFrac)
100 positions.append(cellTetIs[tetI].tet(
mesh).randomPoint(rnd));
102 injectorCells.append(celli);
103 injectorTetFaces.append(cellTetIs[tetI].face());
104 injectorTetPts.append(cellTetIs[tetI].tetPt());
109 globalIndex globalPositions(positions.size());
110 List<vector> allPositions(globalPositions.totalSize(),
point::max);
111 List<label> allInjectorCells(globalPositions.totalSize(), -1);
112 List<label> allInjectorTetFaces(globalPositions.totalSize(), -1);
113 List<label> allInjectorTetPts(globalPositions.totalSize(), -1);
119 globalPositions.localSize(Pstream::myProcNo()),
120 globalPositions.localStart(Pstream::myProcNo())
123 Pstream::listCombineReduce(allPositions, minEqOp<point>());
129 globalPositions.localSize(Pstream::myProcNo()),
130 globalPositions.localStart(Pstream::myProcNo())
135 globalPositions.localSize(Pstream::myProcNo()),
136 globalPositions.localStart(Pstream::myProcNo())
137 ) = injectorTetFaces;
141 globalPositions.localSize(Pstream::myProcNo()),
142 globalPositions.localStart(Pstream::myProcNo())
146 positions_.transfer(allPositions);
147 injectorCells_.transfer(allInjectorCells);
148 injectorTetFaces_.transfer(allInjectorTetFaces);
149 injectorTetPts_.transfer(allInjectorTetPts);
154 OFstream
points(
"points.obj");
165 template<
class CloudType>
170 const word& modelName
174 cellZoneName_(this->coeffDict().
lookup(
"cellZone")),
175 numberDensity_(this->coeffDict().getScalar(
"numberDensity")),
181 U0_(this->coeffDict().
lookup(
"U0")),
186 this->coeffDict().subDict(
"sizeDistribution"), owner.
rndGen()
194 template<
class CloudType>
201 cellZoneName_(im.cellZoneName_),
202 numberDensity_(im.numberDensity_),
203 positions_(im.positions_),
204 injectorCells_(im.injectorCells_),
205 injectorTetFaces_(im.injectorTetFaces_),
206 injectorTetPts_(im.injectorTetPts_),
207 diameters_(im.diameters_),
209 sizeDistribution_(im.sizeDistribution_.clone())
215 template<
class CloudType>
222 template<
class CloudType>
227 const label zoneI =
mesh.cellZones().findZoneID(cellZoneName_);
232 <<
"Unknown cell zone name: " << cellZoneName_
233 <<
". Valid cell zones are: " <<
mesh.cellZones().names()
238 const label nCells = cellZoneCells.size();
239 const scalar nCellsTotal =
returnReduce(nCells, sumOp<label>());
241 const scalar VCellsTotal =
returnReduce(VCells, sumOp<scalar>());
242 Info<<
" cell zone size = " << nCellsTotal <<
endl;
243 Info<<
" cell zone volume = " << VCellsTotal <<
endl;
245 if ((nCellsTotal == 0) || (VCellsTotal*numberDensity_ < 1))
248 <<
"Number of particles to be added to cellZone " << cellZoneName_
249 <<
" is zero" <<
endl;
253 setPositions(cellZoneCells);
255 Info<<
" number density = " << numberDensity_ <<
nl 256 <<
" number of particles = " << positions_.size() <<
endl;
259 diameters_.setSize(positions_.size());
262 diameters_[i] = sizeDistribution_->sample();
271 template<
class CloudType>
279 template<
class CloudType>
286 if ((0.0 >= time0) && (0.0 < time1))
288 return positions_.size();
295 template<
class CloudType>
303 if ((0.0 >= time0) && (0.0 < time1))
305 return this->volumeTotal_;
312 template<
class CloudType>
316 const label nParcels,
324 position = positions_[parcelI];
325 cellOwner = injectorCells_[parcelI];
326 tetFacei = injectorTetFaces_[parcelI];
327 tetPti = injectorTetPts_[parcelI];
331 template<
class CloudType>
344 parcel.d() = diameters_[parcelI];
348 template<
class CloudType>
355 template<
class CloudType>
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
virtual ~CellZoneInjection()
Destructor.
errorManipArg< error, int > exit(error &err, const int errNo=1)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Templated injection model class.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
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 >> &tdf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
Lookup type of boundary radiation properties.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
#define forAll(list, i)
Loop across all elements in list.
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
Injection positions specified by a particle number density within a cell set.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< scalar > scalarList
A List of scalars.
constexpr scalar pi(M_PI)
ParcelType parcelType
Type of parcel the cloud was instantiated for.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
label parcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
T & last()
Access last element of the list, position [size()-1].
CellZoneInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
dimensionedScalar pow3(const dimensionedScalar &ds)
#define WarningInFunction
Report a warning using Foam::Warning.
virtual void setProperties(const label parcelI, const label nParcels, const scalar time, typename CloudType::parcelType &parcel)
Set the parcel properties.
Mesh data needed to do the Finite Volume discretisation.
A library of runtime-selectable doubly-truncated probability distribution models. Returns random samp...
messageStream Info
Information stream (stdout output on master, null elsewhere)
scalar timeEnd() const
Return the end-of-injection time.
virtual void setPositionAndCell(const label parcelI, const label nParcels, const scalar time, vector &position, label &cellOwner, label &tetFacei, label &tetPti)
Set the injection position and owner cell, tetFace and tetPt.
List< label > labelList
A List of labels.
Templated base class for dsmc cloud.
virtual void updateMesh()
Set injector locations when mesh is updated.
scalar volumeToInject(const scalar time0, const scalar time1)
Volume of parcels to introduce relative to SOI.
virtual bool validInjection(const label parcelI)
Return flag to identify whether or not injection of parcelI is.
static constexpr const zero Zero
Global zero (0)