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-2021 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 #include "polygonTriangulate.H"
34 
35 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
36 
38 (
39  const polyMesh& mesh,
40  const word& patchName
41 )
42 :
43  patchName_(patchName),
44  patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
45  patchArea_(0.0),
46  patchNormal_(),
47  cellOwners_(),
48  triFace_(),
49  triToFace_(),
50  triCumulativeMagSf_(),
51  sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
52 {
53  if (patchId_ < 0)
54  {
56  << "Requested patch " << patchName_ << " not found" << nl
57  << "Available patches are: " << mesh.boundaryMesh().names() << nl
58  << exit(FatalError);
59  }
60 
61  updateMesh(mesh);
62 }
63 
64 
66 :
67  patchName_(pib.patchName_),
68  patchId_(pib.patchId_),
69  patchArea_(pib.patchArea_),
70  patchNormal_(pib.patchNormal_),
71  cellOwners_(pib.cellOwners_),
72  triFace_(pib.triFace_),
73  triToFace_(pib.triToFace_),
74  triCumulativeMagSf_(pib.triCumulativeMagSf_),
75  sumTriMagSf_(pib.sumTriMagSf_)
76 {}
77 
78 
79 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
80 
82 {}
83 
84 
85 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
86 
88 {
89  // Set/cache the injector cells
90  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
91  const pointField& points = patch.points();
92 
93  cellOwners_ = patch.faceCells();
94 
95  // Triangulate the patch faces and create addressing
96  DynamicList<label> triToFace(2*patch.size());
97  DynamicList<scalar> triMagSf(2*patch.size());
98  DynamicList<triFace> triFace(2*patch.size());
99 
100  // Set zero value at the start of the tri area list
101  triMagSf.append(0.0);
102 
103  // Create a triangulation engine
104  polygonTriangulate triEngine;
105 
106  forAll(patch, facei)
107  {
108  const face& f = patch[facei];
109 
110  triEngine.triangulate(UIndirectList<point>(points, f));
111 
112  forAll(triEngine.triPoints(), i)
113  {
114  triToFace.append(facei);
115  triFace.append(triEngine.triPoints(i, f));
116  triMagSf.append(triEngine.triPoints(i, f).mag(points));
117  }
118  }
119 
120  forAll(sumTriMagSf_, i)
121  {
122  sumTriMagSf_[i] = 0.0;
123  }
124 
125  sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
126 
129 
130  for (label i = 1; i < triMagSf.size(); i++)
131  {
132  triMagSf[i] += triMagSf[i-1];
133  }
134 
135  // Transfer to persistent storage
137  triToFace_.transfer(triToFace);
138  triCumulativeMagSf_.transfer(triMagSf);
139 
140  // Convert sumTriMagSf_ into cumulative sum of areas per proc
141  for (label i = 1; i < sumTriMagSf_.size(); i++)
142  {
143  sumTriMagSf_[i] += sumTriMagSf_[i-1];
144  }
145 
146  const scalarField magSf(patch.magFaceAreas());
147  patchArea_ = sum(magSf);
148  patchNormal_ = patch.faceAreas()/magSf;
149  reduce(patchArea_, sumOp<scalar>());
150 }
151 
152 
154 (
155  const fvMesh& mesh,
156  Random& rnd,
157  vector& position,
158  label& cellOwner,
159  label& tetFacei,
160  label& tetPti
161 )
162 {
163  scalar areaFraction = rnd.globalScalar01()*patchArea_;
164 
165  if (cellOwners_.size() > 0)
166  {
167  // Determine which processor to inject from
168  label proci = 0;
170  {
171  if (areaFraction >= sumTriMagSf_[i])
172  {
173  proci = i;
174  break;
175  }
176  }
177 
178  if (Pstream::myProcNo() == proci)
179  {
180  // Find corresponding decomposed face triangle
181  label trii = 0;
182  scalar offset = sumTriMagSf_[proci];
184  {
185  if (areaFraction > triCumulativeMagSf_[i] + offset)
186  {
187  trii = i;
188  break;
189  }
190  }
191 
192  // Set cellOwner
193  label facei = triToFace_[trii];
194  cellOwner = cellOwners_[facei];
195 
196  // Find random point in triangle
197  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
198  const pointField& points = patch.points();
199  const face& tf = triFace_[trii];
200  const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
201  const point pf(tri.randomPoint(rnd));
202 
203  // Position perturbed away from face (into domain)
204  const scalar a = rnd.scalarAB(0.1, 0.5);
205  const vector& pc = mesh.cellCentres()[cellOwner];
206  const vector d =
207  mag((pf - pc) & patchNormal_[facei])*patchNormal_[facei];
208 
209  position = pf - a*d;
210 
211  // Try to find tetFacei and tetPti in the current position
212  mesh.findTetFacePt(cellOwner, position, tetFacei, tetPti);
213 
214  // tetFacei and tetPti not found, check if the cell has changed
215  if (tetFacei == -1 ||tetPti == -1)
216  {
217  mesh.findCellFacePt(position, cellOwner, tetFacei, tetPti);
218  }
219 
220  // Both searches failed, choose a random position within
221  // the original cell
222  if (tetFacei == -1 ||tetPti == -1)
223  {
224  // Reset cellOwner
225  cellOwner = cellOwners_[facei];
226  const scalarField& V = mesh.V();
227 
228  // Construct cell tet indices
229  const List<tetIndices> cellTetIs =
231 
232  // Construct cell tet volume fractions
233  scalarList cTetVFrac(cellTetIs.size(), 0.0);
234  for (label teti=1; teti<cellTetIs.size()-1; teti++)
235  {
236  cTetVFrac[teti] =
237  cTetVFrac[teti-1]
238  + cellTetIs[teti].tet(mesh).mag()/V[cellOwner];
239  }
240  cTetVFrac.last() = 1;
241 
242  // Set new particle position
243  const scalar volFrac = rnd.sample01<scalar>();
244  label teti = 0;
245  forAll(cTetVFrac, vfI)
246  {
247  if (cTetVFrac[vfI] > volFrac)
248  {
249  teti = vfI;
250  break;
251  }
252  }
253  position = cellTetIs[teti].tet(mesh).randomPoint(rnd);
254  tetFacei = cellTetIs[teti].face();
255  tetPti = cellTetIs[teti].tetPt();
256  }
257  }
258  else
259  {
260  cellOwner = -1;
261  tetFacei = -1;
262  tetPti = -1;
263 
264  // Dummy position
265  position = pTraits<vector>::max;
266  }
267  }
268  else
269  {
270  cellOwner = -1;
271  tetFacei = -1;
272  tetPti = -1;
273 
274  // Dummy position
275  position = pTraits<vector>::max;
276  }
277 }
278 
279 
280 // ************************************************************************* //
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:323
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:164
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
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:334
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
triFaceList triFace_
Decomposed patch faces as a list of triangles.
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 areas.
Definition: polyPatch.C:303
static const char nl
Definition: Ostream.H:260
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.
const scalarField::subField magFaceAreas() const
Return face area magnitudes.
Definition: polyPatch.C:309
Triangulation of three-dimensional polygons.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A List with indirect addressing.
Definition: fvMatrix.H:106
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
const UList< triFace > & triPoints() const
Get the triangles&#39; points.
virtual ~patchInjectionBase()
Destructor.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342