checkTools.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) 2015-2018 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 "checkTools.H"
27 #include "polyMesh.H"
28 #include "globalMeshData.H"
29 #include "hexMatcher.H"
30 #include "wedgeMatcher.H"
31 #include "prismMatcher.H"
32 #include "pyrMatcher.H"
33 #include "tetWedgeMatcher.H"
34 #include "tetMatcher.H"
35 #include "IOmanip.H"
36 #include "pointSet.H"
37 #include "faceSet.H"
38 #include "cellSet.H"
39 #include "Time.H"
40 #include "surfaceWriter.H"
41 #include "syncTools.H"
42 #include "globalIndex.H"
43 #include "PatchTools.H"
44 
45 
46 void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
47 {
48  Info<< "Mesh stats" << nl
49  << " points: "
50  << returnReduce(mesh.points().size(), sumOp<label>()) << nl;
51 
52  label nInternalPoints = returnReduce
53  (
54  mesh.nInternalPoints(),
55  sumOp<label>()
56  );
57 
58  if (nInternalPoints != -Pstream::nProcs())
59  {
60  Info<< " internal points: " << nInternalPoints << nl;
61 
62  if (returnReduce(mesh.nInternalPoints(), minOp<label>()) == -1)
63  {
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  }
70 
71  if (allTopology && nInternalPoints != -Pstream::nProcs())
72  {
73  label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
74  label nInternalEdges = returnReduce
75  (
76  mesh.nInternalEdges(),
77  sumOp<label>()
78  );
79  label nInternal1Edges = returnReduce
80  (
81  mesh.nInternal1Edges(),
82  sumOp<label>()
83  );
84  label nInternal0Edges = returnReduce
85  (
86  mesh.nInternal0Edges(),
87  sumOp<label>()
88  );
89 
90  Info<< " edges: " << nEdges << nl
91  << " internal edges: " << nInternalEdges << nl
92  << " internal edges using one boundary point: "
93  << nInternal1Edges-nInternal0Edges << nl
94  << " internal edges using two boundary points: "
95  << nInternalEdges-nInternal1Edges << nl;
96  }
97 
98  label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
99  label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
100  label nCells = returnReduce(mesh.cells().size(), sumOp<label>());
101 
102  Info<< " faces: " << nFaces << nl
103  << " internal faces: " << nIntFaces << nl
104  << " cells: " << nCells << nl
105  << " faces per cell: "
106  << scalar(nFaces + nIntFaces)/max(1, nCells) << nl
107  << " boundary patches: " << mesh.boundaryMesh().size() << nl
108  << " point zones: " << mesh.pointZones().size() << nl
109  << " face zones: " << mesh.faceZones().size() << nl
110  << " cell zones: " << mesh.cellZones().size() << nl
111  << endl;
112 
113  // Construct shape recognizers
114  hexMatcher hex;
115  prismMatcher prism;
116  wedgeMatcher wedge;
117  pyrMatcher pyr;
118  tetWedgeMatcher tetWedge;
119  tetMatcher tet;
120 
121  // Counters for different cell types
122  label nHex = 0;
123  label nWedge = 0;
124  label nPrism = 0;
125  label nPyr = 0;
126  label nTet = 0;
127  label nTetWedge = 0;
128  label nUnknown = 0;
129 
130  Map<label> polyhedralFaces;
131 
132  for (label celli = 0; celli < mesh.nCells(); celli++)
133  {
134  if (hex.isA(mesh, celli))
135  {
136  nHex++;
137  }
138  else if (tet.isA(mesh, celli))
139  {
140  nTet++;
141  }
142  else if (pyr.isA(mesh, celli))
143  {
144  nPyr++;
145  }
146  else if (prism.isA(mesh, celli))
147  {
148  nPrism++;
149  }
150  else if (wedge.isA(mesh, celli))
151  {
152  nWedge++;
153  }
154  else if (tetWedge.isA(mesh, celli))
155  {
156  nTetWedge++;
157  }
158  else
159  {
160  nUnknown++;
161  polyhedralFaces(mesh.cells()[celli].size())++;
162  }
163  }
164 
165  reduce(nHex,sumOp<label>());
166  reduce(nPrism,sumOp<label>());
167  reduce(nWedge,sumOp<label>());
168  reduce(nPyr,sumOp<label>());
169  reduce(nTetWedge,sumOp<label>());
170  reduce(nTet,sumOp<label>());
171  reduce(nUnknown,sumOp<label>());
172 
173  Info<< "Overall number of cells of each type:" << nl
174  << " hexahedra: " << nHex << nl
175  << " prisms: " << nPrism << nl
176  << " wedges: " << nWedge << nl
177  << " pyramids: " << nPyr << nl
178  << " tet wedges: " << nTetWedge << nl
179  << " tetrahedra: " << nTet << nl
180  << " polyhedra: " << nUnknown
181  << endl;
182 
183  if (nUnknown > 0)
184  {
185  Pstream::mapCombineGather(polyhedralFaces, plusEqOp<label>());
186 
187  Info<< " Breakdown of polyhedra by number of faces:" << nl
188  << " faces" << " number of cells" << endl;
189 
190  const labelList sortedKeys = polyhedralFaces.sortedToc();
191 
192  forAll(sortedKeys, keyI)
193  {
194  const label nFaces = sortedKeys[keyI];
195 
196  Info<< setf(std::ios::right) << setw(13)
197  << nFaces << " " << polyhedralFaces[nFaces] << nl;
198  }
199  }
200 
201  Info<< endl;
202 }
203 
204 
206 (
207  const polyMesh& mesh,
208  const surfaceWriter& writer,
209  const word& name,
210  const indirectPrimitivePatch setPatch,
211  const fileName& outputDir
212 )
213 {
214  if (Pstream::parRun())
215  {
216  labelList pointToGlobal;
217  labelList uniqueMeshPointLabels;
218  autoPtr<globalIndex> globalPoints;
219  autoPtr<globalIndex> globalFaces;
220  faceList mergedFaces;
221  pointField mergedPoints;
223  (
224  mesh,
225  setPatch.localFaces(),
226  setPatch.meshPoints(),
227  setPatch.meshPointMap(),
228 
229  pointToGlobal,
230  uniqueMeshPointLabels,
231  globalPoints,
232  globalFaces,
233 
234  mergedFaces,
235  mergedPoints
236  );
237 
238  // Write
239  if (Pstream::master())
240  {
241  writer.write(outputDir, name, mergedPoints, mergedFaces);
242  }
243  }
244  else
245  {
246  writer.write
247  (
248  outputDir,
249  name,
250  setPatch.localPoints(),
251  setPatch.localFaces()
252  );
253  }
254 }
255 
256 
258 (
259  const surfaceWriter& writer,
260  const faceSet& set
261 )
262 {
263  const polyMesh& mesh = refCast<const polyMesh>(set.db());
264 
265  const indirectPrimitivePatch setPatch
266  (
267  IndirectList<face>(mesh.faces(), set.sortedToc()),
268  mesh.points()
269  );
270 
271  fileName outputDir
272  (
273  set.time().path()
274  / (Pstream::parRun() ? ".." : "")
275  / "postProcessing"
276  / mesh.pointsInstance()
277  / set.name()
278  );
279  outputDir.clean();
280 
281  mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
282 }
283 
284 
286 (
287  const surfaceWriter& writer,
288  const cellSet& set
289 )
290 {
291  const polyMesh& mesh = refCast<const polyMesh>(set.db());
292  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
293 
294 
295  // Determine faces on outside of cellSet
296  PackedBoolList isInSet(mesh.nCells());
297  forAllConstIter(cellSet, set, iter)
298  {
299  isInSet[iter.key()] = true;
300  }
301 
302 
303  boolList bndInSet(mesh.nFaces()-mesh.nInternalFaces());
304  forAll(pbm, patchi)
305  {
306  const polyPatch& pp = pbm[patchi];
307  const labelList& fc = pp.faceCells();
308  forAll(fc, i)
309  {
310  bndInSet[pp.start()+i-mesh.nInternalFaces()] = isInSet[fc[i]];
311  }
312  }
313  syncTools::swapBoundaryFaceList(mesh, bndInSet);
314 
315 
316  DynamicList<label> outsideFaces(3*set.size());
317  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
318  {
319  bool ownVal = isInSet[mesh.faceOwner()[facei]];
320  bool neiVal = isInSet[mesh.faceNeighbour()[facei]];
321 
322  if (ownVal != neiVal)
323  {
324  outsideFaces.append(facei);
325  }
326  }
327 
328 
329  forAll(pbm, patchi)
330  {
331  const polyPatch& pp = pbm[patchi];
332  const labelList& fc = pp.faceCells();
333  if (pp.coupled())
334  {
335  forAll(fc, i)
336  {
337  label facei = pp.start()+i;
338 
339  bool neiVal = bndInSet[facei-mesh.nInternalFaces()];
340  if (isInSet[fc[i]] && !neiVal)
341  {
342  outsideFaces.append(facei);
343  }
344  }
345  }
346  else
347  {
348  forAll(fc, i)
349  {
350  if (isInSet[fc[i]])
351  {
352  outsideFaces.append(pp.start()+i);
353  }
354  }
355  }
356  }
357 
358 
359  const indirectPrimitivePatch setPatch
360  (
361  IndirectList<face>(mesh.faces(), outsideFaces),
362  mesh.points()
363  );
364 
365  fileName outputDir
366  (
367  set.time().path()
368  / (Pstream::parRun() ? ".." : "")
369  / "postProcessing"
370  / mesh.pointsInstance()
371  / set.name()
372  );
373  outputDir.clean();
374 
375  mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
376 }
377 
378 
380 (
381  const writer<scalar>& writer,
382  const pointSet& set
383 )
384 {
385  const polyMesh& mesh = refCast<const polyMesh>(set.db());
386 
387  pointField mergedPts;
388  labelList mergedIDs;
389 
390  if (Pstream::parRun())
391  {
392  // Note: we explicitly do not merge the points
393  // (mesh.globalData().mergePoints etc) since this might
394  // hide any synchronisation problem
395 
396  globalIndex globalNumbering(mesh.nPoints());
397 
398  mergedPts.setSize(returnReduce(set.size(), sumOp<label>()));
399  mergedIDs.setSize(mergedPts.size());
400 
401  labelList setPointIDs(set.sortedToc());
402 
403  // Get renumbered local data
404  pointField myPoints(mesh.points(), setPointIDs);
405  labelList myIDs(setPointIDs.size());
406  forAll(setPointIDs, i)
407  {
408  myIDs[i] = globalNumbering.toGlobal(setPointIDs[i]);
409  }
410 
411  if (Pstream::master())
412  {
413  // Insert master data first
414  label pOffset = 0;
415  SubList<point>(mergedPts, myPoints.size(), pOffset) = myPoints;
416  SubList<label>(mergedIDs, myIDs.size(), pOffset) = myIDs;
417  pOffset += myPoints.size();
418 
419  // Receive slave ones
420  for (int slave=1; slave<Pstream::nProcs(); slave++)
421  {
422  IPstream fromSlave(Pstream::commsTypes::scheduled, slave);
423 
424  pointField slavePts(fromSlave);
425  labelList slaveIDs(fromSlave);
426 
427  SubList<point>(mergedPts, slavePts.size(), pOffset) = slavePts;
428  SubList<label>(mergedIDs, slaveIDs.size(), pOffset) = slaveIDs;
429  pOffset += slaveIDs.size();
430  }
431  }
432  else
433  {
434  // Construct processor stream with estimate of size. Could
435  // be improved.
436  OPstream toMaster
437  (
438  Pstream::commsTypes::scheduled,
439  Pstream::masterNo(),
440  myPoints.byteSize() + myIDs.byteSize()
441  );
442  toMaster << myPoints << myIDs;
443  }
444  }
445  else
446  {
447  mergedIDs = set.sortedToc();
448  mergedPts = pointField(mesh.points(), mergedIDs);
449  }
450 
451 
452  // Write with scalar pointID
453  if (Pstream::master())
454  {
455  scalarField scalarPointIDs(mergedIDs.size());
456  forAll(mergedIDs, i)
457  {
458  scalarPointIDs[i] = 1.0*mergedIDs[i];
459  }
460 
461  coordSet points(set.name(), "distance", mergedPts, mag(mergedPts));
462 
463  List<const scalarField*> flds(1, &scalarPointIDs);
464 
465  wordList fldNames(1, "pointID");
466 
467  // Output e.g. pointSet p0 to
468  // postProcessing/<time>/p0.vtk
469  fileName outputDir
470  (
471  set.time().path()
472  / (Pstream::parRun() ? ".." : "")
473  / "postProcessing"
474  / mesh.pointsInstance()
475  // set.name()
476  );
477  outputDir.clean();
478  mkDir(outputDir);
479 
480  fileName outputFile(outputDir/writer.getFileName(points, wordList()));
481  // fileName outputFile(outputDir/set.name());
482 
483  OFstream os(outputFile);
484 
485  writer.write(points, fldNames, flds, os);
486  }
487 }
488 
489 
490 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
IOstream & hex(IOstream &io)
Definition: IOstream.H:564
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
const pointField & points
static void gatherAndMerge(const scalar mergeDist, const PrimitivePatch< FaceList, PointField > &p, Field< typename PrimitivePatch< FaceList, PointField >::PointType > &mergedPoints, List< typename PrimitivePatch< FaceList, PointField >::FaceType > &mergedFaces, labelList &pointMergeMap)
Gather points and faces onto master and merge into single patch.
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:164
List< label > labelList
A List of labels.
Definition: labelList.H:56
void mergeAndWrite(const polyMesh &mesh, const surfaceWriter &writer, const word &name, const indirectPrimitivePatch setPatch, const fileName &outputDir)
Generate merged surface on master and write. Needs input patch.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
volScalarField scalarField(fieldObject, mesh)
Istream and Ostream manipulators taking arguments.
static const char nl
Definition: Ostream.H:265
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
Definition: POSIX.C:290
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void printMeshStats(const polyMesh &mesh, const bool allTopology)
List< word > wordList
A List of words.
Definition: fileName.H:54
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
fileName path(UMean.rootPath()/UMean.caseName()/functionObjects::writeFile::outputPrefix/"graphs"/UMean.instance())