meshStructure.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) 2013-2016 OpenFOAM Foundation
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "meshStructure.H"
30#include "FaceCellWave.H"
31#include "topoDistanceData.H"
33#include "PointEdgeWave.H"
34#include "globalIndex.H"
35
36// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37
38namespace Foam
39{
41}
42
43
44// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45
46bool Foam::meshStructure::isStructuredCell
47(
48 const polyMesh& mesh,
49 const label layerI,
50 const label celli
51) const
52{
53 const cell& cFaces = mesh.cells()[celli];
54
55 // Count number of side faces
56 label nSide = 0;
57 forAll(cFaces, i)
58 {
59 if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
60 {
61 nSide++;
62 }
63 }
64
65 if (nSide != cFaces.size()-2)
66 {
67 return false;
68 }
69
70 // Check that side faces have correct point layers
71 forAll(cFaces, i)
72 {
73 if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
74 {
75 const face& f = mesh.faces()[cFaces[i]];
76
77 label nLayer = 0;
78 label nLayerPlus1 = 0;
79 forAll(f, fp)
80 {
81 label pointi = f[fp];
82 if (pointLayer_[pointi] == layerI)
83 {
84 nLayer++;
85 }
86 else if (pointLayer_[pointi] == layerI+1)
87 {
88 nLayerPlus1++;
89 }
90 }
91
92 if (f.size() != 4 || (nLayer+nLayerPlus1 != 4))
93 {
94 return false;
95 }
96 }
97 }
98
99 return true;
100}
101
102
104(
105 const polyMesh& mesh,
106 const uindirectPrimitivePatch& pp,
107 const globalIndex& globalFaces,
108 const globalIndex& globalEdges,
109 const globalIndex& globalPoints
110)
111{
112 // Field on cells and faces.
113 List<topoDistanceData<label>> cellData(mesh.nCells());
114 List<topoDistanceData<label>> faceData(mesh.nFaces());
115
116 {
117 if (debug)
118 {
119 Info<< typeName << " : seeding "
120 << returnReduce(pp.size(), sumOp<label>()) << " patch faces"
121 << nl << endl;
122 }
123
124
125 // Start of changes
126 labelList patchFaces(pp.size());
127 List<topoDistanceData<label>> patchData(pp.size());
128 forAll(pp, patchFacei)
129 {
130 patchFaces[patchFacei] = pp.addressing()[patchFacei];
131 patchData[patchFacei] = topoDistanceData<label>
132 (
133 0, // distance
134 globalFaces.toGlobal(patchFacei) // passive data
135 );
136 }
137
138
139 // Propagate information inwards
140 FaceCellWave<topoDistanceData<label>> distanceCalc
141 (
142 mesh,
143 patchFaces,
144 patchData,
145 faceData,
146 cellData,
148 );
149
150
151 // Determine cells from face-cell-walk
152 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
153
154 cellToPatchFaceAddressing_.setSize(mesh.nCells());
155 cellLayer_.setSize(mesh.nCells());
156 forAll(cellToPatchFaceAddressing_, celli)
157 {
158 cellToPatchFaceAddressing_[celli] = cellData[celli].data();
159 cellLayer_[celli] = cellData[celli].distance();
160 }
161
162
163
164 // Determine faces from face-cell-walk
165 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
166
167 faceToPatchFaceAddressing_.setSize(mesh.nFaces());
168 faceToPatchEdgeAddressing_.setSize(mesh.nFaces());
169 faceToPatchEdgeAddressing_ = labelMin;
170 faceLayer_.setSize(mesh.nFaces());
171
172 forAll(faceToPatchFaceAddressing_, facei)
173 {
174 label own = mesh.faceOwner()[facei];
175 label patchFacei = faceData[facei].data();
176 label patchDist = faceData[facei].distance();
177
178 if (mesh.isInternalFace(facei))
179 {
180 label nei = mesh.faceNeighbour()[facei];
181
182 if (cellData[own].distance() == cellData[nei].distance())
183 {
184 // side face
185 faceToPatchFaceAddressing_[facei] = 0;
186 faceLayer_[facei] = cellData[own].distance();
187 }
188 else if (cellData[own].distance() < cellData[nei].distance())
189 {
190 // unturned face
191 faceToPatchFaceAddressing_[facei] = patchFacei+1;
192 faceToPatchEdgeAddressing_[facei] = -1;
193 faceLayer_[facei] = patchDist;
194 }
195 else
196 {
197 // turned face
198 faceToPatchFaceAddressing_[facei] = -(patchFacei+1);
199 faceToPatchEdgeAddressing_[facei] = -1;
200 faceLayer_[facei] = patchDist;
201 }
202 }
203 else if (patchDist == cellData[own].distance())
204 {
205 // starting face
206 faceToPatchFaceAddressing_[facei] = -(patchFacei+1);
207 faceToPatchEdgeAddressing_[facei] = -1;
208 faceLayer_[facei] = patchDist;
209 }
210 else
211 {
212 // unturned face or side face. Cannot be determined until
213 // we determine the point layers. Problem is that both are
214 // the same number of steps away from the initial seed face.
215 }
216 }
217 }
218
219
220 // Determine points from separate walk on point-edge
221 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
222
223 {
224 pointToPatchPointAddressing_.setSize(mesh.nPoints());
225 pointLayer_.setSize(mesh.nPoints());
226
227 if (debug)
228 {
229 Info<< typeName << " : seeding "
230 << returnReduce(pp.nPoints(), sumOp<label>()) << " patch points"
231 << nl << endl;
232 }
233
234 // Field on edges and points.
235 List<pointTopoDistanceData<label>> edgeData(mesh.nEdges());
236 List<pointTopoDistanceData<label>> pointData(mesh.nPoints());
237
238 // Start of changes
239 labelList patchPoints(pp.nPoints());
240 List<pointTopoDistanceData<label>> patchData(pp.nPoints());
241 forAll(pp.meshPoints(), patchPointi)
242 {
243 patchPoints[patchPointi] = pp.meshPoints()[patchPointi];
244 patchData[patchPointi] = pointTopoDistanceData<label>
245 (
246 0, // distance
247 globalPoints.toGlobal(patchPointi) // passive data
248 );
249 }
250
251
252 // Walk
253 PointEdgeWave<pointTopoDistanceData<label>> distanceCalc
254 (
255 mesh,
256 patchPoints,
257 patchData,
258
259 pointData,
260 edgeData,
261 mesh.globalData().nTotalPoints() // max iterations
262 );
263
264 forAll(pointData, pointi)
265 {
266 pointToPatchPointAddressing_[pointi] = pointData[pointi].data();
267 pointLayer_[pointi] = pointData[pointi].distance();
268 }
269
270
271 // Derive from originating patch points what the patch edges were.
272 EdgeMap<label> pointsToEdge(pp.nEdges());
273 forAll(pp.edges(), edgeI)
274 {
275 const edge& e = pp.edges()[edgeI];
276 edge globalEdge
277 (
278 globalPoints.toGlobal(e[0]),
279 globalPoints.toGlobal(e[1])
280 );
281 pointsToEdge.insert(globalEdge, globalEdges.toGlobal(edgeI));
282 }
283
284 // Look up on faces
285 forAll(faceToPatchEdgeAddressing_, facei)
286 {
287 if (faceToPatchEdgeAddressing_[facei] == labelMin)
288 {
289 // Face not yet done. Check if all points on same level
290 // or if not see what edge it originates from
291
292 const face& f = mesh.faces()[facei];
293
294 label levelI = pointLayer_[f[0]];
295 for (label fp = 1; fp < f.size(); fp++)
296 {
297 if (pointLayer_[f[fp]] != levelI)
298 {
299 levelI = -1;
300 break;
301 }
302 }
303
304 if (levelI != -1)
305 {
306 // All same level
307 //Pout<< "Horizontal boundary face " << facei
308 // << " at:" << mesh.faceCentres()[facei]
309 // << " data:" << faceData[facei]
310 // << " pointDatas:"
311 // << UIndirectList<pointTopoDistanceData<label>>
312 // (pointData, f)
313 // << endl;
314
315 label patchFacei = faceData[facei].data();
316 label patchDist = faceData[facei].distance();
317
318 faceToPatchEdgeAddressing_[facei] = -1;
319 faceToPatchFaceAddressing_[facei] = patchFacei+1;
320 faceLayer_[facei] = patchDist;
321 }
322 else
323 {
324 // Points of face on different levels
325
326 // See if there is any edge
327 forAll(f, fp)
328 {
329 label pointi = f[fp];
330 label nextPointi = f.nextLabel(fp);
331
332 const auto fnd = pointsToEdge.cfind
333 (
334 edge
335 (
336 pointData[pointi].data(),
337 pointData[nextPointi].data()
338 )
339 );
340
341 if (fnd.found())
342 {
343 faceToPatchEdgeAddressing_[facei] = fnd.val();
344 faceToPatchFaceAddressing_[facei] = 0;
345 label own = mesh.faceOwner()[facei];
346 faceLayer_[facei] = cellData[own].distance();
347
348 // Note: could test whether the other edges on the
349 // face are consistent
350 break;
351 }
352 }
353 }
354 }
355 }
356 }
357
358
359
360 // Use maps to find out mesh structure.
361 {
362 label nLayers = gMax(cellLayer_)+1;
363 labelListList layerToCells(invertOneToMany(nLayers, cellLayer_));
364
365 structured_ = true;
366 forAll(layerToCells, layerI)
367 {
368 const labelList& lCells = layerToCells[layerI];
369
370 forAll(lCells, lCelli)
371 {
372 label celli = lCells[lCelli];
373
374 structured_ = isStructuredCell
375 (
376 mesh,
377 layerI,
378 celli
379 );
380
381 if (!structured_)
382 {
383 break;
384 }
385 }
386
387 if (!structured_)
388 {
389 break;
390 }
391 }
392
393 reduce(structured_, andOp<bool>());
394 }
395}
396
397
398// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
399
401(
402 const polyMesh& mesh,
404)
405{
406 correct
407 (
408 mesh,
409 pp,
410 globalIndex(pp.size()),
411 globalIndex(pp.nEdges()),
412 globalIndex(pp.nPoints())
413 );
414}
415
416
418(
419 const polyMesh& mesh,
420 const uindirectPrimitivePatch& pp,
421 const globalIndex& globalFaces,
422 const globalIndex& globalEdges,
424)
425{
426 correct
427 (
428 mesh,
429 pp,
430 globalFaces,
431 globalEdges,
433 );
434}
435
436
437// ************************************************************************* //
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A list of faces which address into the list of points.
label nEdges() const
Number of edges in patch.
label nPoints() const
Number of points supporting patch faces.
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
T * data() noexcept
Return pointer to the underlying array serving as data storage.
Definition: UListI.H:237
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
label nTotalPoints() const noexcept
Return total number of points in decomposed mesh. Not.
label nTotalCells() const noexcept
Return total number of cells in decomposed mesh.
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:103
Detect extruded mesh structure given a set of faces (uindirectPrimitivePatch).
Definition: meshStructure.H:62
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1310
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const cellList & cells() const
label nEdges() const
Number of mesh edges.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
thermo correct()
dynamicFvMesh & mesh
Namespace for OpenFOAM.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
List< label > labelList
A List of labels.
Definition: List.H:66
constexpr label labelMin
Definition: label.H:60
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:114
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333