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 | Copyright (C) 2013-2015 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 "cachedRandom.H"
30 #include "triPointRef.H"
31 
32 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
33 
35 (
36  const polyMesh& mesh,
37  const word& patchName
38 )
39 :
40  patchName_(patchName),
41  patchId_(mesh.boundaryMesh().findPatchID(patchName_)),
42  patchArea_(0.0),
43  patchNormal_(),
44  cellOwners_(),
45  triFace_(),
46  triToFace_(),
47  triCumulativeMagSf_(),
48  sumTriMagSf_(Pstream::nProcs() + 1, 0.0)
49 {
50  if (patchId_ < 0)
51  {
53  (
54  "patchInjectionBase::patchInjectionBase"
55  "("
56  "const polyMesh& mesh, "
57  "const word& patchName"
58  ")"
59  ) << "Requested patch " << patchName_ << " not found" << nl
60  << "Available patches are: " << mesh.boundaryMesh().names() << nl
61  << exit(FatalError);
62  }
63 
64  updateMesh(mesh);
65 }
66 
67 
69 :
70  patchName_(pib.patchName_),
71  patchId_(pib.patchId_),
72  patchArea_(pib.patchArea_),
73  patchNormal_(pib.patchNormal_),
74  cellOwners_(pib.cellOwners_),
75  triFace_(pib.triFace_),
76  triToFace_(pib.triToFace_),
77  triCumulativeMagSf_(pib.triCumulativeMagSf_),
78  sumTriMagSf_(pib.sumTriMagSf_)
79 {}
80 
81 
82 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
83 
85 {}
86 
87 
88 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
89 
91 {
92  // Set/cache the injector cells
93  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
94  const pointField& points = patch.points();
95 
96  cellOwners_ = patch.faceCells();
97 
98  // Triangulate the patch faces and create addressing
99  DynamicList<label> triToFace(2*patch.size());
100  DynamicList<scalar> triMagSf(2*patch.size());
101  DynamicList<face> triFace(2*patch.size());
102  DynamicList<face> tris(5);
103 
104  // Set zero value at the start of the tri area list
105  triMagSf.append(0.0);
106 
107  forAll(patch, faceI)
108  {
109  const face& f = patch[faceI];
110 
111  tris.clear();
112  f.triangles(points, tris);
113 
114  forAll(tris, i)
115  {
116  triToFace.append(faceI);
117  triFace.append(tris[i]);
118  triMagSf.append(tris[i].mag(points));
119  }
120  }
121 
122  forAll(sumTriMagSf_, i)
123  {
124  sumTriMagSf_[i] = 0.0;
125  }
126 
127  sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
128 
131 
132  for (label i = 1; i < triMagSf.size(); i++)
133  {
134  triMagSf[i] += triMagSf[i-1];
135  }
136 
137  // Transfer to persistent storage
139  triToFace_.transfer(triToFace);
140  triCumulativeMagSf_.transfer(triMagSf);
141 
142  // Convert sumTriMagSf_ into cumulative sum of areas per proc
143  for (label i = 1; i < sumTriMagSf_.size(); i++)
144  {
145  sumTriMagSf_[i] += sumTriMagSf_[i-1];
146  }
147 
148  const scalarField magSf(mag(patch.faceAreas()));
149  patchArea_ = sum(magSf);
150  patchNormal_ = patch.faceAreas()/magSf;
152 }
153 
154 
156 (
157  const polyMesh& mesh,
158  cachedRandom& rnd,
159  vector& position,
160  label& cellOwner,
161  label& tetFaceI,
162  label& tetPtI
163 )
164 {
165  scalar areaFraction = 0;
166 
167  if (Pstream::master())
168  {
169  areaFraction = rnd.position<scalar>(0, patchArea_);
170  }
171 
172  Pstream::scatter(areaFraction);
173 
174  if (cellOwners_.size() > 0)
175  {
176  // Determine which processor to inject from
177  label procI = 0;
179  {
180  if (areaFraction >= sumTriMagSf_[i])
181  {
182  procI = i;
183  break;
184  }
185  }
186 
187  if (Pstream::myProcNo() == procI)
188  {
189  // Find corresponding decomposed face triangle
190  label triI = 0;
191  scalar offset = sumTriMagSf_[procI];
193  {
194  if (areaFraction > triCumulativeMagSf_[i] + offset)
195  {
196  triI = i;
197  break;
198  }
199  }
200 
201  // Set cellOwner
202  label faceI = triToFace_[triI];
203  cellOwner = cellOwners_[faceI];
204 
205  // Find random point in triangle
206  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
207  const pointField& points = patch.points();
208  const face& tf = triFace_[triI];
209  const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
210  const point pf(tri.randomPoint(rnd));
211 
212  // Position perturbed away from face (into domain)
213  const scalar a = rnd.position(scalar(0.1), scalar(0.5));
214  const vector& pc = mesh.cellCentres()[cellOwner];
215  const vector d = mag(pf - pc)*patchNormal_[faceI];
216 
217  position = pf - a*d;
218 
219  // The position is between the face and cell centre, which could
220  // be in any tet of the decomposed cell, so arbitrarily choose the
221  // first face of the cell as the tetFace and the first point after
222  // the base point on the face as the tetPt. The tracking will pick
223  // the cell consistent with the motion in the first tracking step
224  tetFaceI = mesh.cells()[cellOwner][0];
225  tetPtI = 1;
226  }
227  else
228  {
229  cellOwner = -1;
230  tetFaceI = -1;
231  tetPtI = -1;
232 
233  // Dummy position
234  position = pTraits<vector>::max;
235  }
236  }
237  else
238  {
239  cellOwner = -1;
240  tetFaceI = -1;
241  tetPtI = -1;
242 
243  // Dummy position
244  position = pTraits<vector>::max;
245  }
246 }
247 
248 
249 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Point randomPoint(Random &rndGen) const
Return a random point on the triangle from a uniform.
Definition: triangleI.H:232
const pointField & points
scalar patchArea_
Patch area - total across all processors.
labelList triToFace_
Addressing from per triangle to patch face.
const vectorField::subField faceAreas() const
Return face normals.
Definition: polyPatch.C:290
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList f(nPoints)
patchInjectionBase(const polyMesh &mesh, const word &patchName)
Construct from mesh and patch name.
const label patchId_
Patch ID.
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:316
faceList triFace_
Decomposed patch faces as a list of triangles.
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:242
A class for handling words, derived from string.
Definition: word.H:59
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const cellList & cells() const
face triFace(3)
const vectorField & cellCentres() const
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:57
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual ~patchInjectionBase()
Destructor.
wordList names() const
Return a list of patch names.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char nl
Definition: Ostream.H:260
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
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)
Definition: UList.H:421
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
virtual void setPositionAndCell(const polyMesh &mesh, cachedRandom &rnd, vector &position, label &cellOwner, label &tetFaceI, label &tetPtI)
Set the injection position and owner cell, tetFace and tetPt.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
scalarList triCumulativeMagSf_
Cumulative triangle area per triangle face.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:404
const tensorField & tf
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
#define forAllReverse(list, i)
Definition: UList.H:424
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Traits class for primitives.
Definition: pTraits.H:50
error FatalError
Random number generator.
Definition: cachedRandom.H:63
vectorList patchNormal_
Patch face normal directions.
virtual void updateMesh(const polyMesh &mesh)
Update patch geometry and derived info for injection locations.
scalarList sumTriMagSf_
Cumulative area fractions per processor.
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:68
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:398
Type position(const Type &start, const Type &end)
Return a sample between start and end.
label findPatchID(const word &patchName) const
Find patch index given a name.
label triangles(const pointField &points, label &triI, faceList &triFaces) const
Split into triangles using existing points.
Definition: face.C:835
const Field< PointType > & points() const
Return reference to global points.
labelList cellOwners_
List of cell labels corresponding to injector positions.