patchInjectionBase.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2018 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 "patchInjectionBase.H"
27 #include "polyMesh.H"
28 #include "SubField.H"
29 #include "Random.H"
30 #include "triPointRef.H"
31 #include "volFields.H"
33 
34 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
35 
37 (
38  const polyMesh& mesh,
39  const word& patchName
40 )
41 :
42  patchName_(patchName),
43  patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
44  patchArea_(0.0),
45  patchNormal_(),
46  cellOwners_(),
47  triFace_(),
48  triToFace_(),
49  triCumulativeMagSf_(),
50  sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
51 {
52  if (patchId_ < 0)
53  {
55  << "Requested patch " << patchName_ << " not found" << nl
56  << "Available patches are: " << mesh.boundaryMesh().names() << nl
57  << exit(FatalError);
58  }
59 
60  updateMesh(mesh);
61 }
62 
63 
65 :
66  patchName_(pib.patchName_),
67  patchId_(pib.patchId_),
68  patchArea_(pib.patchArea_),
69  patchNormal_(pib.patchNormal_),
70  cellOwners_(pib.cellOwners_),
71  triFace_(pib.triFace_),
72  triToFace_(pib.triToFace_),
73  triCumulativeMagSf_(pib.triCumulativeMagSf_),
74  sumTriMagSf_(pib.sumTriMagSf_)
75 {}
76 
77 
78 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
79 
81 {}
82 
83 
84 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
85 
87 {
88  // Set/cache the injector cells
89  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
90  const pointField& points = patch.points();
91 
92  cellOwners_ = patch.faceCells();
93 
94  // Triangulate the patch faces and create addressing
95  DynamicList<label> triToFace(2*patch.size());
96  DynamicList<scalar> triMagSf(2*patch.size());
97  DynamicList<face> triFace(2*patch.size());
98  DynamicList<face> tris(5);
99 
100  // Set zero value at the start of the tri area list
101  triMagSf.append(0.0);
102 
103  forAll(patch, facei)
104  {
105  const face& f = patch[facei];
106 
107  tris.clear();
108  f.triangles(points, tris);
109 
110  forAll(tris, i)
111  {
112  triToFace.append(facei);
113  triFace.append(tris[i]);
114  triMagSf.append(tris[i].mag(points));
115  }
116  }
117 
118  forAll(sumTriMagSf_, i)
119  {
120  sumTriMagSf_[i] = 0.0;
121  }
122 
123  sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
124 
127 
128  for (label i = 1; i < triMagSf.size(); i++)
129  {
130  triMagSf[i] += triMagSf[i-1];
131  }
132 
133  // Transfer to persistent storage
135  triToFace_.transfer(triToFace);
136  triCumulativeMagSf_.transfer(triMagSf);
137 
138  // Convert sumTriMagSf_ into cumulative sum of areas per proc
139  for (label i = 1; i < sumTriMagSf_.size(); i++)
140  {
141  sumTriMagSf_[i] += sumTriMagSf_[i-1];
142  }
143 
144  const scalarField magSf(mag(patch.faceAreas()));
145  patchArea_ = sum(magSf);
146  patchNormal_ = patch.faceAreas()/magSf;
148 }
149 
150 
152 (
153  const fvMesh& mesh,
154  Random& rnd,
155  vector& position,
156  label& cellOwner,
157  label& tetFacei,
158  label& tetPti
159 )
160 {
161  scalar areaFraction = rnd.globalScalar01()*patchArea_;
162 
163  if (cellOwners_.size() > 0)
164  {
165  // Determine which processor to inject from
166  label proci = 0;
168  {
169  if (areaFraction >= sumTriMagSf_[i])
170  {
171  proci = i;
172  break;
173  }
174  }
175 
176  if (Pstream::myProcNo() == proci)
177  {
178  // Find corresponding decomposed face triangle
179  label trii = 0;
180  scalar offset = sumTriMagSf_[proci];
182  {
183  if (areaFraction > triCumulativeMagSf_[i] + offset)
184  {
185  trii = i;
186  break;
187  }
188  }
189 
190  // Set cellOwner
191  label facei = triToFace_[trii];
192  cellOwner = cellOwners_[facei];
193 
194  // Find random point in triangle
195  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
196  const pointField& points = patch.points();
197  const face& tf = triFace_[trii];
198  const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
199  const point pf(tri.randomPoint(rnd));
200 
201  // Position perturbed away from face (into domain)
202  const scalar a = rnd.scalarAB(0.1, 0.5);
203  const vector& pc = mesh.cellCentres()[cellOwner];
204  const vector d =
205  mag((pf - pc) & patchNormal_[facei])*patchNormal_[facei];
206 
207  position = pf - a*d;
208 
209  // Try to find tetFacei and tetPti in the current position
210  mesh.findTetFacePt(cellOwner, position, tetFacei, tetPti);
211 
212  // tetFacei and tetPti not found, check if the cell has changed
213  if (tetFacei == -1 ||tetPti == -1)
214  {
215  mesh.findCellFacePt(position, cellOwner, tetFacei, tetPti);
216  }
217 
218  // Both searches failed, choose a random position within
219  // the original cell
220  if (tetFacei == -1 ||tetPti == -1)
221  {
222  // Reset cellOwner
223  cellOwner = cellOwners_[facei];
224  const scalarField& V = mesh.V();
225 
226  // Construct cell tet indices
227  const List<tetIndices> cellTetIs =
229 
230  // Construct cell tet volume fractions
231  scalarList cTetVFrac(cellTetIs.size(), 0.0);
232  for (label teti=1; teti<cellTetIs.size()-1; teti++)
233  {
234  cTetVFrac[teti] =
235  cTetVFrac[teti-1]
236  + cellTetIs[teti].tet(mesh).mag()/V[cellOwner];
237  }
238  cTetVFrac.last() = 1;
239 
240  // Set new particle position
241  const scalar volFrac = rnd.sample01<scalar>();
242  label teti = 0;
243  forAll(cTetVFrac, vfI)
244  {
245  if (cTetVFrac[vfI] > volFrac)
246  {
247  teti = vfI;
248  break;
249  }
250  }
251  position = cellTetIs[teti].tet(mesh).randomPoint(rnd);
252  tetFacei = cellTetIs[teti].face();
253  tetPti = cellTetIs[teti].tetPt();
254  }
255  }
256  else
257  {
258  cellOwner = -1;
259  tetFacei = -1;
260  tetPti = -1;
261 
262  // Dummy position
263  position = pTraits<vector>::max;
264  }
265  }
266  else
267  {
268  cellOwner = -1;
269  tetFacei = -1;
270  tetPti = -1;
271 
272  // Dummy position
273  position = pTraits<vector>::max;
274  }
275 }
276 
277 
278 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
labelList triToFace_
Addressing from per triangle to patch face.
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
virtual void setPositionAndCell(const fvMesh &mesh, Random &rnd, vector &position, label &cellOwner, label &tetFacei, label &tetPti)
Set the injection position and owner cell, tetFace and tetPt.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
scalar scalarAB(const scalar a, const scalar b)
Advance the state and return a scalar sample from a uniform.
Definition: RandomI.H:63
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
Traits class for primitives.
Definition: pTraits.H:50
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:446
label findPatchID(const word &patchName) const
Find patch index given a name.
patchInjectionBase(const polyMesh &mesh, const word &patchName)
Construct from mesh and patch name.
virtual void updateMesh(const polyMesh &mesh)
Update patch geometry and derived info for injection locations.
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
const tensorField & tf
label triangles(const pointField &points, label &triI, faceList &triFaces) const
Split into triangles using existing points.
Definition: face.C:837
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
vectorList patchNormal_
Patch face normal directions.
Type sample01()
Advance the state and return a sample of a given type from a.
Definition: RandomI.H:70
const pointField & points
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:315
A class for handling words, derived from string.
Definition: word.H:59
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
void findCellFacePt(const point &p, label &celli, label &tetFacei, label &tetPti) const
Find the cell, tetFacei and tetPti for point p.
Definition: polyMesh.C:1455
scalarList sumTriMagSf_
Cumulative area fractions per processor.
face triFace(3)
wordList names() const
Return a list of patch names.
Base class for patch-based injection models.
const Field< PointType > & points() const
Return reference to global points.
const label patchId_
Patch ID.
scalar globalScalar01()
Advance the state and return a scalar sample from a uniform.
Definition: Random.C:94
scalarList triCumulativeMagSf_
Cumulative triangle area per triangle face.
const vectorField & cellCentres() const
Random number generator.
Definition: Random.H:57
const vectorField::subField faceAreas() const
Return face normals.
Definition: polyPatch.C:290
static const char nl
Definition: Ostream.H:265
labelList f(nPoints)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
scalar patchArea_
Patch area - total across all processors.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
labelList cellOwners_
List of cell labels corresponding to injector positions.
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:1480
Point randomPoint(Random &rndGen) const
Return a random point on the triangle from a uniform.
Definition: triangleI.H:241
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
T & last()
Return the last element of the list.
Definition: UListI.H:128
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
virtual ~patchInjectionBase()
Destructor.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
faceList triFace_
Decomposed patch faces as a list of triangles.