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-2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
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 "pointSet.H"
40 #include "faceSet.H"
41 #include "cellSet.H"
42 #include "Time.H"
43 #include "surfaceWriter.H"
44 #include "syncTools.H"
45 #include "globalIndex.H"
46 #include "PatchTools.H"
47 #include "functionObject.H"
48 
49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50 
51 void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
52 {
53  Info<< "Mesh stats" << nl
54  << " points: "
55  << returnReduce(mesh.points().size(), sumOp<label>()) << nl;
56 
57  // Count number of internal points (-1 if not sorted; 0 if no internal
58  // points)
59  const label minInt = returnReduce(mesh.nInternalPoints(), minOp<label>());
60  const label maxInt = returnReduce(mesh.nInternalPoints(), maxOp<label>());
61 
62  if (minInt == -1 && maxInt > 0)
63  {
64  WarningInFunction
65  << "Some processors have their points sorted into internal"
66  << " and external and some do not." << endl
67  << " This can cause problems later on." << endl;
68  }
69  else if (minInt != -1)
70  {
71  // Assume all sorted
72  label nInternalPoints = returnReduce
73  (
74  mesh.nInternalPoints(),
75  sumOp<label>()
76  );
77  Info<< " internal points: " << nInternalPoints << nl;
78  }
79 
80  if (allTopology && (minInt != -1))
81  {
82  label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
83  label nInternalEdges = returnReduce
84  (
85  mesh.nInternalEdges(),
86  sumOp<label>()
87  );
88  label nInternal1Edges = returnReduce
89  (
90  mesh.nInternal1Edges(),
91  sumOp<label>()
92  );
93  label nInternal0Edges = returnReduce
94  (
95  mesh.nInternal0Edges(),
96  sumOp<label>()
97  );
98 
99  Info<< " edges: " << nEdges << nl
100  << " internal edges: " << nInternalEdges << nl
101  << " internal edges using one boundary point: "
102  << nInternal1Edges-nInternal0Edges << nl
103  << " internal edges using two boundary points: "
104  << nInternalEdges-nInternal1Edges << nl;
105  }
106 
107  label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
108  label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
109  label nCells = returnReduce(mesh.cells().size(), sumOp<label>());
110  label nPatches = mesh.boundaryMesh().size();
111 
112  Info<< " faces: " << nFaces << nl
113  << " internal faces: " << nIntFaces << nl
114  << " cells: " << nCells << nl
115  << " faces per cell: "
116  << (scalar(nFaces) + scalar(nIntFaces))/max(1, nCells) << nl
117  << " boundary patches: ";
118 
119  if (Pstream::parRun())
120  {
121  // Number of global patches and min-max range of total patches
122  Info<< mesh.boundaryMesh().nNonProcessor() << ' '
123  << returnReduce(labelMinMax(nPatches), minMaxOp<label>()) << nl;
124  }
125  else
126  {
127  Info<< nPatches << nl;
128  }
129 
130  Info<< " point zones: " << mesh.pointZones().size() << nl
131  << " face zones: " << mesh.faceZones().size() << nl
132  << " cell zones: " << mesh.cellZones().size() << nl
133  << endl;
134 
135  // Construct shape recognizers
136  hexMatcher hex;
137  prismMatcher prism;
138  wedgeMatcher wedge;
139  pyrMatcher pyr;
140  tetWedgeMatcher tetWedge;
141  tetMatcher tet;
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 (hex.isA(mesh, celli))
157  {
158  nHex++;
159  }
160  else if (tet.isA(mesh, celli))
161  {
162  nTet++;
163  }
164  else if (pyr.isA(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 
227 void Foam::mergeAndWrite
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;
244  Foam::PatchTools::gatherAndMerge
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 
291 void Foam::mergeAndWrite
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();
313 
314  mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
315 }
316 
317 
318 void Foam::mergeAndWrite
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();
406 
407  mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
408 }
409 
410 
411 void Foam::mergeAndWrite
412 (
413  const writer<scalar>& writer,
414  const pointSet& set
415 )
416 {
417  const polyMesh& mesh = refCast<const polyMesh>(set.db());
418 
419  pointField mergedPts;
420  labelList 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  globalIndex globalNumbering(mesh.nPoints());
429 
430  mergedPts.setSize(returnReduce(set.size(), sumOp<label>()));
431  mergedIDs.setSize(mergedPts.size());
432 
433  labelList setPointIDs(set.sortedToc());
434 
435  // Get renumbered local data
436  pointField myPoints(mesh.points(), setPointIDs);
437  labelList myIDs(globalNumbering.toGlobal(setPointIDs));
438 
439  if (Pstream::master())
440  {
441  // Insert master data first
442  label pOffset = 0;
443  SubList<point>(mergedPts, myPoints.size(), pOffset) = myPoints;
444  SubList<label>(mergedIDs, myIDs.size(), pOffset) = myIDs;
445  pOffset += myPoints.size();
446 
447  // Receive slave ones
448  for (int slave=1; slave<Pstream::nProcs(); slave++)
449  {
450  IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
451 
452  pointField slavePts(fromSlave);
453  labelList slaveIDs(fromSlave);
454 
455  SubList<point>(mergedPts, slavePts.size(), pOffset) = slavePts;
456  SubList<label>(mergedIDs, slaveIDs.size(), pOffset) = slaveIDs;
457  pOffset += slaveIDs.size();
458  }
459  }
460  else
461  {
462  // Construct processor stream with estimate of size. Could
463  // be improved.
464  OPstream toMaster
465  (
466  Pstream::commsTypes::scheduled,
467  Pstream::masterNo(),
468  myPoints.byteSize() + myIDs.byteSize()
469  );
470  toMaster << myPoints << myIDs;
471  }
472  }
473  else
474  {
475  mergedIDs = set.sortedToc();
476  mergedPts = pointField(mesh.points(), mergedIDs);
477  }
478 
479 
480  // Write with scalar pointID
481  if (Pstream::master())
482  {
483  scalarField scalarPointIDs(mergedIDs.size());
484  forAll(mergedIDs, i)
485  {
486  scalarPointIDs[i] = 1.0*mergedIDs[i];
487  }
488 
489  coordSet points(set.name(), "distance", mergedPts, mag(mergedPts));
490 
491  List<const scalarField*> flds(1, &scalarPointIDs);
492 
493  wordList fldNames(1, "pointID");
494 
495  // Output e.g. pointSet p0 to
496  // postProcessing/<time>/p0.vtk
497  fileName outputDir
498  (
499  set.time().globalPath()
500  / functionObject::outputPrefix
501  / mesh.pointsInstance()
502  // set.name()
503  );
504  outputDir.clean();
505  mkDir(outputDir);
506 
507  fileName outputFile(outputDir/writer.getFileName(points, wordList()));
508  //fileName outputFile(outputDir/set.name());
509 
510  OFstream os(outputFile);
511 
512  writer.write(points, fldNames, flds, os);
513  }
514 }
515 
516 
517 // ************************************************************************* //
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:71
Foam::setf
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:169
hex
const cellModel & hex
Definition: createBlockMesh.H:1
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
wedgeMatcher.H
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
nPatches
label nPatches
Definition: readKivaGrid.H:396
PatchTools.H
globalMeshData.H
Foam::PatchTools::gatherAndMerge
static void gatherAndMerge(const scalar mergeDist, const PrimitivePatch< FaceList, PointField > &p, Field< typename PrimitivePatch< FaceList, PointField >::point_type > &mergedPoints, List< typename PrimitivePatch< FaceList, PointField >::face_type > &mergedFaces, labelList &pointMergeMap)
Gather points and faces onto master and merge into single patch.
Definition: PatchToolsGatherAndMerge.C:38
globalIndex.H
hexMatcher.H
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:69
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
surfaceWriter.H
tetWedgeMatcher.H
polyMesh.H
syncTools.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::wordList
List< word > wordList
A List of words.
Definition: fileName.H:59
prismMatcher.H
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
Foam::name
word name(const complex &c)
Return string representation of complex.
Definition: complex.C:76
faceSet.H
IOmanip.H
Istream and Ostream manipulators taking arguments.
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
reduce
reduce(hasMovingMesh, orOp< bool >())
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::labelMinMax
MinMax< label > labelMinMax
A label min/max range.
Definition: MinMax.H:113
Foam::setw
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
Foam::indirectPrimitivePatch
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
Definition: indirectPrimitivePatch.H:49
Time.H
Foam::nl
constexpr char nl
Definition: Ostream.H:385
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::printMeshStats
void printMeshStats(const polyMesh &mesh, const bool allTopology)
tetMatcher.H
Foam::mergeAndWrite
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.
cellSet.H
Foam::List::setSize
void setSize(const label newSize)
Alias for resize(const label)
Definition: ListI.H:146
functionObject.H
Foam::mkDir
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
Definition: MSwindows.C:507
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:298
checkTools.H
pyrMatcher.H
pointSet.H