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-2023 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  sumProcArea_(),
46  sumFaceArea_(),
47  sumFaceTriArea_()
48 {
49  if (patchId_ < 0)
50  {
52  << "Requested patch " << patchName_ << " not found" << nl
53  << "Available patches are: " << mesh.boundaryMesh().names() << nl
54  << exit(FatalError);
55  }
56 
57  topoChange(mesh);
58 }
59 
60 
62 :
63  patchName_(pib.patchName_),
64  patchId_(pib.patchId_),
65  sumProcArea_(pib.sumProcArea_),
66  sumFaceArea_(pib.sumFaceArea_),
67  sumFaceTriArea_(pib.sumFaceTriArea_)
68 {}
69 
70 
71 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
72 
74 {}
75 
76 
77 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
78 
80 {
81  // Set/cache the injector cells
82  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
83 
84  // Initialise
85  sumProcArea_.resize(Pstream::nProcs());
86  sumProcArea_ = 0;
87  sumFaceArea_.resize(patch.size());
88  sumFaceArea_ = 0;
89  sumFaceTriArea_.resize(patch.size());
90  forAll(patch, patchFacei)
91  {
92  sumFaceTriArea_[patchFacei].resize(patch[patchFacei].nTriangles());
93  sumFaceTriArea_[patchFacei] = 0;
94  }
95 
96  // Cumulatively sum up the areas through the local patch
97  scalar patchArea = 0;
98  forAll(patch, patchFacei)
99  {
100  const label facei = patchFacei + patch.start();
101  const label celli = patch.faceCells()[patchFacei];
102 
103  scalar patchFaceArea = 0;
104  for
105  (
106  label patchFaceTrii = 0;
107  patchFaceTrii < patch[patchFacei].nTriangles();
108  ++ patchFaceTrii
109  )
110  {
111  const tetIndices tet(celli, facei, patchFaceTrii + 1);
112 
113  patchFaceArea += tet.faceTri(mesh).mag();
114  sumFaceTriArea_[patchFacei][patchFaceTrii] = patchFaceArea;
115  }
116 
117  patchArea += patchFaceArea;
118  sumFaceArea_[patchFacei] = patchArea;
119  }
120 
121  // Cumulatively sum the total areas across the processors
122  sumProcArea_[Pstream::myProcNo()] = patchArea;
123  Pstream::listCombineGather(sumProcArea_, maxEqOp<scalar>());
124  Pstream::listCombineScatter(sumProcArea_);
125  for (label proci = 1; proci < Pstream::nProcs(); proci ++)
126  {
127  sumProcArea_[proci] += sumProcArea_[proci - 1];
128  }
129 }
130 
131 
133 (
134  const fvMesh& mesh,
135  Random& rnd,
136  barycentric& coordinates,
137  label& celli,
138  label& tetFacei,
139  label& tetPti,
140  label& facei
141 )
142 {
143  // Linear searching for area fractions
144  auto findArea = []
145  (
146  const scalarList& areas,
147  scalar& area
148  )
149  {
150  for (label i = areas.size(); i > 0; i --)
151  {
152  if (area >= areas[i - 1])
153  {
154  area -= areas[i - 1];
155  return i;
156  }
157  }
158 
159  return label(0);
160  };
161 
162  const polyPatch& patch = mesh.boundaryMesh()[patchId_];
163 
164  scalar area = rnd.globalScalar01()*sumProcArea_.last();
165 
166  if (patch.size() > 0)
167  {
168  // Determine which processor to inject from
169  const label proci = findArea(sumProcArea_, area);
170 
171  // If this processor...
172  if (Pstream::myProcNo() == proci)
173  {
174  // Determine the face to inject from
175  const label patchFacei = findArea(sumFaceArea_, area);
176 
177  // Determine the face-tri to inject from
178  const label patchFaceTrii =
179  findArea(sumFaceTriArea_[patchFacei], area);
180 
181  // Set the topology
182  const barycentric2D r = barycentric2D01(rnd);
183  coordinates = barycentric(0, r.a(), r.b(), r.c());
184  celli = mesh.faceOwner()[patch.start() + patchFacei];
185  tetFacei = patch.start() + patchFacei;
186  tetPti = patchFaceTrii + 1;
187  facei = patch.start() + patchFacei;
188  }
189  else
190  {
191  coordinates = barycentric::uniform(NaN);
192  celli = -1;
193  tetFacei = -1;
194  tetPti = -1;
195  facei = -1;
196  }
197  }
198  else
199  {
200  coordinates = barycentric::uniform(NaN);
201  celli = -1;
202  tetFacei = -1;
203  tetPti = -1;
204  facei = -1;
205  }
206 }
207 
208 
209 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Templated 2D Barycentric derived from VectorSpace. Has 3 components, one of which is redundant.
Definition: Barycentric2D.H:53
const Cmpt & a() const
const Cmpt & b() const
const Cmpt & c() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
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.
void resize(const label)
Alias for setSize(const label)
Definition: PtrListI.H:32
Random number generator.
Definition: Random.H:58
scalar globalScalar01()
Advance the state and return a scalar sample from a uniform.
Definition: Random.C:94
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:101
Base class for patch-based injection models.
virtual ~patchInjectionBase()
Destructor.
const word patchName_
Patch name.
const label patchId_
Patch ID.
patchInjectionBase(const polyMesh &mesh, const word &patchName)
Construct from mesh and patch name.
virtual void setPositionAndCell(const fvMesh &mesh, Random &rnd, barycentric &coordinates, label &celli, label &tetFacei, label &tetPti, label &facei)
Set the injection position and owner cell, tetFace and tetPt.
virtual void topoChange(const polyMesh &mesh)
Update patch geometry and derived info for injection locations.
wordList names() const
Return a list of patch names.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1387
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:403
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:70
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:280
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:313
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:82
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the face for.
Definition: tetIndicesI.H:123
scalar mag() const
Return scalar magnitude.
Definition: triangleI.H:103
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
static Type NaN()
Return a primitive with all components set to NaN.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Barycentric< scalar > barycentric
A scalar version of the templated Barycentric.
Definition: barycentric.H:45
barycentric2D barycentric2D01(Random &rndGen)
Generate a random barycentric coordinate within the unit triangle.
Definition: barycentric2D.C:31
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
error FatalError
static const char nl
Definition: Ostream.H:260