renumberMesh.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) 2016-2021 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
27Application
28 renumberMesh
29
30Group
31 grpMeshManipulationUtilities
32
33Description
34 Renumbers the cell list in order to reduce the bandwidth, reading and
35 renumbering all fields from all the time directories.
36
37 By default uses bandCompression (CuthillMcKee) but will
38 read system/renumberMeshDict if -dict option is present
39
40\*---------------------------------------------------------------------------*/
41
42#include "argList.H"
43#include "IOobjectList.H"
44#include "fvMesh.H"
45#include "polyTopoChange.H"
46#include "ReadFields.H"
47#include "volFields.H"
48#include "surfaceFields.H"
49#include "SortableList.H"
50#include "decompositionMethod.H"
51#include "renumberMethod.H"
54#include "fvMeshSubset.H"
55#include "cellSet.H"
56#include "faceSet.H"
57#include "pointSet.H"
58#include "processorMeshes.H"
59#include "hexRef8Data.H"
60#include "regionProperties.H"
61#include "polyMeshTools.H"
62
63#ifdef HAVE_ZOLTAN
64 #include "zoltanRenumber.H"
65#endif
66
67
68using namespace Foam;
69
70
71// Create named field from labelList for post-processing
72tmp<volScalarField> createScalarField
73(
74 const fvMesh& mesh,
75 const word& name,
76 const labelList& elems
77)
78{
80 (
82 (
84 (
85 name,
86 mesh.time().timeName(),
87 mesh,
88 IOobject::NO_READ,
89 IOobject::AUTO_WRITE,
90 false
91 ),
92 mesh,
93 dimensionedScalar(dimless, Zero),
95 )
96 );
97 volScalarField& fld = tfld.ref();
98
99 forAll(fld, celli)
100 {
101 fld[celli] = elems[celli];
102 }
103
104 return tfld;
105}
106
107
108// Calculate band of matrix
109label getBand(const labelList& owner, const labelList& neighbour)
110{
111 label band = 0;
112
113 forAll(neighbour, facei)
114 {
115 label diff = neighbour[facei] - owner[facei];
116
117 if (diff > band)
118 {
119 band = diff;
120 }
121 }
122 return band;
123}
124
125
126// Calculate band of matrix
127void getBand
128(
129 const bool calculateIntersect,
130 const label nCells,
131 const labelList& owner,
132 const labelList& neighbour,
133 label& bandwidth,
134 scalar& profile, // scalar to avoid overflow
135 scalar& sumSqrIntersect // scalar to avoid overflow
136)
137{
138 labelList cellBandwidth(nCells, Zero);
139 scalarField nIntersect(nCells, Zero);
140
141 forAll(neighbour, facei)
142 {
143 label own = owner[facei];
144 label nei = neighbour[facei];
145
146 // Note: mag not necessary for correct (upper-triangular) ordering.
147 label diff = nei-own;
148 cellBandwidth[nei] = max(cellBandwidth[nei], diff);
149 }
150
151 bandwidth = max(cellBandwidth);
152
153 // Do not use field algebra because of conversion label to scalar
154 profile = 0.0;
155 forAll(cellBandwidth, celli)
156 {
157 profile += 1.0*cellBandwidth[celli];
158 }
159
160 sumSqrIntersect = 0.0;
161 if (calculateIntersect)
162 {
163 forAll(nIntersect, celli)
164 {
165 for (label colI = celli-cellBandwidth[celli]; colI <= celli; colI++)
166 {
167 nIntersect[colI] += 1.0;
168 }
169 }
170
171 sumSqrIntersect = sum(Foam::sqr(nIntersect));
172 }
173}
174
175
176// Determine upper-triangular face order
177labelList getFaceOrder
178(
179 const primitiveMesh& mesh,
180 const labelList& cellOrder // New to old cell
181)
182{
183 labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
184
185 labelList oldToNewFace(mesh.nFaces(), -1);
186
187 label newFacei = 0;
188
189 labelList nbr;
190 labelList order;
191
192 forAll(cellOrder, newCelli)
193 {
194 label oldCelli = cellOrder[newCelli];
195
196 const cell& cFaces = mesh.cells()[oldCelli];
197
198 // Neighbouring cells
199 nbr.setSize(cFaces.size());
200
201 forAll(cFaces, i)
202 {
203 label facei = cFaces[i];
204
205 if (mesh.isInternalFace(facei))
206 {
207 // Internal face. Get cell on other side.
208 label nbrCelli = reverseCellOrder[mesh.faceNeighbour()[facei]];
209 if (nbrCelli == newCelli)
210 {
211 nbrCelli = reverseCellOrder[mesh.faceOwner()[facei]];
212 }
213
214 if (newCelli < nbrCelli)
215 {
216 // Celli is master
217 nbr[i] = nbrCelli;
218 }
219 else
220 {
221 // nbrCell is master. Let it handle this face.
222 nbr[i] = -1;
223 }
224 }
225 else
226 {
227 // External face. Do later.
228 nbr[i] = -1;
229 }
230 }
231
232 sortedOrder(nbr, order);
233
234 for (const label index : order)
235 {
236 if (nbr[index] != -1)
237 {
238 oldToNewFace[cFaces[index]] = newFacei++;
239 }
240 }
241 }
242
243 // Leave patch faces intact.
244 for (label facei = newFacei; facei < mesh.nFaces(); facei++)
245 {
246 oldToNewFace[facei] = facei;
247 }
248
249
250 // Check done all faces.
251 forAll(oldToNewFace, facei)
252 {
253 if (oldToNewFace[facei] == -1)
254 {
256 << "Did not determine new position" << " for face " << facei
257 << abort(FatalError);
258 }
259 }
260
261 return invert(mesh.nFaces(), oldToNewFace);
262}
263
264
265// Determine face order such that inside region faces are sorted
266// upper-triangular but inbetween region faces are handled like boundary faces.
267labelList getRegionFaceOrder
268(
269 const primitiveMesh& mesh,
270 const labelList& cellOrder, // New to old cell
271 const labelList& cellToRegion // Old cell to region
272)
273{
274 labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
275
276 labelList oldToNewFace(mesh.nFaces(), -1);
277
278 label newFacei = 0;
279
280 label prevRegion = -1;
281
282 forAll(cellOrder, newCelli)
283 {
284 label oldCelli = cellOrder[newCelli];
285
286 if (cellToRegion[oldCelli] != prevRegion)
287 {
288 prevRegion = cellToRegion[oldCelli];
289 }
290
291 const cell& cFaces = mesh.cells()[oldCelli];
292
293 SortableList<label> nbr(cFaces.size());
294
295 forAll(cFaces, i)
296 {
297 label facei = cFaces[i];
298
299 if (mesh.isInternalFace(facei))
300 {
301 // Internal face. Get cell on other side.
302 label nbrCelli = reverseCellOrder[mesh.faceNeighbour()[facei]];
303 if (nbrCelli == newCelli)
304 {
305 nbrCelli = reverseCellOrder[mesh.faceOwner()[facei]];
306 }
307
308 if (cellToRegion[oldCelli] != cellToRegion[cellOrder[nbrCelli]])
309 {
310 // Treat like external face. Do later.
311 nbr[i] = -1;
312 }
313 else if (newCelli < nbrCelli)
314 {
315 // Celli is master
316 nbr[i] = nbrCelli;
317 }
318 else
319 {
320 // nbrCell is master. Let it handle this face.
321 nbr[i] = -1;
322 }
323 }
324 else
325 {
326 // External face. Do later.
327 nbr[i] = -1;
328 }
329 }
330
331 nbr.sort();
332
333 forAll(nbr, i)
334 {
335 if (nbr[i] != -1)
336 {
337 oldToNewFace[cFaces[nbr.indices()[i]]] = newFacei++;
338 }
339 }
340 }
341
342 // Do region interfaces
343 label nRegions = max(cellToRegion)+1;
344 {
345 // Sort in increasing region
346 SortableList<label> sortKey(mesh.nFaces(), labelMax);
347
348 for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
349 {
350 label ownRegion = cellToRegion[mesh.faceOwner()[facei]];
351 label neiRegion = cellToRegion[mesh.faceNeighbour()[facei]];
352
353 if (ownRegion != neiRegion)
354 {
355 sortKey[facei] =
356 min(ownRegion, neiRegion)*nRegions
357 +max(ownRegion, neiRegion);
358 }
359 }
360 sortKey.sort();
361
362 // Extract.
363 label prevKey = -1;
364 forAll(sortKey, i)
365 {
366 label key = sortKey[i];
367
368 if (key == labelMax)
369 {
370 break;
371 }
372
373 if (prevKey != key)
374 {
375 prevKey = key;
376 }
377
378 oldToNewFace[sortKey.indices()[i]] = newFacei++;
379 }
380 }
381
382 // Leave patch faces intact.
383 for (label facei = newFacei; facei < mesh.nFaces(); facei++)
384 {
385 oldToNewFace[facei] = facei;
386 }
387
388
389 // Check done all faces.
390 forAll(oldToNewFace, facei)
391 {
392 if (oldToNewFace[facei] == -1)
393 {
395 << "Did not determine new position"
396 << " for face " << facei
397 << abort(FatalError);
398 }
399 }
400
401 return invert(mesh.nFaces(), oldToNewFace);
402}
403
404
405// cellOrder: old cell for every new cell
406// faceOrder: old face for every new face. Ordering of boundary faces not
407// changed.
408autoPtr<mapPolyMesh> reorderMesh
409(
410 polyMesh& mesh,
411 const labelList& cellOrder,
412 const labelList& faceOrder
413)
414{
415 labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
416 labelList reverseFaceOrder(invert(faceOrder.size(), faceOrder));
417
418 faceList newFaces(reorder(reverseFaceOrder, mesh.faces()));
419 labelList newOwner
420 (
421 renumber
422 (
423 reverseCellOrder,
424 reorder(reverseFaceOrder, mesh.faceOwner())
425 )
426 );
427 labelList newNeighbour
428 (
429 renumber
430 (
431 reverseCellOrder,
432 reorder(reverseFaceOrder, mesh.faceNeighbour())
433 )
434 );
435
436 // Check if any faces need swapping.
437 labelHashSet flipFaceFlux(newOwner.size());
438 forAll(newNeighbour, facei)
439 {
440 label own = newOwner[facei];
441 label nei = newNeighbour[facei];
442
443 if (nei < own)
444 {
445 newFaces[facei].flip();
446 std::swap(newOwner[facei], newNeighbour[facei]);
447 flipFaceFlux.insert(facei);
448 }
449 }
450
451 const polyBoundaryMesh& patches = mesh.boundaryMesh();
452 labelList patchSizes(patches.size());
453 labelList patchStarts(patches.size());
454 labelList oldPatchNMeshPoints(patches.size());
455 labelListList patchPointMap(patches.size());
456
457 forAll(patches, patchi)
458 {
459 patchSizes[patchi] = patches[patchi].size();
460 patchStarts[patchi] = patches[patchi].start();
461 oldPatchNMeshPoints[patchi] = patches[patchi].nPoints();
462 patchPointMap[patchi] = identity(patches[patchi].nPoints());
463 }
464
465 mesh.resetPrimitives
466 (
467 autoPtr<pointField>(), // <- null: leaves points untouched
468 autoPtr<faceList>::New(std::move(newFaces)),
469 autoPtr<labelList>::New(std::move(newOwner)),
470 autoPtr<labelList>::New(std::move(newNeighbour)),
471 patchSizes,
472 patchStarts,
473 true
474 );
475
476
477 // Re-do the faceZones
478 {
479 faceZoneMesh& faceZones = mesh.faceZones();
480 faceZones.clearAddressing();
481 forAll(faceZones, zoneI)
482 {
483 faceZone& fZone = faceZones[zoneI];
484 labelList newAddressing(fZone.size());
485 boolList newFlipMap(fZone.size());
486 forAll(fZone, i)
487 {
488 label oldFacei = fZone[i];
489 newAddressing[i] = reverseFaceOrder[oldFacei];
490 if (flipFaceFlux.found(newAddressing[i]))
491 {
492 newFlipMap[i] = !fZone.flipMap()[i];
493 }
494 else
495 {
496 newFlipMap[i] = fZone.flipMap()[i];
497 }
498 }
499 labelList newToOld(sortedOrder(newAddressing));
500 fZone.resetAddressing
501 (
502 labelUIndList(newAddressing, newToOld)(),
503 boolUIndList(newFlipMap, newToOld)()
504 );
505 }
506 }
507 // Re-do the cellZones
508 {
509 cellZoneMesh& cellZones = mesh.cellZones();
510 cellZones.clearAddressing();
511 forAll(cellZones, zoneI)
512 {
513 cellZones[zoneI] = labelUIndList
514 (
515 reverseCellOrder,
516 cellZones[zoneI]
517 )();
518 Foam::sort(cellZones[zoneI]);
519 }
520 }
521
522
524 (
525 mesh, // const polyMesh& mesh,
526 mesh.nPoints(), // nOldPoints,
527 mesh.nFaces(), // nOldFaces,
528 mesh.nCells(), // nOldCells,
529 identity(mesh.nPoints()), // pointMap,
530 List<objectMap>(), // pointsFromPoints,
531 faceOrder, // faceMap,
532 List<objectMap>(), // facesFromPoints,
533 List<objectMap>(), // facesFromEdges,
534 List<objectMap>(), // facesFromFaces,
535 cellOrder, // cellMap,
536 List<objectMap>(), // cellsFromPoints,
537 List<objectMap>(), // cellsFromEdges,
538 List<objectMap>(), // cellsFromFaces,
539 List<objectMap>(), // cellsFromCells,
540 identity(mesh.nPoints()), // reversePointMap,
541 reverseFaceOrder, // reverseFaceMap,
542 reverseCellOrder, // reverseCellMap,
543 flipFaceFlux, // flipFaceFlux,
544 patchPointMap, // patchPointMap,
545 labelListList(), // pointZoneMap,
546 labelListList(), // faceZonePointMap,
547 labelListList(), // faceZoneFaceMap,
548 labelListList(), // cellZoneMap,
549 pointField(), // preMotionPoints,
550 patchStarts, // oldPatchStarts,
551 oldPatchNMeshPoints, // oldPatchNMeshPoints
552 autoPtr<scalarField>() // oldCellVolumes
553 );
554}
555
556
557// Return new to old cell numbering
558labelList regionRenumber
559(
560 const renumberMethod& method,
561 const fvMesh& mesh,
562 const labelList& cellToRegion
563)
564{
565 Info<< "Determining cell order:" << endl;
566
567 labelList cellOrder(cellToRegion.size());
568
569 label nRegions = max(cellToRegion)+1;
570
571 labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
572
573 label celli = 0;
574
575 forAll(regionToCells, regioni)
576 {
577 Info<< " region " << regioni << " starts at " << celli << endl;
578
579 // Make sure no parallel comms
580 const bool oldParRun = UPstream::parRun(false);
581
582 // Per region do a reordering.
583 fvMeshSubset subsetter(mesh, regioni, cellToRegion);
584
585 const fvMesh& subMesh = subsetter.subMesh();
586
587 labelList subCellOrder = method.renumber
588 (
589 subMesh,
590 subMesh.cellCentres()
591 );
592
593 UPstream::parRun(oldParRun); // Restore parallel state
594
595
596 const labelList& cellMap = subsetter.cellMap();
597
598 forAll(subCellOrder, i)
599 {
600 cellOrder[celli++] = cellMap[subCellOrder[i]];
601 }
602 }
603 Info<< endl;
604
605 return cellOrder;
606}
607
608
609// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
610
611int main(int argc, char *argv[])
612{
613 argList::addNote
614 (
615 "Renumber mesh cells to reduce the bandwidth"
616 );
617
618 #include "addAllRegionOptions.H"
619 #include "addOverwriteOption.H"
620 #include "addTimeOptions.H"
621
622 argList::addOption("dict", "file", "Alternative renumberMeshDict");
623
624 argList::addBoolOption
625 (
626 "frontWidth",
627 "Calculate the rms of the front-width"
628 );
629
630 argList::noFunctionObjects(); // Never use function objects
631
632 #include "setRootCase.H"
633 #include "createTime.H"
634 #include "getAllRegionOptions.H"
635
636 // Force linker to include zoltan symbols. This section is only needed since
637 // Zoltan is a static library
638 #ifdef HAVE_ZOLTAN
639 Info<< "renumberMesh built with zoltan support." << nl << endl;
640 (void)zoltanRenumber::typeName;
641 #endif
642
643 // Get times list
644 instantList Times = runTime.times();
645
646 // Set startTime and endTime depending on -time and -latestTime options
647 #include "checkTimeOptions.H"
648
649 runTime.setTime(Times[startTime], startTime);
650
651 #include "createNamedMeshes.H"
652
653 const bool readDict = args.found("dict");
654 const bool doFrontWidth = args.found("frontWidth");
655 const bool overwrite = args.found("overwrite");
656
657
658 for (fvMesh& mesh : meshes)
659 {
660 // Reset time in case of multiple meshes and not overwrite
661 runTime.setTime(Times[startTime], startTime);
662
663 const word oldInstance = mesh.pointsInstance();
664
665 label band;
666 scalar profile;
667 scalar sumSqrIntersect;
668 getBand
669 (
670 doFrontWidth,
671 mesh.nCells(),
672 mesh.faceOwner(),
673 mesh.faceNeighbour(),
674 band,
675 profile,
676 sumSqrIntersect
677 );
678
679 reduce(band, maxOp<label>());
680 reduce(profile, sumOp<scalar>());
681 scalar rmsFrontwidth = Foam::sqrt
682 (
683 returnReduce
684 (
685 sumSqrIntersect,
687 )/mesh.globalData().nTotalCells()
688 );
689
690 Info<< "Mesh " << mesh.name()
691 << " size: " << mesh.globalData().nTotalCells() << nl
692 << "Before renumbering :" << nl
693 << " band : " << band << nl
694 << " profile : " << profile << nl;
695
696 if (doFrontWidth)
697 {
698 Info<< " rms frontwidth : " << rmsFrontwidth << nl;
699 }
700
701 Info<< endl;
702
703 bool sortCoupledFaceCells = false;
704 bool writeMaps = false;
705 bool orderPoints = false;
706 label blockSize = 0;
707
708 // Construct renumberMethod
709 autoPtr<IOdictionary> renumberDictPtr;
710 autoPtr<renumberMethod> renumberPtr;
711
712 if (readDict)
713 {
714 const word dictName("renumberMeshDict");
716
717 Info<< "Renumber according to " << dictIO.name() << nl << endl;
718
719 renumberDictPtr.reset(new IOdictionary(dictIO));
720 const IOdictionary& renumberDict = renumberDictPtr();
721
722 renumberPtr = renumberMethod::New(renumberDict);
723
724 sortCoupledFaceCells = renumberDict.getOrDefault
725 (
726 "sortCoupledFaceCells",
727 false
728 );
729 if (sortCoupledFaceCells)
730 {
731 Info<< "Sorting cells on coupled boundaries to be last." << nl
732 << endl;
733 }
734
735 blockSize = renumberDict.getOrDefault("blockSize", 0);
736 if (blockSize > 0)
737 {
738 Info<< "Ordering cells into regions of size " << blockSize
739 << " (using decomposition);"
740 << " ordering faces into region-internal"
741 << " and region-external."
742 << nl << endl;
743
744 if (blockSize < 0 || blockSize >= mesh.nCells())
745 {
747 << "Block size " << blockSize
748 << " should be positive integer"
749 << " and less than the number of cells in the mesh."
750 << exit(FatalError);
751 }
752 }
753
754 orderPoints = renumberDict.getOrDefault("orderPoints", false);
755 if (orderPoints)
756 {
757 Info<< "Ordering points into internal and boundary points."
758 << nl
759 << endl;
760 }
761
762 renumberDict.readEntry("writeMaps", writeMaps);
763 if (writeMaps)
764 {
765 Info<< "Writing renumber maps (new to old) to polyMesh." << nl
766 << endl;
767 }
768 }
769 else
770 {
771 Info<< "Using default renumberMethod." << nl << endl;
772 dictionary renumberDict;
773 renumberPtr.reset(new CuthillMcKeeRenumber(renumberDict));
774 }
775
776 Info<< "Selecting renumberMethod " << renumberPtr().type() << nl
777 << endl;
778
779
780
781 // Read parallel reconstruct maps
782 labelIOList cellProcAddressing
783 (
785 (
786 "cellProcAddressing",
787 mesh.facesInstance(),
788 polyMesh::meshSubDir,
789 mesh,
790 IOobject::READ_IF_PRESENT,
791 IOobject::AUTO_WRITE
792 ),
793 labelList(0)
794 );
795
796 labelIOList faceProcAddressing
797 (
799 (
800 "faceProcAddressing",
801 mesh.facesInstance(),
802 polyMesh::meshSubDir,
803 mesh,
804 IOobject::READ_IF_PRESENT,
805 IOobject::AUTO_WRITE
806 ),
807 labelList(0)
808 );
809 labelIOList pointProcAddressing
810 (
812 (
813 "pointProcAddressing",
814 mesh.pointsInstance(),
815 polyMesh::meshSubDir,
816 mesh,
817 IOobject::READ_IF_PRESENT,
818 IOobject::AUTO_WRITE
819 ),
820 labelList(0)
821 );
822 labelIOList boundaryProcAddressing
823 (
825 (
826 "boundaryProcAddressing",
827 mesh.pointsInstance(),
828 polyMesh::meshSubDir,
829 mesh,
830 IOobject::READ_IF_PRESENT,
831 IOobject::AUTO_WRITE
832 ),
833 labelList(0)
834 );
835
836
837 // Read objects in time directory
838 IOobjectList objects(mesh, runTime.timeName());
839
840
841 // Read vol fields.
842
844 ReadFields(mesh, objects, vsFlds);
845
847 ReadFields(mesh, objects, vvFlds);
848
850 ReadFields(mesh, objects, vstFlds);
851
853 ReadFields(mesh, objects, vsymtFlds);
854
856 ReadFields(mesh, objects, vtFlds);
857
858
859 // Read surface fields.
860
862 ReadFields(mesh, objects, ssFlds);
863
865 ReadFields(mesh, objects, svFlds);
866
868 ReadFields(mesh, objects, sstFlds);
869
871 ReadFields(mesh, objects, ssymtFlds);
872
874 ReadFields(mesh, objects, stFlds);
875
876
877 // Read point fields.
878
880 ReadFields(pointMesh::New(mesh), objects, psFlds);
881
883 ReadFields(pointMesh::New(mesh), objects, pvFlds);
884
886 ReadFields(pointMesh::New(mesh), objects, pstFlds);
887
889 ReadFields(pointMesh::New(mesh), objects, psymtFlds);
890
892 ReadFields(pointMesh::New(mesh), objects, ptFlds);
893
894
895 // Read sets
896 PtrList<cellSet> cellSets;
897 PtrList<faceSet> faceSets;
898 PtrList<pointSet> pointSets;
899 {
900 // Read sets
901 IOobjectList objects(mesh, mesh.facesInstance(), "polyMesh/sets");
902 ReadFields(objects, cellSets);
903 ReadFields(objects, faceSets);
904 ReadFields(objects, pointSets);
905 }
906
907
908 Info<< endl;
909
910 // From renumbering:
911 // - from new cell/face back to original cell/face
912 labelList cellOrder;
913 labelList faceOrder;
914
915 if (blockSize > 0)
916 {
917 // Renumbering in two phases. Should be done in one so mapping of
918 // fields is done correctly!
919
920 label nBlocks = mesh.nCells()/blockSize;
921 Info<< "nBlocks = " << nBlocks << endl;
922
923 // Read decompositionMethod dictionary
924 dictionary decomposeDict(renumberDictPtr().subDict("blockCoeffs"));
925 decomposeDict.set("numberOfSubdomains", nBlocks);
926
927 const bool oldParRun = UPstream::parRun(false);
928
929 autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
930 (
931 decomposeDict
932 );
933
934 labelList cellToRegion
935 (
936 decomposePtr().decompose
937 (
938 mesh,
939 mesh.cellCentres()
940 )
941 );
942
943 UPstream::parRun(oldParRun); // Restore parallel state
944
945
946 // For debugging: write out region
947 createScalarField
948 (
949 mesh,
950 "cellDist",
951 cellToRegion
952 )().write();
953
954 Info<< nl << "Written decomposition as volScalarField to "
955 << "cellDist for use in postprocessing."
956 << nl << endl;
957
958
959 cellOrder = regionRenumber(renumberPtr(), mesh, cellToRegion);
960
961 // Determine new to old face order with new cell numbering
962 faceOrder = getRegionFaceOrder
963 (
964 mesh,
965 cellOrder,
966 cellToRegion
967 );
968 }
969 else
970 {
971 // Determines sorted back to original cell ordering
972 cellOrder = renumberPtr().renumber
973 (
974 mesh,
975 mesh.cellCentres()
976 );
977
978 if (sortCoupledFaceCells)
979 {
980 // Change order so all coupled patch faceCells are at the end.
981 const polyBoundaryMesh& pbm = mesh.boundaryMesh();
982
983 // Collect all boundary cells on coupled patches
984 label nBndCells = 0;
985 forAll(pbm, patchi)
986 {
987 if (pbm[patchi].coupled())
988 {
989 nBndCells += pbm[patchi].size();
990 }
991 }
992
993 labelList reverseCellOrder = invert(mesh.nCells(), cellOrder);
994
995 labelList bndCells(nBndCells);
996 labelList bndCellMap(nBndCells);
997 nBndCells = 0;
998 forAll(pbm, patchi)
999 {
1000 if (pbm[patchi].coupled())
1001 {
1002 const labelUList& faceCells = pbm[patchi].faceCells();
1003 forAll(faceCells, i)
1004 {
1005 label celli = faceCells[i];
1006
1007 if (reverseCellOrder[celli] != -1)
1008 {
1009 bndCells[nBndCells] = celli;
1010 bndCellMap[nBndCells++] =
1011 reverseCellOrder[celli];
1012 reverseCellOrder[celli] = -1;
1013 }
1014 }
1015 }
1016 }
1017 bndCells.setSize(nBndCells);
1018 bndCellMap.setSize(nBndCells);
1019
1020 // Sort
1021 labelList order(sortedOrder(bndCellMap));
1022
1023 // Redo newReverseCellOrder
1024 labelList newReverseCellOrder(mesh.nCells(), -1);
1025
1026 label sortedI = mesh.nCells();
1027 forAllReverse(order, i)
1028 {
1029 label origCelli = bndCells[order[i]];
1030 newReverseCellOrder[origCelli] = --sortedI;
1031 }
1032
1033 Info<< "Ordered all " << nBndCells
1034 << " cells with a coupled face"
1035 << " to the end of the cell list, starting at " << sortedI
1036 << endl;
1037
1038 // Compact
1039 sortedI = 0;
1040 forAll(cellOrder, newCelli)
1041 {
1042 label origCelli = cellOrder[newCelli];
1043 if (newReverseCellOrder[origCelli] == -1)
1044 {
1045 newReverseCellOrder[origCelli] = sortedI++;
1046 }
1047 }
1048
1049 // Update sorted back to original (unsorted) map
1050 cellOrder = invert(mesh.nCells(), newReverseCellOrder);
1051 }
1052
1053
1054 // Determine new to old face order with new cell numbering
1055 faceOrder = getFaceOrder
1056 (
1057 mesh,
1058 cellOrder // New to old cell
1059 );
1060 }
1061
1062
1063 if (!overwrite)
1064 {
1065 ++runTime;
1066 }
1067
1068
1069 // Change the mesh.
1070 autoPtr<mapPolyMesh> map = reorderMesh(mesh, cellOrder, faceOrder);
1071
1072
1073 if (orderPoints)
1074 {
1075 polyTopoChange meshMod(mesh);
1076 autoPtr<mapPolyMesh> pointOrderMap = meshMod.changeMesh
1077 (
1078 mesh,
1079 false, // inflate
1080 true, // syncParallel
1081 false, // orderCells
1082 orderPoints // orderPoints
1083 );
1084
1085 // Combine point reordering into map.
1086 const_cast<labelList&>(map().pointMap()) = labelUIndList
1087 (
1088 map().pointMap(),
1089 pointOrderMap().pointMap()
1090 )();
1091
1093 (
1094 pointOrderMap().reversePointMap(),
1095 const_cast<labelList&>(map().reversePointMap())
1096 );
1097 }
1098
1099
1100 // Update fields
1101 mesh.updateMesh(map());
1102
1103 // Update proc maps
1104 if (cellProcAddressing.headerOk())
1105 {
1106 bool localOk = (cellProcAddressing.size() == mesh.nCells());
1107
1108 if (returnReduce(localOk, andOp<bool>()))
1109 {
1110 Info<< "Renumbering processor cell decomposition map "
1111 << cellProcAddressing.name() << endl;
1112
1113 cellProcAddressing = labelList
1114 (
1115 labelUIndList(cellProcAddressing, map().cellMap())
1116 );
1117 }
1118 else
1119 {
1120 Info<< "Not writing inconsistent processor cell decomposition"
1121 << " map " << cellProcAddressing.filePath() << endl;
1122 cellProcAddressing.writeOpt(IOobject::NO_WRITE);
1123 }
1124 }
1125 else
1126 {
1127 cellProcAddressing.writeOpt(IOobject::NO_WRITE);
1128 }
1129
1130 if (faceProcAddressing.headerOk())
1131 {
1132 bool localOk = (faceProcAddressing.size() == mesh.nFaces());
1133
1134 if (returnReduce(localOk, andOp<bool>()))
1135 {
1136 Info<< "Renumbering processor face decomposition map "
1137 << faceProcAddressing.name() << endl;
1138
1139 faceProcAddressing = labelList
1140 (
1141 labelUIndList(faceProcAddressing, map().faceMap())
1142 );
1143
1144 // Detect any flips.
1145 const labelHashSet& fff = map().flipFaceFlux();
1146 for (const label facei : fff)
1147 {
1148 label masterFacei = faceProcAddressing[facei];
1149
1150 faceProcAddressing[facei] = -masterFacei;
1151
1152 if (masterFacei == 0)
1153 {
1155 << " masterFacei:" << masterFacei
1156 << exit(FatalError);
1157 }
1158 }
1159 }
1160 else
1161 {
1162 Info<< "Not writing inconsistent processor face decomposition"
1163 << " map " << faceProcAddressing.filePath() << endl;
1164 faceProcAddressing.writeOpt(IOobject::NO_WRITE);
1165 }
1166 }
1167 else
1168 {
1169 faceProcAddressing.writeOpt(IOobject::NO_WRITE);
1170 }
1171
1172 if (pointProcAddressing.headerOk())
1173 {
1174 bool localOk = (pointProcAddressing.size() == mesh.nPoints());
1175
1176 if (returnReduce(localOk, andOp<bool>()))
1177 {
1178 Info<< "Renumbering processor point decomposition map "
1179 << pointProcAddressing.name() << endl;
1180
1181 pointProcAddressing = labelList
1182 (
1183 labelUIndList(pointProcAddressing, map().pointMap())
1184 );
1185 }
1186 else
1187 {
1188 Info<< "Not writing inconsistent processor point decomposition"
1189 << " map " << pointProcAddressing.filePath() << endl;
1190 pointProcAddressing.writeOpt(IOobject::NO_WRITE);
1191 }
1192 }
1193 else
1194 {
1195 pointProcAddressing.writeOpt(IOobject::NO_WRITE);
1196 }
1197
1198 if (boundaryProcAddressing.headerOk())
1199 {
1200 bool localOk =
1201 (
1202 boundaryProcAddressing.size()
1203 == mesh.boundaryMesh().size()
1204 );
1205 if (returnReduce(localOk, andOp<bool>()))
1206 {
1207 // No renumbering needed
1208 }
1209 else
1210 {
1211 Info<< "Not writing inconsistent processor patch decomposition"
1212 << " map " << boundaryProcAddressing.filePath() << endl;
1213 boundaryProcAddressing.writeOpt(IOobject::NO_WRITE);
1214 }
1215 }
1216 else
1217 {
1218 boundaryProcAddressing.writeOpt(IOobject::NO_WRITE);
1219 }
1220
1221
1222
1223
1224 // Move mesh (since morphing might not do this)
1225 if (map().hasMotionPoints())
1226 {
1227 mesh.movePoints(map().preMotionPoints());
1228 }
1229
1230
1231 {
1232 label band;
1233 scalar profile;
1234 scalar sumSqrIntersect;
1235 getBand
1236 (
1237 doFrontWidth,
1238 mesh.nCells(),
1239 mesh.faceOwner(),
1240 mesh.faceNeighbour(),
1241 band,
1242 profile,
1243 sumSqrIntersect
1244 );
1245 reduce(band, maxOp<label>());
1246 reduce(profile, sumOp<scalar>());
1247 scalar rmsFrontwidth = Foam::sqrt
1248 (
1249 returnReduce
1250 (
1251 sumSqrIntersect,
1253 )/mesh.globalData().nTotalCells()
1254 );
1255
1256 Info<< "After renumbering :" << nl
1257 << " band : " << band << nl
1258 << " profile : " << profile << nl;
1259
1260 if (doFrontWidth)
1261 {
1262
1263 Info<< " rms frontwidth : " << rmsFrontwidth << nl;
1264 }
1265
1266 Info<< endl;
1267 }
1268
1269 if (orderPoints)
1270 {
1271 // Force edge calculation (since only reason that points would
1272 // need to be sorted)
1273 (void)mesh.edges();
1274
1275 label nTotPoints = returnReduce
1276 (
1277 mesh.nPoints(),
1278 sumOp<label>()
1279 );
1280 label nTotIntPoints = returnReduce
1281 (
1282 mesh.nInternalPoints(),
1283 sumOp<label>()
1284 );
1285
1286 label nTotEdges = returnReduce
1287 (
1288 mesh.nEdges(),
1289 sumOp<label>()
1290 );
1291 label nTotIntEdges = returnReduce
1292 (
1293 mesh.nInternalEdges(),
1294 sumOp<label>()
1295 );
1296 label nTotInt0Edges = returnReduce
1297 (
1298 mesh.nInternal0Edges(),
1299 sumOp<label>()
1300 );
1301 label nTotInt1Edges = returnReduce
1302 (
1303 mesh.nInternal1Edges(),
1304 sumOp<label>()
1305 );
1306
1307 Info<< "Points:" << nl
1308 << " total : " << nTotPoints << nl
1309 << " internal: " << nTotIntPoints << nl
1310 << " boundary: " << nTotPoints-nTotIntPoints << nl
1311 << "Edges:" << nl
1312 << " total : " << nTotEdges << nl
1313 << " internal: " << nTotIntEdges << nl
1314 << " internal using 0 boundary points: "
1315 << nTotInt0Edges << nl
1316 << " internal using 1 boundary points: "
1317 << nTotInt1Edges-nTotInt0Edges << nl
1318 << " internal using 2 boundary points: "
1319 << nTotIntEdges-nTotInt1Edges << nl
1320 << " boundary: " << nTotEdges-nTotIntEdges << nl
1321 << endl;
1322 }
1323
1324
1325 if (overwrite)
1326 {
1327 mesh.setInstance(oldInstance);
1328 }
1329 else
1330 {
1331 mesh.setInstance(runTime.timeName());
1332 }
1333
1334
1335 Info<< "Writing mesh to " << mesh.facesInstance() << endl;
1336
1337 // Remove old procAddressing files
1338 processorMeshes::removeFiles(mesh);
1339
1340 // Update refinement data
1341 hexRef8Data refData
1342 (
1343 IOobject
1344 (
1345 "dummy",
1346 mesh.facesInstance(),
1347 polyMesh::meshSubDir,
1348 mesh,
1349 IOobject::READ_IF_PRESENT,
1350 IOobject::NO_WRITE,
1351 false
1352 )
1353 );
1354 refData.updateMesh(map());
1355 refData.write();
1356
1357 // Update sets
1358 topoSet::updateMesh(mesh.facesInstance(), map(), cellSets);
1359 topoSet::updateMesh(mesh.facesInstance(), map(), faceSets);
1360 topoSet::updateMesh(mesh.facesInstance(), map(), pointSets);
1361
1362 mesh.write();
1363
1364 if (writeMaps)
1365 {
1366 // For debugging: write out region
1367 createScalarField
1368 (
1369 mesh,
1370 "origCellID",
1371 map().cellMap()
1372 )().write();
1373
1374 createScalarField
1375 (
1376 mesh,
1377 "cellID",
1378 identity(mesh.nCells())
1379 )().write();
1380
1381 Info<< nl
1382 << "Written current cellID and origCellID as volScalarField"
1383 << " for use in postprocessing." << nl << endl;
1384
1386 (
1387 IOobject
1388 (
1389 "cellMap",
1390 mesh.facesInstance(),
1391 polyMesh::meshSubDir,
1392 mesh,
1393 IOobject::NO_READ,
1394 IOobject::NO_WRITE,
1395 false
1396 ),
1397 map().cellMap()
1398 ).write();
1399
1401 (
1402 IOobject
1403 (
1404 "faceMap",
1405 mesh.facesInstance(),
1406 polyMesh::meshSubDir,
1407 mesh,
1408 IOobject::NO_READ,
1409 IOobject::NO_WRITE,
1410 false
1411 ),
1412 map().faceMap()
1413 ).write();
1414
1416 (
1417 IOobject
1418 (
1419 "pointMap",
1420 mesh.facesInstance(),
1421 polyMesh::meshSubDir,
1422 mesh,
1423 IOobject::NO_READ,
1424 IOobject::NO_WRITE,
1425 false
1426 ),
1427 map().pointMap()
1428 ).write();
1429 }
1430 }
1431
1432 Info<< "End\n" << endl;
1433
1434 return 0;
1435}
1436
1437
1438// ************************************************************************* //
Field reading functions for post-processing utilities.
Y[inertIndex] max(0.0)
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
reduce(hasMovingMesh, orOp< bool >())
Foam::label startTime
Cuthill-McKee renumbering.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:59
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableList.H:63
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
void clearAddressing()
Clear addressing.
Definition: ZoneMesh.C:715
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
void reset(autoPtr< T > &&other) noexcept
Delete managed object and set to new given pointer.
Definition: autoPtrI.H:117
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:59
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
virtual void resetAddressing(const labelUList &addr, const bool flipMapValue)
Reset addressing - use uniform flip map value.
Definition: faceZone.C:444
const boolList & flipMap() const noexcept
Return face flip map.
Definition: faceZone.H:272
Holds a reference to the original mesh (the baseMesh) and optionally to a subset of that mesh (the su...
Definition: fvMeshSubset.H:80
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Various for reading/decomposing/reconstructing/distributing refinement data.
Definition: hexRef8Data.H:61
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
UPtrList< const labelUList > faceCells() const
Return a list of faceCells for each patch.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Direct mesh changes based on v1.3 polyTopoChange syntax.
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:79
const vectorField & cellCentres() const
virtual bool write(const bool valid=true) const
Write using setting from DB.
Abstract base class for renumbering.
virtual labelList renumber(const pointField &) const
A class for managing temporary objects.
Definition: tmp.H:65
A class for handling words, derived from Foam::string.
Definition: word.H:68
const polyBoundaryMesh & patches
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
Foam::PtrList< Foam::fvMesh > meshes(regionNames.size())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const word dictName("faMeshDefinition")
label nPoints
void inplaceRenumber(const labelUList &oldToNew, IntListType &lists)
Inplace renumber the values (not the indices) of a list of lists.
Definition: ensightOutput.H:90
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
Namespace for OpenFOAM.
wordList ReadFields(const typename GeoMesh::Mesh &mesh, const IOobjectList &objects, PtrList< GeometricField< Type, PatchField, GeoMesh > > &fields, const bool syncPar=true, const bool readOldTime=false)
Read Geometric fields of templated type.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:378
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
void sort(UList< T > &list)
Sort the list.
Definition: UList.C:342
errorManip< error > abort(error &err)
Definition: errorManip.H:144
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
ListType reorder(const labelUList &oldToNew, const ListType &input, const bool prune=false)
Reorder the elements of a list.
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
runTime write()
IOobject dictIO
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:346
static const char *const typeName
The type name used in ensight case files.
Foam::surfaceFields.