CellZoneInjection.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "CellZoneInjection.H"
27 #include "mathematicalConstants.H"
29 #include "globalIndex.H"
30 #include "Pstream.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 template<class CloudType>
36 (
37  const labelList& cellZoneCells
38 )
39 {
40  const fvMesh& mesh = this->owner().mesh();
41  const scalarField& V = mesh.V();
42  const label nCells = cellZoneCells.size();
43  cachedRandom& rnd = this->owner().rndGen();
44 
45  DynamicList<vector> positions(nCells); // initial size only
46  DynamicList<label> injectorCells(nCells); // initial size only
47  DynamicList<label> injectorTetFaces(nCells); // initial size only
48  DynamicList<label> injectorTetPts(nCells); // initial size only
49 
50  scalar newParticlesTotal = 0.0;
51  label addParticlesTotal = 0;
52 
53  forAll(cellZoneCells, i)
54  {
55  const label celli = cellZoneCells[i];
56 
57  // Calc number of particles to add
58  const scalar newParticles = V[celli]*numberDensity_;
59  newParticlesTotal += newParticles;
60  label addParticles = floor(newParticles);
61  addParticlesTotal += addParticles;
62 
63  const scalar diff = newParticlesTotal - addParticlesTotal;
64  if (diff > 1)
65  {
66  label corr = floor(diff);
67  addParticles += corr;
68  addParticlesTotal += corr;
69  }
70 
71  // Construct cell tet indices
72  const List<tetIndices> cellTetIs =
73  polyMeshTetDecomposition::cellTetIndices(mesh, celli);
74 
75  // Construct cell tet volume fractions
76  scalarList cTetVFrac(cellTetIs.size(), 0.0);
77  for (label tetI = 1; tetI < cellTetIs.size() - 1; tetI++)
78  {
79  cTetVFrac[tetI] =
80  cTetVFrac[tetI-1] + cellTetIs[tetI].tet(mesh).mag()/V[celli];
81  }
82  cTetVFrac.last() = 1.0;
83 
84  // Set new particle position and cellId
85  for (label pI = 0; pI < addParticles; pI++)
86  {
87  const scalar volFrac = rnd.sample01<scalar>();
88  label tetI = 0;
89  forAll(cTetVFrac, vfI)
90  {
91  if (cTetVFrac[vfI] > volFrac)
92  {
93  tetI = vfI;
94  break;
95  }
96  }
97  positions.append(cellTetIs[tetI].tet(mesh).randomPoint(rnd));
98 
99  injectorCells.append(celli);
100  injectorTetFaces.append(cellTetIs[tetI].face());
101  injectorTetPts.append(cellTetIs[tetI].tetPt());
102  }
103  }
104 
105  // Parallel operation manipulations
106  globalIndex globalPositions(positions.size());
107  List<vector> allPositions(globalPositions.size(), point::max);
108  List<label> allInjectorCells(globalPositions.size(), -1);
109  List<label> allInjectorTetFaces(globalPositions.size(), -1);
110  List<label> allInjectorTetPts(globalPositions.size(), -1);
111 
112  // Gather all positions on to all processors
113  SubList<vector>
114  (
115  allPositions,
116  globalPositions.localSize(Pstream::myProcNo()),
117  globalPositions.offset(Pstream::myProcNo())
118  ) = positions;
119 
120  Pstream::listCombineGather(allPositions, minEqOp<point>());
121  Pstream::listCombineScatter(allPositions);
122 
123  // Gather local cell tet and tet-point Ids, but leave non-local ids set -1
124  SubList<label>
125  (
126  allInjectorCells,
127  globalPositions.localSize(Pstream::myProcNo()),
128  globalPositions.offset(Pstream::myProcNo())
129  ) = injectorCells;
130  SubList<label>
131  (
132  allInjectorTetFaces,
133  globalPositions.localSize(Pstream::myProcNo()),
134  globalPositions.offset(Pstream::myProcNo())
135  ) = injectorTetFaces;
136  SubList<label>
137  (
138  allInjectorTetPts,
139  globalPositions.localSize(Pstream::myProcNo()),
140  globalPositions.offset(Pstream::myProcNo())
141  ) = injectorTetPts;
142 
143  // Transfer data
144  positions_.transfer(allPositions);
145  injectorCells_.transfer(allInjectorCells);
146  injectorTetFaces_.transfer(allInjectorTetFaces);
147  injectorTetPts_.transfer(allInjectorTetPts);
148 
149 
150  if (debug)
151  {
152  OFstream points("points.obj");
153  forAll(positions_, i)
154  {
155  meshTools::writeOBJ(points, positions_[i]);
156  }
157  }
158 }
159 
160 
161 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
162 
163 template<class CloudType>
165 (
166  const dictionary& dict,
167  CloudType& owner,
168  const word& modelName
169 )
170 :
171  InjectionModel<CloudType>(dict, owner, modelName, typeName),
172  cellZoneName_(this->coeffDict().lookup("cellZone")),
173  numberDensity_(readScalar(this->coeffDict().lookup("numberDensity"))),
174  positions_(),
175  injectorCells_(),
176  injectorTetFaces_(),
177  injectorTetPts_(),
178  diameters_(),
179  U0_(this->coeffDict().lookup("U0")),
180  sizeDistribution_
181  (
183  (
184  this->coeffDict().subDict("sizeDistribution"), owner.rndGen()
185  )
186  )
187 {
188  updateMesh();
189 }
190 
191 
192 template<class CloudType>
194 (
196 )
197 :
199  cellZoneName_(im.cellZoneName_),
200  numberDensity_(im.numberDensity_),
201  positions_(im.positions_),
202  injectorCells_(im.injectorCells_),
203  injectorTetFaces_(im.injectorTetFaces_),
204  injectorTetPts_(im.injectorTetPts_),
205  diameters_(im.diameters_),
206  U0_(im.U0_),
207  sizeDistribution_(im.sizeDistribution_().clone().ptr())
208 {}
209 
210 
211 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
212 
213 template<class CloudType>
215 {}
216 
217 
218 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
219 
220 template<class CloudType>
222 {
223  // Set/cache the injector cells
224  const fvMesh& mesh = this->owner().mesh();
225  const label zoneI = mesh.cellZones().findZoneID(cellZoneName_);
226 
227  if (zoneI < 0)
228  {
230  << "Unknown cell zone name: " << cellZoneName_
231  << ". Valid cell zones are: " << mesh.cellZones().names()
232  << nl << exit(FatalError);
233  }
234 
235  const labelList& cellZoneCells = mesh.cellZones()[zoneI];
236  const label nCells = cellZoneCells.size();
237  const scalar nCellsTotal = returnReduce(nCells, sumOp<label>());
238  const scalar VCells = sum(scalarField(mesh.V(), cellZoneCells));
239  const scalar VCellsTotal = returnReduce(VCells, sumOp<scalar>());
240  Info<< " cell zone size = " << nCellsTotal << endl;
241  Info<< " cell zone volume = " << VCellsTotal << endl;
242 
243  if ((nCellsTotal == 0) || (VCellsTotal*numberDensity_ < 1))
244  {
246  << "Number of particles to be added to cellZone " << cellZoneName_
247  << " is zero" << endl;
248  }
249  else
250  {
251  setPositions(cellZoneCells);
252 
253  Info<< " number density = " << numberDensity_ << nl
254  << " number of particles = " << positions_.size() << endl;
255 
256  // Construct parcel diameters
257  diameters_.setSize(positions_.size());
258  forAll(diameters_, i)
259  {
260  diameters_[i] = sizeDistribution_->sample();
261  }
262  }
263 
264  // Determine volume of particles to inject
265  this->volumeTotal_ = sum(pow3(diameters_))*constant::mathematical::pi/6.0;
266 }
267 
268 
269 template<class CloudType>
271 {
272  // Not used
273  return this->SOI_;
274 }
275 
276 
277 template<class CloudType>
279 (
280  const scalar time0,
281  const scalar time1
282 )
283 {
284  if ((0.0 >= time0) && (0.0 < time1))
285  {
286  return positions_.size();
287  }
288  else
289  {
290  return 0;
291  }
292 }
293 
294 
295 template<class CloudType>
297 (
298  const scalar time0,
299  const scalar time1
300 )
301 {
302  // All parcels introduced at SOI
303  if ((0.0 >= time0) && (0.0 < time1))
304  {
305  return this->volumeTotal_;
306  }
307  else
308  {
309  return 0.0;
310  }
311 }
312 
313 
314 template<class CloudType>
316 (
317  const label parcelI,
318  const label nParcels,
319  const scalar time,
320  vector& position,
321  label& cellOwner,
322  label& tetFacei,
323  label& tetPtI
324 )
325 {
326  position = positions_[parcelI];
327  cellOwner = injectorCells_[parcelI];
328  tetFacei = injectorTetFaces_[parcelI];
329  tetPtI = injectorTetPts_[parcelI];
330 }
331 
332 
333 template<class CloudType>
335 (
336  const label parcelI,
337  const label,
338  const scalar,
339  typename CloudType::parcelType& parcel
340 )
341 {
342  // set particle velocity
343  parcel.U() = U0_;
344 
345  // set particle diameter
346  parcel.d() = diameters_[parcelI];
347 }
348 
349 
350 template<class CloudType>
352 {
353  return false;
354 }
355 
356 
357 template<class CloudType>
359 {
360  return true;
361 }
362 
363 
364 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
virtual ~CellZoneInjection()
Destructor.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Templated injection model class.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: ZoneMesh.C:341
stressControl lookup("compactNormalStress") >> compactNormalStress
const pointField & points
Injection positions specified by a particle number density within a cell set.
A class for handling words, derived from string.
Definition: word.H:59
wordList names() const
Return a list of zone names.
Definition: ZoneMesh.C:256
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
List< label > labelList
A List of labels.
Definition: labelList.H:56
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:217
volScalarField scalarField(fieldObject, mesh)
static const char nl
Definition: Ostream.H:262
label parcelsToInject(const scalar time0, const scalar time1)
Number of parcels to introduce relative to SOI.
CellZoneInjection(const dictionary &dict, CloudType &owner, const word &modelName)
Construct from dictionary.
dimensionedScalar pow3(const dimensionedScalar &ds)
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:469
#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.
Definition: fvMesh.H:78
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
virtual bool fullyDescribed() const
Flag to identify whether model fully describes the parcel.
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.
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
scalar timeEnd() const
Return the end-of-injection time.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
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.