patchInjectionBase.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2013-2017 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "patchInjectionBase.H"
29#include "polyMesh.H"
30#include "SubField.H"
31#include "Random.H"
32#include "triPointRef.H"
33#include "volFields.H"
35
36// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
37
39(
40 const polyMesh& mesh,
41 const word& patchName
42)
43:
44 patchName_(patchName),
45 patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
46 patchArea_(0.0),
47 patchNormal_(),
48 cellOwners_(),
49 triFace_(),
50 triToFace_(),
51 triCumulativeMagSf_(),
52 sumTriMagSf_(Pstream::nProcs() + 1, Zero)
53{
54 if (patchId_ < 0)
55 {
57 << "Requested patch " << patchName_ << " not found" << nl
58 << "Available patches are: " << mesh.boundaryMesh().names() << nl
59 << exit(FatalError);
60 }
61
63}
64
65
67:
68 patchName_(pib.patchName_),
69 patchId_(pib.patchId_),
70 patchArea_(pib.patchArea_),
71 patchNormal_(pib.patchNormal_),
72 cellOwners_(pib.cellOwners_),
73 triFace_(pib.triFace_),
74 triToFace_(pib.triToFace_),
75 triCumulativeMagSf_(pib.triCumulativeMagSf_),
76 sumTriMagSf_(pib.sumTriMagSf_)
77{}
78
79
80// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
81
83{
84 // Set/cache the injector cells
85 const polyPatch& patch = mesh.boundaryMesh()[patchId_];
86 const pointField& points = patch.points();
87
88 cellOwners_ = patch.faceCells();
89
90 // Triangulate the patch faces and create addressing
91 DynamicList<label> triToFace(2*patch.size());
92 DynamicList<scalar> triMagSf(2*patch.size());
93 DynamicList<face> triFace(2*patch.size());
94 DynamicList<face> tris(5);
95
96 // Set zero value at the start of the tri area list
97 triMagSf.append(0.0);
98
99 forAll(patch, facei)
100 {
101 const face& f = patch[facei];
102
103 tris.clear();
104 f.triangles(points, tris);
105
106 forAll(tris, i)
107 {
108 triToFace.append(facei);
109 triFace.append(tris[i]);
110 triMagSf.append(tris[i].mag(points));
111 }
112 }
113
114 sumTriMagSf_ = Zero;
115 sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
118
119 for (label i = 1; i < triMagSf.size(); i++)
120 {
121 triMagSf[i] += triMagSf[i-1];
122 }
123
124 // Transfer to persistent storage
125 triFace_.transfer(triFace);
126 triToFace_.transfer(triToFace);
127 triCumulativeMagSf_.transfer(triMagSf);
128
129 // Convert sumTriMagSf_ into cumulative sum of areas per proc
130 for (label i = 1; i < sumTriMagSf_.size(); i++)
131 {
132 sumTriMagSf_[i] += sumTriMagSf_[i-1];
133 }
135 const scalarField magSf(mag(patch.faceAreas()));
136 patchArea_ = sum(magSf);
137 patchNormal_ = patch.faceAreas()/magSf;
138 reduce(patchArea_, sumOp<scalar>());
139}
140
141
143(
144 const fvMesh& mesh,
145 const scalar fraction01,
146 Random& rnd,
147 vector& position,
148 label& cellOwner,
149 label& tetFacei,
150 label& tetPti
151)
152{
153 label facei = -1;
154
155 if (cellOwners_.size() > 0)
156 {
157 // Determine which processor to inject from
158 const label proci = whichProc(fraction01);
159
160 if (Pstream::myProcNo() == proci)
161 {
162 const scalar areaFraction = fraction01*patchArea_;
163
164 // Find corresponding decomposed face triangle
165 label trii = 0;
166 scalar offset = sumTriMagSf_[proci];
167 forAllReverse(triCumulativeMagSf_, i)
168 {
169 if (areaFraction > triCumulativeMagSf_[i] + offset)
170 {
171 trii = i;
172 break;
173 }
174 }
175
176 // Set cellOwner
177 facei = triToFace_[trii];
178 cellOwner = cellOwners_[facei];
179
180 // Find random point in triangle
181 const polyPatch& patch = mesh.boundaryMesh()[patchId_];
182 const pointField& points = patch.points();
183 const face& tf = triFace_[trii];
184 const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
185 const point pf(tri.randomPoint(rnd));
186
187 // Position perturbed away from face (into domain)
188 const scalar a = rnd.position(scalar(0.1), scalar(0.5));
189 const vector& pc = mesh.cellCentres()[cellOwner];
190 const vector d =
191 mag((pf - pc) & patchNormal_[facei])*patchNormal_[facei];
192
193 position = pf - a*d;
194
195 // Try to find tetFacei and tetPti in the current position
196 mesh.findTetFacePt(cellOwner, position, tetFacei, tetPti);
197
198 // tetFacei and tetPti not found, check if the cell has changed
199 if (tetFacei == -1 ||tetPti == -1)
200 {
201 mesh.findCellFacePt(position, cellOwner, tetFacei, tetPti);
202 }
203
204 // Both searches failed, choose a random position within
205 // the original cell
206 if (tetFacei == -1 ||tetPti == -1)
207 {
208 // Reset cellOwner
209 cellOwner = cellOwners_[facei];
210 const scalarField& V = mesh.V();
211
212 // Construct cell tet indices
213 const List<tetIndices> cellTetIs =
215
216 // Construct cell tet volume fractions
217 scalarList cTetVFrac(cellTetIs.size(), Zero);
218 for (label teti=1; teti<cellTetIs.size()-1; teti++)
219 {
220 cTetVFrac[teti] =
221 cTetVFrac[teti-1]
222 + cellTetIs[teti].tet(mesh).mag()/V[cellOwner];
223 }
224 cTetVFrac.last() = 1;
225
226 // Set new particle position
227 const scalar volFrac = rnd.sample01<scalar>();
228 label teti = 0;
229 forAll(cTetVFrac, vfI)
230 {
231 if (cTetVFrac[vfI] > volFrac)
232 {
233 teti = vfI;
234 break;
235 }
236 }
237 position = cellTetIs[teti].tet(mesh).randomPoint(rnd);
238 tetFacei = cellTetIs[teti].face();
239 tetPti = cellTetIs[teti].tetPt();
240 }
241 }
242 else
243 {
244 cellOwner = -1;
245 tetFacei = -1;
246 tetPti = -1;
247
248 // Dummy position
249 position = pTraits<vector>::max;
250 }
251 }
252 else
253 {
254 cellOwner = -1;
255 tetFacei = -1;
256 tetPti = -1;
257
258 // Dummy position
259 position = pTraits<vector>::max;
260 }
261
262 return facei;
263}
264
265
267(
268 const fvMesh& mesh,
269 Random& rnd,
270 vector& position,
271 label& cellOwner,
272 label& tetFacei,
273 label& tetPti
274)
275{
276 scalar fraction01 = rnd.globalSample01<scalar>();
277
278 return setPositionAndCell
279 (
280 mesh,
281 fraction01,
282 rnd,
283 position,
284 cellOwner,
285 tetFacei,
286 tetPti
287 );
288}
289
290
291Foam::label Foam::patchInjectionBase::whichProc(const scalar fraction01) const
292{
293 const scalar areaFraction = fraction01*patchArea_;
294
295 // Determine which processor to inject from
296 forAllReverse(sumTriMagSf_, i)
297 {
298 if (areaFraction >= sumTriMagSf_[i])
299 {
300 return i;
301 }
302 }
303
304 return 0;
305}
306
307
308// ************************************************************************* //
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:391
void transfer(List< T > &list)
Transfer contents of the argument List into this.
Definition: DynamicListI.H:467
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
Inter-processor communications stream.
Definition: Pstream.H:63
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
Random number generator.
Definition: Random.H:60
Type position(const Type &start, const Type &end)
Return a sample on the interval [start,end].
Type sample01()
Return a sample whose components lie in the range [0,1].
Type globalSample01()
Return a sample whose components lie in the range [0,1].
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
void updateMesh()
Update for new mesh topology.
A traits class, which is primarily used for primitives.
Definition: pTraits.H:59
const word patchName_
Patch name.
const label patchId_
Patch ID.
virtual void updateMesh(const polyMesh &mesh)
Update patch geometry and derived info for injection locations.
label whichProc(const scalar fraction01) const
Return the processor that has the location specified by the fraction.
label setPositionAndCell(const fvMesh &mesh, const scalar fraction01, Random &rnd, vector &position, label &cellOwner, label &tetFacei, label &tetPti)
Set the injection position and owner cell, tetFace and tetPt.
wordList names() const
Return a list of patch names.
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
void findCellFacePt(const point &p, label &celli, label &tetFacei, label &tetPti) const
Find the cell, tetFacei and tetPti for point p.
Definition: polyMesh.C:1371
void findTetFacePt(const label celli, const point &p, label &tetFacei, label &tetPti) const
Find the tetFacei and tetPti for point p in celli.
Definition: polyMesh.C:1396
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
const vectorField & cellCentres() const
int myProcNo() const noexcept
Return processor number.
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:72
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
Point randomPoint(Random &rndGen) const
Return a random point on the triangle from a uniform.
Definition: triangleI.H:254
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:346