channelIndex.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) 2011-2016 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 "channelIndex.H"
29#include "boolList.H"
30#include "syncTools.H"
31#include "OFstream.H"
32#include "meshTools.H"
33#include "Time.H"
34#include "SortableList.H"
35
36// * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * * //
37
38const Foam::Enum
39<
41>
42Foam::channelIndex::vectorComponentsNames_
43({
44 { vector::components::X, "x" },
45 { vector::components::Y, "y" },
46 { vector::components::Z, "z" },
47});
48
49
50// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
51
52// Determines face blocking
53void Foam::channelIndex::walkOppositeFaces
54(
55 const polyMesh& mesh,
56 const labelList& startFaces,
57 boolList& blockedFace
58)
59{
60 const cellList& cells = mesh.cells();
61 const faceList& faces = mesh.faces();
62 const label nBnd = mesh.nBoundaryFaces();
63
64 DynamicList<label> frontFaces(startFaces);
65 forAll(frontFaces, i)
66 {
67 label facei = frontFaces[i];
68 blockedFace[facei] = true;
69 }
70
71 while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
72 {
73 // Transfer across.
74 boolList isFrontBndFace(nBnd, false);
75 forAll(frontFaces, i)
76 {
77 label facei = frontFaces[i];
78
79 if (!mesh.isInternalFace(facei))
80 {
81 isFrontBndFace[facei-mesh.nInternalFaces()] = true;
82 }
83 }
85
86 // Add
87 forAll(isFrontBndFace, i)
88 {
89 label facei = mesh.nInternalFaces()+i;
90 if (isFrontBndFace[i] && !blockedFace[facei])
91 {
92 blockedFace[facei] = true;
93 frontFaces.append(facei);
94 }
95 }
96
97 // Transfer across cells
98 DynamicList<label> newFrontFaces(frontFaces.size());
99
100 forAll(frontFaces, i)
101 {
102 label facei = frontFaces[i];
103
104 {
105 const cell& ownCell = cells[mesh.faceOwner()[facei]];
106
107 label oppositeFacei = ownCell.opposingFaceLabel(facei, faces);
108
109 if (oppositeFacei == -1)
110 {
112 << "Face:" << facei << " owner cell:" << ownCell
113 << " is not a hex?" << abort(FatalError);
114 }
115 else
116 {
117 if (!blockedFace[oppositeFacei])
118 {
119 blockedFace[oppositeFacei] = true;
120 newFrontFaces.append(oppositeFacei);
121 }
122 }
123 }
124
125 if (mesh.isInternalFace(facei))
126 {
127 const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[facei]];
128
129 label oppositeFacei = neiCell.opposingFaceLabel(facei, faces);
130
131 if (oppositeFacei == -1)
132 {
134 << "Face:" << facei << " neighbour cell:" << neiCell
135 << " is not a hex?" << abort(FatalError);
136 }
137 else
138 {
139 if (!blockedFace[oppositeFacei])
140 {
141 blockedFace[oppositeFacei] = true;
142 newFrontFaces.append(oppositeFacei);
143 }
144 }
145 }
146 }
147
148 frontFaces.transfer(newFrontFaces);
149 }
150}
151
152
153// Calculate regions.
154void Foam::channelIndex::calcLayeredRegions
155(
156 const polyMesh& mesh,
157 const labelList& startFaces
158)
159{
160 boolList blockedFace(mesh.nFaces(), false);
161 walkOppositeFaces
162 (
163 mesh,
164 startFaces,
165 blockedFace
166 );
167
168
169 if (false)
170 {
171 OFstream str(mesh.time().path()/"blockedFaces.obj");
172 label vertI = 0;
173 forAll(blockedFace, facei)
174 {
175 if (blockedFace[facei])
176 {
177 const face& f = mesh.faces()[facei];
178 forAll(f, fp)
179 {
180 meshTools::writeOBJ(str, mesh.points()[f[fp]]);
181 }
182 str<< 'f';
183 forAll(f, fp)
184 {
185 str << ' ' << vertI+fp+1;
186 }
187 str << nl;
188 vertI += f.size();
189 }
190 }
191 }
192
193
194 // Do analysis for connected regions
195 cellRegion_.reset(new regionSplit(mesh, blockedFace));
196
197 Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
198
199 // Sum number of entries per region
200 regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
201
202 // Average cell centres to determine ordering.
203 pointField regionCc
204 (
206 / regionCount_
207 );
208
209 SortableList<scalar> sortComponent(regionCc.component(dir_));
210
211 sortMap_ = sortComponent.indices();
212
213 y_ = sortComponent;
214
215 if (symmetric_)
216 {
217 y_.setSize(cellRegion_().nRegions()/2);
218 }
219}
220
221
222// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
223
225(
226 const polyMesh& mesh,
227 const dictionary& dict
228)
229:
230 symmetric_(dict.get<bool>("symmetric")),
231 dir_(vectorComponentsNames_.get("component", dict))
232{
233 const polyBoundaryMesh& patches = mesh.boundaryMesh();
234
235 const wordList patchNames(dict.get<wordList>("patches"));
236
237 label nFaces = 0;
238
240 {
241 const label patchi = patches.findPatchID(patchNames[i]);
242
243 if (patchi == -1)
244 {
246 << "Illegal patch " << patchNames[i]
247 << ". Valid patches are " << patches.name()
248 << exit(FatalError);
249 }
250
251 nFaces += patches[patchi].size();
252 }
253
254 labelList startFaces(nFaces);
255 nFaces = 0;
256
258 {
259 const polyPatch& pp = patches[patchNames[i]];
260
261 forAll(pp, j)
262 {
263 startFaces[nFaces++] = pp.start()+j;
264 }
265 }
266
267 // Calculate regions.
268 calcLayeredRegions(mesh, startFaces);
269}
270
271
273(
274 const polyMesh& mesh,
275 const labelList& startFaces,
276 const bool symmetric,
277 const direction dir
278)
279:
280 symmetric_(symmetric),
281 dir_(dir)
282{
283 // Calculate regions.
284 calcLayeredRegions(mesh, startFaces);
285}
286
287
288// ************************************************************************* //
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
fileName path() const
Return path.
Definition: Time.H:358
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
components
Component labeling enumeration.
Definition: Vector.H:81
Does averaging of fields over layers of cells. Assumes layered mesh.
Definition: channelIndex.H:55
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
label nInternalFaces() const noexcept
Number of internal faces.
const vectorField & cellCentres() const
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const cellList & cells() const
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:445
bool
Definition: EEqn.H:20
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const cellShapeList & cells
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
List< word > wordList
A List of words.
Definition: fileName.H:63
List< label > labelList
A List of labels.
Definition: List.H:66
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static Map< Type > regionSum(const regionSplit &regions, const Field< Type > &fld)
List< bool > boolList
A List of bools.
Definition: List.H:64
error FatalError
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
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
wordList patchNames(nPatches)
labelList f(nPoints)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333