checkTools.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) 2015-2017 OpenFOAM Foundation
9 Copyright (C) 2015-2022 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 "checkTools.H"
30#include "polyMesh.H"
31#include "globalMeshData.H"
32#include "hexMatcher.H"
33#include "wedgeMatcher.H"
34#include "prismMatcher.H"
35#include "pyrMatcher.H"
36#include "tetWedgeMatcher.H"
37#include "tetMatcher.H"
38#include "IOmanip.H"
39#include "OFstream.H"
40#include "pointSet.H"
41#include "faceSet.H"
42#include "cellSet.H"
43#include "Time.H"
44#include "coordSetWriter.H"
45#include "surfaceWriter.H"
46#include "syncTools.H"
47#include "globalIndex.H"
48#include "PatchTools.H"
49#include "functionObject.H"
50
51// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
52
53void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
54{
55 Info<< "Mesh stats " << mesh.regionName() << nl
56 << " points: "
57 << returnReduce(mesh.points().size(), sumOp<label>()) << nl;
58
59
60 // Count number of internal points (-1 if not sorted; 0 if no internal
61 // points)
62 const label minInt = returnReduce(mesh.nInternalPoints(), minOp<label>());
63 const label maxInt = returnReduce(mesh.nInternalPoints(), maxOp<label>());
64
65 if (minInt == -1 && maxInt > 0)
66 {
68 << "Some processors have their points sorted into internal"
69 << " and external and some do not." << endl
70 << " This can cause problems later on." << endl;
71 }
72 else if (minInt != -1)
73 {
74 // Assume all sorted
75 label nInternalPoints = returnReduce
76 (
77 mesh.nInternalPoints(),
78 sumOp<label>()
79 );
80 Info<< " internal points: " << nInternalPoints << nl;
81 }
82
83 if (allTopology && (minInt != -1))
84 {
85 label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
86 label nInternalEdges = returnReduce
87 (
88 mesh.nInternalEdges(),
89 sumOp<label>()
90 );
91 label nInternal1Edges = returnReduce
92 (
93 mesh.nInternal1Edges(),
94 sumOp<label>()
95 );
96 label nInternal0Edges = returnReduce
97 (
98 mesh.nInternal0Edges(),
99 sumOp<label>()
100 );
101
102 Info<< " edges: " << nEdges << nl
103 << " internal edges: " << nInternalEdges << nl
104 << " internal edges using one boundary point: "
105 << nInternal1Edges-nInternal0Edges << nl
106 << " internal edges using two boundary points: "
107 << nInternalEdges-nInternal1Edges << nl;
108 }
109
110 label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
111 label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
112 label nCells = returnReduce(mesh.cells().size(), sumOp<label>());
113 label nPatches = mesh.boundaryMesh().size();
114
115 Info<< " faces: " << nFaces << nl
116 << " internal faces: " << nIntFaces << nl
117 << " cells: " << nCells << nl
118 << " faces per cell: "
119 << (scalar(nFaces) + scalar(nIntFaces))/max(1, nCells) << nl
120 << " boundary patches: ";
121
122 if (Pstream::parRun())
123 {
124 // Number of global patches and min-max range of total patches
125 Info<< mesh.boundaryMesh().nNonProcessor() << ' '
126 << returnReduce(labelMinMax(nPatches), minMaxOp<label>()) << nl;
127 }
128 else
129 {
130 Info<< nPatches << nl;
131 }
132
133 Info<< " point zones: " << mesh.pointZones().size() << nl
134 << " face zones: " << mesh.faceZones().size() << nl
135 << " cell zones: " << mesh.cellZones().size() << nl
136 << endl;
137
138 // Construct shape recognizers
139 prismMatcher prism;
140 wedgeMatcher wedge;
141 tetWedgeMatcher tetWedge;
142
143 // Counters for different cell types
144 label nHex = 0;
145 label nWedge = 0;
146 label nPrism = 0;
147 label nPyr = 0;
148 label nTet = 0;
149 label nTetWedge = 0;
150 label nUnknown = 0;
151
152 Map<label> polyhedralFaces;
153
154 for (label celli = 0; celli < mesh.nCells(); celli++)
155 {
156 if (hexMatcher::test(mesh, celli))
157 {
158 nHex++;
159 }
160 else if (tetMatcher::test(mesh, celli))
161 {
162 nTet++;
163 }
164 else if (pyrMatcher::test(mesh, celli))
165 {
166 nPyr++;
167 }
168 else if (prism.isA(mesh, celli))
169 {
170 nPrism++;
171 }
172 else if (wedge.isA(mesh, celli))
173 {
174 nWedge++;
175 }
176 else if (tetWedge.isA(mesh, celli))
177 {
178 nTetWedge++;
179 }
180 else
181 {
182 nUnknown++;
183 polyhedralFaces(mesh.cells()[celli].size())++;
184 }
185 }
186
187 reduce(nHex,sumOp<label>());
188 reduce(nPrism,sumOp<label>());
189 reduce(nWedge,sumOp<label>());
190 reduce(nPyr,sumOp<label>());
191 reduce(nTetWedge,sumOp<label>());
192 reduce(nTet,sumOp<label>());
193 reduce(nUnknown,sumOp<label>());
194
195 Info<< "Overall number of cells of each type:" << nl
196 << " hexahedra: " << nHex << nl
197 << " prisms: " << nPrism << nl
198 << " wedges: " << nWedge << nl
199 << " pyramids: " << nPyr << nl
200 << " tet wedges: " << nTetWedge << nl
201 << " tetrahedra: " << nTet << nl
202 << " polyhedra: " << nUnknown
203 << endl;
204
205 if (nUnknown > 0)
206 {
207 Pstream::mapCombineGather(polyhedralFaces, plusEqOp<label>());
208
209 Info<< " Breakdown of polyhedra by number of faces:" << nl
210 << " faces" << " number of cells" << endl;
211
212 const labelList sortedKeys = polyhedralFaces.sortedToc();
213
214 forAll(sortedKeys, keyi)
215 {
216 const label nFaces = sortedKeys[keyi];
217
218 Info<< setf(std::ios::right) << setw(13)
219 << nFaces << " " << polyhedralFaces[nFaces] << nl;
220 }
221 }
222
223 Info<< endl;
224}
225
226
228(
229 const polyMesh& mesh,
230 surfaceWriter& writer,
231 const word& name,
232 const indirectPrimitivePatch& setPatch,
233 const fileName& outputDir
234)
235{
236 if (Pstream::parRun())
237 {
238 labelList pointToGlobal;
239 labelList uniqueMeshPointLabels;
240 autoPtr<globalIndex> globalPoints;
241 autoPtr<globalIndex> globalFaces;
242 faceList mergedFaces;
243 pointField mergedPoints;
245 (
246 mesh,
247 setPatch.localFaces(),
248 setPatch.meshPoints(),
249 setPatch.meshPointMap(),
250
251 pointToGlobal,
252 uniqueMeshPointLabels,
253 globalPoints,
254 globalFaces,
255
256 mergedFaces,
257 mergedPoints
258 );
259
260 // Write
261 if (Pstream::master())
262 {
263 writer.open
264 (
265 mergedPoints,
266 mergedFaces,
267 (outputDir / name),
268 false // serial - already merged
269 );
270
271 writer.write();
272 writer.clear();
273 }
274 }
275 else
276 {
277 writer.open
278 (
279 setPatch.localPoints(),
280 setPatch.localFaces(),
281 (outputDir / name),
282 false // serial - already merged
283 );
284
285 writer.write();
286 writer.clear();
287 }
288}
289
290
292(
293 surfaceWriter& writer,
294 const faceSet& set
295)
296{
297 const polyMesh& mesh = refCast<const polyMesh>(set.db());
298
299 const indirectPrimitivePatch setPatch
300 (
301 IndirectList<face>(mesh.faces(), set.sortedToc()),
302 mesh.points()
303 );
304
305 fileName outputDir
306 (
307 set.time().globalPath()
308 / functionObject::outputPrefix
309 / mesh.pointsInstance()
310 / set.name()
311 );
312 outputDir.clean(); // Remove unneeded ".."
313
314 mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
315}
316
317
319(
320 surfaceWriter& writer,
321 const cellSet& set
322)
323{
324 const polyMesh& mesh = refCast<const polyMesh>(set.db());
325 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
326
327
328 // Determine faces on outside of cellSet
329 bitSet isInSet(mesh.nCells());
330 for (const label celli : set)
331 {
332 isInSet.set(celli);
333 }
334
335
336 boolList bndInSet(mesh.nBoundaryFaces());
337 forAll(pbm, patchi)
338 {
339 const polyPatch& pp = pbm[patchi];
340 const labelList& fc = pp.faceCells();
341 forAll(fc, i)
342 {
343 bndInSet[pp.start()+i-mesh.nInternalFaces()] = isInSet[fc[i]];
344 }
345 }
346 syncTools::swapBoundaryFaceList(mesh, bndInSet);
347
348
349 DynamicList<label> outsideFaces(3*set.size());
350 for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
351 {
352 const bool ownVal = isInSet[mesh.faceOwner()[facei]];
353 const bool neiVal = isInSet[mesh.faceNeighbour()[facei]];
354
355 if (ownVal != neiVal)
356 {
357 outsideFaces.append(facei);
358 }
359 }
360
361
362 forAll(pbm, patchi)
363 {
364 const polyPatch& pp = pbm[patchi];
365 const labelList& fc = pp.faceCells();
366 if (pp.coupled())
367 {
368 forAll(fc, i)
369 {
370 label facei = pp.start()+i;
371
372 const bool neiVal = bndInSet[facei-mesh.nInternalFaces()];
373 if (isInSet[fc[i]] && !neiVal)
374 {
375 outsideFaces.append(facei);
376 }
377 }
378 }
379 else
380 {
381 forAll(fc, i)
382 {
383 if (isInSet[fc[i]])
384 {
385 outsideFaces.append(pp.start()+i);
386 }
387 }
388 }
389 }
390
391
392 const indirectPrimitivePatch setPatch
393 (
394 IndirectList<face>(mesh.faces(), outsideFaces),
395 mesh.points()
396 );
397
398 fileName outputDir
399 (
400 set.time().globalPath()
401 / functionObject::outputPrefix
402 / mesh.pointsInstance()
403 / set.name()
404 );
405 outputDir.clean(); // Remove unneeded ".."
406
407 mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
408}
409
410
412(
413 coordSetWriter& writer,
414 const pointSet& set
415)
416{
417 const polyMesh& mesh = refCast<const polyMesh>(set.db());
418
419 labelField mergedIDs(set.sortedToc());
420 pointField mergedPts(mesh.points(), mergedIDs);
421
422 if (Pstream::parRun())
423 {
424 // Note: we explicitly do not merge the points
425 // (mesh.globalData().mergePoints etc) since this might
426 // hide any synchronisation problem
427
428 // Renumber local ids -> global ids
429 globalIndex(mesh.nPoints()).inplaceToGlobal(mergedIDs);
430
431 globalIndex gatherer(mergedIDs.size(), globalIndex::gatherOnly{});
432 gatherer.gatherInplace(mergedIDs);
433 gatherer.gatherInplace(mergedPts);
434 }
435
436
437 // Write with pointID
438 if (Pstream::master())
439 {
440 coordSet coords(set.name(), "distance", mergedPts, mag(mergedPts));
441
442 // Output. E.g. pointSet p0 -> postProcessing/<time>/p0.vtk
443
444 fileName outputPath
445 (
446 set.time().globalPath()
447 / functionObject::outputPrefix
448 / mesh.pointsInstance()
449 / set.name()
450 );
451 outputPath.clean(); // Remove unneeded ".."
452
453 writer.open(coords, outputPath);
454 writer.nFields(1);
455 writer.write("pointID", mergedIDs);
456 writer.close(true);
457 }
458}
459
460
461// ************************************************************************* //
Istream and Ostream manipulators taking arguments.
Y[inertIndex] max(0.0)
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
reduce(hasMovingMesh, orOp< bool >())
static void gatherAndMerge(const scalar mergeDist, const PrimitivePatch< FaceList, PointField > &pp, Field< typename PrimitivePatch< FaceList, PointField >::point_type > &mergedPoints, List< typename PrimitivePatch< FaceList, PointField >::face_type > &mergedFaces, labelList &pointMergeMap, const bool useLocal=false)
Gather points and faces onto master and merge into single patch.
dynamicFvMesh & mesh
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const label nPatches
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
#define WarningInFunction
Report a warning using Foam::Warning.
void set(List< bool > &bools, const labelUList &locations)
Set the listed locations (assign 'true').
Definition: BitOps.C:38
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:169
List< label > labelList
A List of labels.
Definition: List.H:66
messageStream Info
Information stream (stdout output on master, null elsewhere)
MinMax< label > labelMinMax
A label min/max range.
Definition: MinMax.H:116
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void mergeAndWrite(const polyMesh &mesh, surfaceWriter &writer, const word &name, const indirectPrimitivePatch &setPatch, const fileName &outputDir)
Generate merged surface on master and write. Needs input patch.
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:52
List< bool > boolList
A List of bools.
Definition: List.H:64
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
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
void printMeshStats(const polyMesh &mesh, const bool allTopology)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333