checkTopology.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 Copyright (C) 2017-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 "checkTopology.H"
30#include "polyMesh.H"
31#include "Time.H"
32#include "regionSplit.H"
33#include "cellSet.H"
34#include "faceSet.H"
35#include "pointSet.H"
36#include "IOmanip.H"
37#include "emptyPolyPatch.H"
38#include "processorPolyPatch.H"
39#include "vtkCoordSetWriter.H"
40#include "vtkSurfaceWriter.H"
41#include "checkTools.H"
42#include "treeBoundBox.H"
43
44// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
45
46template<class PatchType>
48(
49 const bool allGeometry,
50 const word& name,
51 const PatchType& pp,
52 pointSet& points
53)
54{
55 Info<< " "
56 << setw(20) << name
57 << setw(9) << returnReduce(pp.size(), sumOp<label>())
58 << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
59
60 if (!Pstream::parRun())
61 {
62 typedef typename PatchType::surfaceTopo TopoType;
63 TopoType pTyp = pp.surfaceType();
64
65 if (pp.empty())
66 {
67 Info<< setw(34) << "ok (empty)";
68 }
69 else if (pTyp == TopoType::MANIFOLD)
70 {
71 if (pp.checkPointManifold(true, &points))
72 {
73 Info<< setw(34)
74 << "multiply connected (shared point)";
75 }
76 else
77 {
78 Info<< setw(34) << "ok (closed singly connected)";
79 }
80
81 // Add points on non-manifold edges to make set complete
82 pp.checkTopology(false, &points);
83 }
84 else
85 {
86 pp.checkTopology(false, &points);
87
88 if (pTyp == TopoType::OPEN)
89 {
90 Info<< setw(34)
91 << "ok (non-closed singly connected)";
92 }
93 else
94 {
95 Info<< setw(34)
96 << "multiply connected (shared edge)";
97 }
98 }
99 }
100
101 if (allGeometry)
102 {
103 const labelList& mp = pp.meshPoints();
104
105 if (returnReduce(mp.size(), sumOp<label>()) > 0)
106 {
107 boundBox bb(pp.points(), mp, true); // reduce
108 Info<< ' ' << bb;
109 }
110 }
111}
112
113
114Foam::label Foam::checkTopology
115(
116 const polyMesh& mesh,
117 const bool allTopology,
118 const bool allGeometry,
119 autoPtr<surfaceWriter>& surfWriter,
120 autoPtr<coordSetWriter>& setWriter
121)
122{
123 label noFailedChecks = 0;
124
125 Info<< "Checking topology..." << endl;
126
127 // Check if the boundary definition is unique
128 mesh.boundaryMesh().checkDefinition(true);
129
130 // Check that empty patches cover all sides of the mesh
131 {
132 label nEmpty = 0;
133 forAll(mesh.boundaryMesh(), patchi)
134 {
135 if (isA<emptyPolyPatch>(mesh.boundaryMesh()[patchi]))
136 {
137 nEmpty += mesh.boundaryMesh()[patchi].size();
138 }
139 }
140 reduce(nEmpty, sumOp<label>());
141 label nTotCells = returnReduce(mesh.cells().size(), sumOp<label>());
142
143 // These are actually warnings, not errors.
144 if (nTotCells && (nEmpty % nTotCells))
145 {
146 Info<< " ***Total number of faces on empty patches"
147 << " is not divisible by the number of cells in the mesh."
148 << " Hence this mesh is not 1D or 2D."
149 << endl;
150 }
151 }
152
153 // Check if the boundary processor patches are correct
154 mesh.boundaryMesh().checkParallelSync(true);
155
156 // Check names of zones are equal
157 mesh.cellZones().checkDefinition(true);
158 if (mesh.cellZones().checkParallelSync(true))
159 {
160 noFailedChecks++;
161 }
162 mesh.faceZones().checkDefinition(true);
163 if (mesh.faceZones().checkParallelSync(true))
164 {
165 noFailedChecks++;
166 }
167 mesh.pointZones().checkDefinition(true);
168 if (mesh.pointZones().checkParallelSync(true))
169 {
170 noFailedChecks++;
171 }
172
173
174 {
175 cellSet cells(mesh, "illegalCells", mesh.nCells()/100);
176
177 forAll(mesh.cells(), celli)
178 {
179 const cell& cFaces = mesh.cells()[celli];
180
181 if (cFaces.size() <= 3)
182 {
183 cells.insert(celli);
184 }
185 for (const label facei : cFaces)
186 {
187 if (facei < 0 || facei >= mesh.nFaces())
188 {
189 cells.insert(celli);
190 break;
191 }
192 }
193 }
194 const label nCells = returnReduce(cells.size(), sumOp<label>());
195
196 if (nCells > 0)
197 {
198 Info<< " Illegal cells (less than 4 faces or out of range faces)"
199 << " found, number of cells: " << nCells << endl;
200 noFailedChecks++;
201
202 Info<< " <<Writing " << nCells
203 << " illegal cells to set " << cells.name() << endl;
204 cells.instance() = mesh.pointsInstance();
205 cells.write();
206 if (surfWriter && surfWriter->enabled())
207 {
208 mergeAndWrite(*surfWriter, cells);
209 }
210 }
211 else
212 {
213 Info<< " Cell to face addressing OK." << endl;
214 }
215 }
216
217
218 {
219 pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
220 if (mesh.checkPoints(true, &points))
221 {
222 noFailedChecks++;
223
224 label nPoints = returnReduce(points.size(), sumOp<label>());
225
226 Info<< " <<Writing " << nPoints
227 << " unused points to set " << points.name() << endl;
228 points.instance() = mesh.pointsInstance();
229 points.write();
230 if (setWriter && setWriter->enabled())
231 {
232 mergeAndWrite(*setWriter, points);
233 }
234 }
235 }
236
237 {
238 faceSet faces(mesh, "upperTriangularFace", mesh.nFaces()/100);
239 if (mesh.checkUpperTriangular(true, &faces))
240 {
241 noFailedChecks++;
242 }
243
244 const label nFaces = returnReduce(faces.size(), sumOp<label>());
245
246 if (nFaces > 0)
247 {
248 Info<< " <<Writing " << nFaces
249 << " unordered faces to set " << faces.name() << endl;
250 faces.instance() = mesh.pointsInstance();
251 faces.write();
252 if (surfWriter && surfWriter->enabled())
253 {
254 mergeAndWrite(*surfWriter, faces);
255 }
256 }
257 }
258
259 {
260 faceSet faces(mesh, "outOfRangeFaces", mesh.nFaces()/100);
261 if (mesh.checkFaceVertices(true, &faces))
262 {
263 noFailedChecks++;
264
265 const label nFaces = returnReduce(faces.size(), sumOp<label>());
266
267 Info<< " <<Writing " << nFaces
268 << " faces with out-of-range or duplicate vertices to set "
269 << faces.name() << endl;
270 faces.instance() = mesh.pointsInstance();
271 faces.write();
272 if (surfWriter && surfWriter->enabled())
273 {
274 mergeAndWrite(*surfWriter, faces);
275 }
276 }
277 }
278
279 if (allTopology)
280 {
281 cellSet cells(mesh, "zipUpCells", mesh.nCells()/100);
282 if (mesh.checkCellsZipUp(true, &cells))
283 {
284 noFailedChecks++;
285
286 const label nCells = returnReduce(cells.size(), sumOp<label>());
287
288 Info<< " <<Writing " << nCells
289 << " cells with over used edges to set " << cells.name()
290 << endl;
291 cells.instance() = mesh.pointsInstance();
292 cells.write();
293 if (surfWriter && surfWriter->enabled())
294 {
295 mergeAndWrite(*surfWriter, cells);
296 }
297 }
298 }
299
300 if (allTopology)
301 {
302 faceSet faces(mesh, "edgeFaces", mesh.nFaces()/100);
303 if (mesh.checkFaceFaces(true, &faces))
304 {
305 noFailedChecks++;
306 }
307
308 const label nFaces = returnReduce(faces.size(), sumOp<label>());
309 if (nFaces > 0)
310 {
311 Info<< " <<Writing " << nFaces
312 << " faces with non-standard edge connectivity to set "
313 << faces.name() << endl;
314 faces.instance() = mesh.pointsInstance();
315 faces.write();
316 if (surfWriter && surfWriter->enabled())
317 {
318 mergeAndWrite(*surfWriter, faces);
319 }
320 }
321 }
322
323 if (allTopology)
324 {
325 labelList nInternalFaces(mesh.nCells(), Zero);
326
327 for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
328 {
329 nInternalFaces[mesh.faceOwner()[facei]]++;
330 nInternalFaces[mesh.faceNeighbour()[facei]]++;
331 }
332 const polyBoundaryMesh& patches = mesh.boundaryMesh();
333 forAll(patches, patchi)
334 {
335 if (patches[patchi].coupled())
336 {
337 const labelUList& owners = patches[patchi].faceCells();
338
339 for (const label facei : owners)
340 {
341 nInternalFaces[facei]++;
342 }
343 }
344 }
345
346 cellSet oneCells(mesh, "oneInternalFaceCells", mesh.nCells()/100);
347 cellSet twoCells(mesh, "twoInternalFacesCells", mesh.nCells()/100);
348
349 forAll(nInternalFaces, celli)
350 {
351 if (nInternalFaces[celli] <= 1)
352 {
353 oneCells.insert(celli);
354 }
355 else if (nInternalFaces[celli] == 2)
356 {
357 twoCells.insert(celli);
358 }
359 }
360
361 label nOneCells = returnReduce(oneCells.size(), sumOp<label>());
362
363 if (nOneCells > 0)
364 {
365 Info<< " <<Writing " << nOneCells
366 << " cells with zero or one non-boundary face to set "
367 << oneCells.name()
368 << endl;
369 oneCells.instance() = mesh.pointsInstance();
370 oneCells.write();
371 if (surfWriter && surfWriter->enabled())
372 {
373 mergeAndWrite(*surfWriter, oneCells);
374 }
375 }
376
377 label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
378
379 if (nTwoCells > 0)
380 {
381 Info<< " <<Writing " << nTwoCells
382 << " cells with two non-boundary faces to set "
383 << twoCells.name()
384 << endl;
385 twoCells.instance() = mesh.pointsInstance();
386 twoCells.write();
387 if (surfWriter && surfWriter->enabled())
388 {
389 mergeAndWrite(*surfWriter, twoCells);
390 }
391 }
392 }
393
394 {
395 regionSplit rs(mesh);
396
397 if (rs.nRegions() <= 1)
398 {
399 Info<< " Number of regions: " << rs.nRegions() << " (OK)."
400 << endl;
401
402 }
403 else
404 {
405 Info<< " *Number of regions: " << rs.nRegions() << endl;
406
407 Info<< " The mesh has multiple regions which are not connected "
408 "by any face." << endl
409 << " <<Writing region information to "
410 << mesh.time().timeName()/"cellToRegion"
411 << endl;
412
413 labelIOList ctr
414 (
415 IOobject
416 (
417 "cellToRegion",
418 mesh.time().timeName(),
419 mesh,
420 IOobject::NO_READ,
421 IOobject::NO_WRITE
422 ),
423 rs
424 );
425 ctr.write();
426
427
428 // Points in multiple regions
429 pointSet points
430 (
431 mesh,
432 "multiRegionPoints",
433 mesh.nPoints()/1000
434 );
435
436 // Is region disconnected
437 boolList regionDisconnected(rs.nRegions(), true);
438 if (allTopology)
439 {
440 // -1 : not assigned
441 // -2 : multiple regions
442 // >= 0 : single region
443 labelList pointToRegion(mesh.nPoints(), -1);
444
445 for
446 (
447 label facei = mesh.nInternalFaces();
448 facei < mesh.nFaces();
449 ++facei
450 )
451 {
452 const label regioni = rs[mesh.faceOwner()[facei]];
453 const face& f = mesh.faces()[facei];
454 for (const label verti : f)
455 {
456 label& pRegion = pointToRegion[verti];
457 if (pRegion == -1)
458 {
459 pRegion = regioni;
460 }
461 else if (pRegion == -2)
462 {
463 // Already marked
464 regionDisconnected[regioni] = false;
465 }
466 else if (pRegion != regioni)
467 {
468 // Multiple regions
469 regionDisconnected[regioni] = false;
470 regionDisconnected[pRegion] = false;
471 pRegion = -2;
472 points.insert(verti);
473 }
474 }
475 }
476
477 Pstream::listCombineAllGather
478 (
479 regionDisconnected,
480 andEqOp<bool>()
481 );
482 }
483
484
485
486 // write cellSet for each region
487 PtrList<cellSet> cellRegions(rs.nRegions());
488 for (label i = 0; i < rs.nRegions(); i++)
489 {
490 cellRegions.set
491 (
492 i,
493 new cellSet
494 (
495 mesh,
496 "region" + Foam::name(i),
497 mesh.nCells()/100
498 )
499 );
500 }
501
502 forAll(rs, i)
503 {
504 cellRegions[rs[i]].insert(i);
505 }
506
507 for (label i = 0; i < rs.nRegions(); i++)
508 {
509 Info<< " <<Writing region " << i;
510 if (allTopology)
511 {
512 if (regionDisconnected[i])
513 {
514 Info<< " (fully disconnected)";
515 }
516 else
517 {
518 Info<< " (point connected)";
519 }
520 }
521 Info<< " with "
522 << returnReduce(cellRegions[i].size(), sumOp<scalar>())
523 << " cells to cellSet " << cellRegions[i].name() << endl;
524
525 cellRegions[i].write();
526 }
527
528 label nPoints = returnReduce(points.size(), sumOp<label>());
529 if (nPoints)
530 {
531 Info<< " <<Writing " << nPoints
532 << " points that are in multiple regions to set "
533 << points.name() << endl;
534 points.write();
535 if (setWriter && setWriter->enabled())
536 {
537 mergeAndWrite(*setWriter, points);
538 }
539 }
540 }
541 }
542
543 // Non-manifold points
544 pointSet points
545 (
546 mesh,
547 "nonManifoldPoints",
548 mesh.nPoints()/1000
549 );
550
551 {
552 if (!Pstream::parRun())
553 {
554 Info<< "\nChecking patch topology for multiply connected"
555 << " surfaces..." << endl;
556 }
557 else
558 {
559 Info<< "\nChecking basic patch addressing..." << endl;
560 }
561
562 const polyBoundaryMesh& patches = mesh.boundaryMesh();
563
564 Pout.setf(ios_base::left);
565
566 Info<< " "
567 << setw(20) << "Patch"
568 << setw(9) << "Faces"
569 << setw(9) << "Points";
570 if (!Pstream::parRun())
571 {
572 Info<< setw(34) << "Surface topology";
573 }
574 if (allGeometry)
575 {
576 Info<< " Bounding box";
577 }
578 Info<< endl;
579
580 forAll(patches, patchi)
581 {
582 const polyPatch& pp = patches[patchi];
583
584 if (!isA<processorPolyPatch>(pp))
585 {
586 checkPatch(allGeometry, pp.name(), pp, points);
587 Info<< endl;
588 }
589 }
590
591 //Info.setf(ios_base::right);
592 }
593
594 {
595 if (!Pstream::parRun())
596 {
597 Info<< "\nChecking faceZone topology for multiply connected"
598 << " surfaces..." << endl;
599 }
600 else
601 {
602 Info<< "\nChecking basic faceZone addressing..." << endl;
603 }
604
605 Pout.setf(ios_base::left);
606
607 const faceZoneMesh& faceZones = mesh.faceZones();
608
609 if (faceZones.size())
610 {
611 Info<< " "
612 << setw(20) << "FaceZone"
613 << setw(9) << "Faces"
614 << setw(9) << "Points";
615
616 if (!Pstream::parRun())
617 {
618 Info<< setw(34) << "Surface topology";
619 }
620 if (allGeometry)
621 {
622 Info<< " Bounding box";
623 }
624 Info<< endl;
625
626 for (const faceZone& fz : faceZones)
627 {
628 checkPatch(allGeometry, fz.name(), fz(), points);
629 Info<< endl;
630 }
631 }
632 else
633 {
634 Info<< " No faceZones found."<<endl;
635 }
636 }
637
638 const label nPoints = returnReduce(points.size(), sumOp<label>());
639
640 if (nPoints)
641 {
642 Info<< " <<Writing " << nPoints
643 << " conflicting points to set " << points.name() << endl;
644 points.instance() = mesh.pointsInstance();
645 points.write();
646 if (setWriter && setWriter->enabled())
647 {
648 mergeAndWrite(*setWriter, points);
649 }
650 }
651
652 {
653 Info<< "\nChecking basic cellZone addressing..." << endl;
654
655 Pout.setf(ios_base::left);
656
657 const cellZoneMesh& cellZones = mesh.cellZones();
658
659 if (cellZones.size())
660 {
661 Info<< " "
662 << setw(20) << "CellZone"
663 << setw(13) << "Cells"
664 << setw(13) << "Points"
665 << setw(13) << "Volume"
666 << "BoundingBox" << endl;
667
668 const cellList& cells = mesh.cells();
669 const faceList& faces = mesh.faces();
670 const scalarField& cellVolumes = mesh.cellVolumes();
671
672 bitSet isZonePoint(mesh.nPoints());
673
674 for (const cellZone& cZone : cellZones)
675 {
676 boundBox bb;
677 isZonePoint.reset(); // clears all bits (reset count)
678 scalar v = 0.0;
679
680 for (const label celli : cZone)
681 {
682 v += cellVolumes[celli];
683 for (const label facei : cells[celli])
684 {
685 const face& f = faces[facei];
686 for (const label verti : f)
687 {
688 if (isZonePoint.set(verti))
689 {
690 bb.add(mesh.points()[verti]);
691 }
692 }
693 }
694 }
695
696 bb.reduce(); // Global min/max
697
698 Info<< " "
699 << setw(19) << cZone.name()
700 << ' ' << setw(12)
701 << returnReduce(cZone.size(), sumOp<label>())
702 << ' ' << setw(12)
703 << returnReduce(isZonePoint.count(), sumOp<label>())
704 << ' ' << setw(12)
705 << returnReduce(v, sumOp<scalar>())
706 << ' ' << bb << endl;
707 }
708 }
709 else
710 {
711 Info<< " No cellZones found."<<endl;
712 }
713 }
714
715 // Force creation of all addressing if requested.
716 // Errors will be reported as required
717 if (allTopology)
718 {
719 mesh.cells();
720 mesh.faces();
721 mesh.edges();
722 mesh.points();
723 mesh.faceOwner();
724 mesh.faceNeighbour();
725 mesh.cellCells();
726 mesh.edgeCells();
727 mesh.pointCells();
728 mesh.edgeFaces();
729 mesh.pointFaces();
730 mesh.cellEdges();
731 mesh.faceEdges();
732 mesh.pointEdges();
733 }
734
735 return noFailedChecks;
736}
737
738
739// ************************************************************************* //
Istream and Ostream manipulators taking arguments.
reduce(hasMovingMesh, orOp< bool >())
ios_base::fmtflags setf(const ios_base::fmtflags f)
Set flags of stream.
Definition: IOstream.H:378
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
const polyBoundaryMesh & patches
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
dynamicFvMesh & mesh
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const pointField & points
label nPoints
const cellShapeList & cells
const dimensionedScalar mp
Proton mass.
List< label > labelList
A List of labels.
Definition: List.H:66
ZoneMesh< cellZone, polyMesh > cellZoneMesh
A ZoneMesh with the type cellZone.
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
void checkPatch(const bool allGeometry, const word &name, const PatchType &pp, pointSet &points)
label checkTopology(const polyMesh &mesh, const bool allTopology, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, autoPtr< coordSetWriter > &setWriter)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
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.
IOList< label > labelIOList
Label container classes.
Definition: labelIOList.H:44
List< bool > boolList
A List of bools.
Definition: List.H:64
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
labelList f(nPoints)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333