reconstructParMesh.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-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
27Application
28 reconstructParMesh
29
30Group
31 grpParallelUtilities
32
33Description
34 Reconstructs a mesh using geometric information only.
35
36 Writes point/face/cell procAddressing so afterwards reconstructPar can be
37 used to reconstruct fields.
38
39Usage
40 \b reconstructParMesh [OPTION]
41
42 Options:
43 - \par -fullMatch
44 Does geometric matching on all boundary faces. Assumes no point
45 ordering
46
47 - \par -procMatch
48 Assumes processor patches already in face order but not point order.
49 This is the pre v2106 default behaviour but might be removed if the new
50 topological method works well
51
52 - \par -mergeTol <tol>
53 Specifies non-default merge tolerance (fraction of mesh bounding box)
54 for above options
55
56 The default is to assume all processor boundaries are correctly ordered
57 (both faces and points) in which case no merge tolerance is needed.
58
59\*---------------------------------------------------------------------------*/
60
61#include "argList.H"
62#include "timeSelector.H"
63
64#include "IOobjectList.H"
65#include "labelIOList.H"
66#include "processorPolyPatch.H"
67#include "mapAddedPolyMesh.H"
68#include "polyMeshAdder.H"
69#include "faceCoupleInfo.H"
70#include "fvMeshAdder.H"
71#include "polyTopoChange.H"
73#include "topoSet.H"
74#include "regionProperties.H"
75#include "fvMeshTools.H"
76
77using namespace Foam;
78
79// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
80
81// Tolerance (as fraction of the bounding box). Needs to be fairly lax since
82// usually meshes get written with limited precision (6 digits)
83static const scalar defaultMergeTol = 1e-7;
84
85
86// Determine which faces are coupled. Uses geometric merge distance.
87// Looks either at all boundaryFaces (fullMatch) or only at the
88// procBoundaries for proci. Assumes that masterMesh contains already merged
89// all the processors < proci.
90autoPtr<faceCoupleInfo> determineCoupledFaces
91(
92 const bool fullMatch,
93 const label masterMeshProcStart,
94 const label masterMeshProcEnd,
95 const polyMesh& masterMesh,
96 const label meshToAddProcStart,
97 const label meshToAddProcEnd,
98 const polyMesh& meshToAdd,
99 const scalar mergeDist
100)
101{
102 if (fullMatch || masterMesh.nCells() == 0)
103 {
105 (
106 masterMesh,
107 meshToAdd,
108 mergeDist, // Absolute merging distance
109 true // Matching faces identical
110 );
111 }
112 else
113 {
114 // Pick up all patches on masterMesh ending in "toDDD" where DDD is
115 // the processor number proci.
116
117 const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
118
119
120 DynamicList<label> masterFaces
121 (
122 masterMesh.nFaces()
123 - masterMesh.nInternalFaces()
124 );
125
126
127 forAll(masterPatches, patchi)
128 {
129 const polyPatch& pp = masterPatches[patchi];
130
131 if (isA<processorPolyPatch>(pp))
132 {
133 for
134 (
135 label proci=meshToAddProcStart;
136 proci<meshToAddProcEnd;
137 proci++
138 )
139 {
140 const string toProcString("to" + name(proci));
141 if (
142 pp.name().rfind(toProcString)
143 == (pp.name().size()-toProcString.size())
144 )
145 {
146 label meshFacei = pp.start();
147 forAll(pp, i)
148 {
149 masterFaces.append(meshFacei++);
150 }
151 break;
152 }
153 }
154
155 }
156 }
157 masterFaces.shrink();
158
159
160 // Pick up all patches on meshToAdd ending in "procBoundaryDDDtoYYY"
161 // where DDD is the processor number proci and YYY is < proci.
162
163 const polyBoundaryMesh& addPatches = meshToAdd.boundaryMesh();
164
165 DynamicList<label> addFaces
166 (
167 meshToAdd.nFaces()
168 - meshToAdd.nInternalFaces()
169 );
170
171 forAll(addPatches, patchi)
172 {
173 const polyPatch& pp = addPatches[patchi];
174
175 if (isA<processorPolyPatch>(pp))
176 {
177 bool isConnected = false;
178
179 for
180 (
181 label mergedProci=masterMeshProcStart;
182 !isConnected && (mergedProci < masterMeshProcEnd);
183 mergedProci++
184 )
185 {
186 for
187 (
188 label proci = meshToAddProcStart;
189 proci < meshToAddProcEnd;
190 proci++
191 )
192 {
193 const word fromProcString
194 (
195 processorPolyPatch::newName(proci, mergedProci)
196 );
197
198 if (pp.name() == fromProcString)
199 {
200 isConnected = true;
201 break;
202 }
203 }
204 }
205
206 if (isConnected)
207 {
208 label meshFacei = pp.start();
209 forAll(pp, i)
210 {
211 addFaces.append(meshFacei++);
212 }
213 }
214 }
215 }
216 addFaces.shrink();
217
219 (
220 masterMesh,
221 masterFaces,
222 meshToAdd,
223 addFaces,
224 mergeDist, // Absolute merging distance
225 true, // Matching faces identical?
226 false, // If perfect match are faces already ordered
227 // (e.g. processor patches)
228 false // are faces each on separate patch?
229 );
230 }
231}
232
233
234autoPtr<mapPolyMesh> mergeSharedPoints
235(
236 const scalar mergeDist,
237 polyMesh& mesh,
238 labelListList& pointProcAddressing
239)
240{
241 // Find out which sets of points get merged and create a map from
242 // mesh point to unique point.
243 Map<label> pointToMaster
244 (
245 fvMeshAdder::findSharedPoints
246 (
247 mesh,
248 mergeDist
249 )
250 );
251
252 Info<< "mergeSharedPoints : detected " << pointToMaster.size()
253 << " points that are to be merged." << endl;
254
255 if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
256 {
257 return nullptr;
258 }
259
260 polyTopoChange meshMod(mesh);
261
262 fvMeshAdder::mergePoints(mesh, pointToMaster, meshMod);
263
264 // Change the mesh (no inflation). Note: parallel comms allowed.
265 autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
266
267 // Update fields. No inflation, parallel sync.
268 mesh.updateMesh(map());
269
270 // pointProcAddressing give indices into the master mesh so adapt them
271 // for changed point numbering.
272
273 // Adapt constructMaps for merged points.
274 forAll(pointProcAddressing, proci)
275 {
276 labelList& constructMap = pointProcAddressing[proci];
277
278 forAll(constructMap, i)
279 {
280 label oldPointi = constructMap[i];
281
282 // New label of point after changeMesh.
283 label newPointi = map().reversePointMap()[oldPointi];
284
285 if (newPointi < -1)
286 {
287 constructMap[i] = -newPointi-2;
288 }
289 else if (newPointi >= 0)
290 {
291 constructMap[i] = newPointi;
292 }
293 else
294 {
296 << "Problem. oldPointi:" << oldPointi
297 << " newPointi:" << newPointi << abort(FatalError);
298 }
299 }
300 }
301
302 return map;
303}
304
305
306boundBox procBounds
307(
308 const PtrList<Time>& databases,
309 const word& regionDir
310)
311{
312 boundBox bb;
313
314 for (const Time& procDb : databases)
315 {
316 fileName pointsInstance
317 (
318 procDb.findInstance
319 (
320 regionDir/polyMesh::meshSubDir,
321 "points"
322 )
323 );
324
325 if (pointsInstance != procDb.timeName())
326 {
328 << "Your time was specified as " << procDb.timeName()
329 << " but there is no polyMesh/points in that time." << nl
330 << "(points file at " << pointsInstance << ')' << nl
331 << "Please rerun with the correct time specified"
332 << " (through the -constant, -time or -latestTime "
333 << "(at your option)."
334 << endl << exit(FatalError);
335 }
336
337 Info<< "Reading points from "
338 << procDb.caseName()
339 << " for time = " << procDb.timeName()
340 << nl << endl;
341
343 (
345 (
346 "points",
347 procDb.findInstance
348 (
349 regionDir/polyMesh::meshSubDir,
350 "points"
351 ),
352 regionDir/polyMesh::meshSubDir,
353 procDb,
354 IOobject::MUST_READ,
355 IOobject::NO_WRITE,
356 false
357 )
358 );
359
360 bb.add(points);
361 }
362
363 return bb;
364}
365
366
367void writeDistribution
368(
369 Time& runTime,
370 const fvMesh& masterMesh,
371 const labelListList& cellProcAddressing
372)
373{
374 // Write the decomposition as labelList for use with 'manual'
375 // decomposition method.
376 labelIOList cellDecomposition
377 (
379 (
380 "cellDecomposition",
381 masterMesh.facesInstance(),
382 masterMesh,
383 IOobject::NO_READ,
384 IOobject::NO_WRITE,
385 false
386 ),
387 masterMesh.nCells()
388 );
389
390 forAll(cellProcAddressing, proci)
391 {
392 const labelList& pCells = cellProcAddressing[proci];
393 labelUIndList(cellDecomposition, pCells) = proci;
394 }
395
396 cellDecomposition.write();
397
398 Info<< nl << "Wrote decomposition to "
399 << cellDecomposition.objectRelPath()
400 << " for use in manual decomposition." << endl;
401
402 // Write as volScalarField for postprocessing.
403 // Change time to 0 if was 'constant'
404 {
405 const scalar oldTime = runTime.value();
406 const label oldIndex = runTime.timeIndex();
407 if (runTime.timeName() == runTime.constant() && oldIndex == 0)
408 {
409 runTime.setTime(0, oldIndex+1);
410 }
411
412 volScalarField cellDist
413 (
415 (
416 "cellDist",
417 runTime.timeName(),
418 masterMesh,
419 IOobject::NO_READ,
420 IOobject::NO_WRITE,
421 false
422 ),
423 masterMesh,
424 dimensionedScalar("cellDist", dimless, -1),
426 );
427
428 forAll(cellDecomposition, celli)
429 {
430 cellDist[celli] = cellDecomposition[celli];
431 }
432
433 cellDist.correctBoundaryConditions();
434 cellDist.write();
435
436 Info<< nl << "Wrote decomposition to "
437 << cellDist.objectRelPath()
438 << " (volScalarField) for visualization."
439 << endl;
440
441 // Restore time
442 runTime.setTime(oldTime, oldIndex);
443 }
444}
445
446
447void writeMesh
448(
449 const fvMesh& mesh,
450 const labelListList& cellProcAddressing
451)
452{
453 const fileName outputDir
454 (
455 mesh.time().path()
456 / mesh.time().timeName()
457 / mesh.regionName()
458 / polyMesh::meshSubDir
459 );
460
461 Info<< nl << "Writing merged mesh to "
462 << mesh.time().relativePath(outputDir) << nl << endl;
463
464 if (!mesh.write())
465 {
467 << "Failed writing polyMesh."
468 << exit(FatalError);
469 }
470 topoSet::removeFiles(mesh);
471}
472
473
474void writeMaps
475(
476 const label masterInternalFaces,
477 const labelUList& masterOwner,
478 const polyMesh& procMesh,
479 const labelUList& cellProcAddressing,
480 const labelUList& faceProcAddressing,
481 const labelUList& pointProcAddressing,
482 const labelUList& boundProcAddressing
483)
484{
485 const fileName outputDir
486 (
487 procMesh.time().caseName()
488 / procMesh.facesInstance()
489 / procMesh.regionName()
490 / polyMesh::meshSubDir
491 );
492
493 IOobject ioAddr
494 (
495 "procAddressing",
496 procMesh.facesInstance(),
497 polyMesh::meshSubDir,
498 procMesh,
499 IOobject::NO_READ,
500 IOobject::NO_WRITE,
501 false // Do not register
502 );
503
504
505 Info<< "Writing addressing : " << outputDir << nl;
506
507 // From processor point to reconstructed mesh point
508
509 Info<< " pointProcAddressing" << endl;
510 ioAddr.rename("pointProcAddressing");
511 IOListRef<label>(ioAddr, pointProcAddressing).write();
512
513 // From processor face to reconstructed mesh face
514 Info<< " faceProcAddressing" << endl;
515 ioAddr.rename("faceProcAddressing");
516 labelIOList faceProcAddr(ioAddr, faceProcAddressing);
517
518 // Now add turning index to faceProcAddressing.
519 // See reconstructPar for meaning of turning index.
520 forAll(faceProcAddr, procFacei)
521 {
522 const label masterFacei = faceProcAddr[procFacei];
523
524 if
525 (
526 !procMesh.isInternalFace(procFacei)
527 && masterFacei < masterInternalFaces
528 )
529 {
530 // proc face is now external but used to be internal face.
531 // Check if we have owner or neighbour.
532
533 label procOwn = procMesh.faceOwner()[procFacei];
534 label masterOwn = masterOwner[masterFacei];
535
536 if (cellProcAddressing[procOwn] == masterOwn)
537 {
538 // No turning. Offset by 1.
539 faceProcAddr[procFacei]++;
540 }
541 else
542 {
543 // Turned face.
544 faceProcAddr[procFacei] = -1 - faceProcAddr[procFacei];
545 }
546 }
547 else
548 {
549 // No turning. Offset by 1.
550 faceProcAddr[procFacei]++;
551 }
552 }
553
554 faceProcAddr.write();
555
556
557 // From processor cell to reconstructed mesh cell
558 Info<< " cellProcAddressing" << endl;
559 ioAddr.rename("cellProcAddressing");
560 IOListRef<label>(ioAddr, cellProcAddressing).write();
561
562
563 // From processor patch to reconstructed mesh patch
564 Info<< " boundaryProcAddressing" << endl;
565 ioAddr.rename("boundaryProcAddressing");
566 IOListRef<label>(ioAddr, boundProcAddressing).write();
567
568 Info<< endl;
569}
570
571
572void printWarning()
573{
574 Info<<
575"Merge individual processor meshes back into one master mesh.\n"
576"Use if the original master mesh has been deleted or the processor meshes\n"
577"have been modified (topology change).\n"
578"This tool will write the resulting mesh to a new time step and construct\n"
579"xxxxProcAddressing files in the processor meshes so reconstructPar can be\n"
580"used to regenerate the fields on the master mesh.\n\n"
581"Not well tested & use at your own risk!\n\n";
582}
583
584
585int main(int argc, char *argv[])
586{
587 argList::addNote
588 (
589 "Reconstruct a mesh using geometric/topological information only"
590 );
591
592 // Enable -constant ... if someone really wants it
593 // Enable -withZero to prevent accidentally trashing the initial fields
594 timeSelector::addOptions(true, true); // constant(true), zero(true)
595
596 argList::noParallel();
597
598 argList::addVerboseOption("Additional verbosity");
599 argList::addBoolOption
600 (
601 "addressing-only",
602 "Create procAddressing only without overwriting the mesh"
603 );
604 argList::addOption
605 (
606 "mergeTol",
607 "scalar",
608 "The merge distance relative to the bounding box size (default 1e-7)"
609 );
610 argList::addBoolOption
611 (
612 "fullMatch",
613 "Do (slower) geometric matching on all boundary faces"
614 );
615 argList::addBoolOption
616 (
617 "procMatch",
618 "Do matching on processor faces only"
619 );
620 argList::addBoolOption
621 (
622 "cellDist",
623 "Write cell distribution as a labelList - for use with 'manual' "
624 "decomposition method or as a volScalarField for post-processing."
625 );
626
627 #include "addAllRegionOptions.H"
628
629 #include "setRootCase.H"
630 #include "createTime.H"
631
632 printWarning();
633
634 const bool fullMatch = args.found("fullMatch");
635 const bool procMatch = args.found("procMatch");
636 const bool writeCellDist = args.found("cellDist");
637
638 const bool writeAddrOnly = args.found("addressing-only");
639
640 const scalar mergeTol =
641 args.getOrDefault<scalar>("mergeTol", defaultMergeTol);
642
643 if (fullMatch)
644 {
645 Info<< "Use geometric matching on all boundary faces." << nl << endl;
646 }
647 else if (procMatch)
648 {
649 Info<< "Use geometric matching on correct procBoundaries only." << nl
650 << "This assumes a correct decomposition." << endl;
651 }
652 else
653 {
654 Info<< "Merge assuming correct, fully matched procBoundaries." << nl
655 << endl;
656 }
657
658 if (fullMatch || procMatch)
659 {
660 const scalar writeTol =
661 Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
662
663 Info<< "Merge tolerance : " << mergeTol << nl
664 << "Write tolerance : " << writeTol << endl;
665
666 if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
667 {
669 << "Your current settings specify ASCII writing with "
670 << IOstream::defaultPrecision() << " digits precision." << endl
671 << "Your merging tolerance (" << mergeTol << ")"
672 << " is finer than this." << endl
673 << "Please change your writeFormat to binary"
674 << " or increase the writePrecision" << endl
675 << "or adjust the merge tolerance (-mergeTol)."
676 << exit(FatalError);
677 }
678 }
679
680 // Get region names
681 #include "getAllRegionOptions.H"
682
683 // Determine the processor count
684 label nProcs{0};
685
686 if (regionNames.empty())
687 {
689 << "No regions specified or detected."
690 << exit(FatalError);
691 }
692 else if (regionNames[0] == polyMesh::defaultRegion)
693 {
694 nProcs = fileHandler().nProcs(args.path());
695 }
696 else
697 {
698 nProcs = fileHandler().nProcs(args.path(), regionNames[0]);
699
700 if (regionNames.size() == 1)
701 {
702 Info<< "Using region: " << regionNames[0] << nl << endl;
703 }
704 }
705
706 if (!nProcs)
707 {
709 << "No processor* directories found"
710 << exit(FatalError);
711 }
712
713 // Read all time databases
714 PtrList<Time> databases(nProcs);
715
716 Info<< "Found " << nProcs << " processor directories" << endl;
717 forAll(databases, proci)
718 {
719 Info<< " Reading database "
720 << args.caseName()/("processor" + Foam::name(proci))
721 << endl;
722
723 databases.set
724 (
725 proci,
726 new Time
727 (
728 Time::controlDictName,
729 args.rootPath(),
730 args.caseName()/("processor" + Foam::name(proci))
731 )
732 );
733 }
734 Info<< endl;
735
736 // Use the times list from the master processor
737 // and select a subset based on the command-line options
738 instantList timeDirs = timeSelector::select
739 (
740 databases[0].times(),
741 args
742 );
743
744 // Loop over all times
745 forAll(timeDirs, timei)
746 {
747 // Set time for global database
748 runTime.setTime(timeDirs[timei], timei);
749
750 // Set time for all databases
751 forAll(databases, proci)
752 {
753 databases[proci].setTime(timeDirs[timei], timei);
754 }
755
756 Info<< "Time = " << runTime.timeName() << endl;
757
758 // Check for any mesh changes
759 label nMeshChanged = 0;
760 boolList hasRegionMesh(regionNames.size(), false);
761 forAll(regionNames, regioni)
762 {
763 const word& regionName = regionNames[regioni];
764
765 IOobject facesIO
766 (
767 "faces",
768 databases[0].timeName(),
769 polyMesh::regionName(regionName)/polyMesh::meshSubDir,
770 databases[0],
771 IOobject::NO_READ,
772 IOobject::NO_WRITE
773 );
774
775 // Problem: faceCompactIOList recognises both 'faceList' and
776 // 'faceCompactList' so we should be lenient when doing
777 // typeHeaderOk
778
779 const bool ok = facesIO.typeHeaderOk<faceCompactIOList>(false);
780
781 if (ok)
782 {
783 hasRegionMesh[regioni] = true;
784 ++nMeshChanged;
785 }
786 }
787
788 // Check for any mesh changes
789 if (!nMeshChanged)
790 {
791 Info<< "No meshes" << nl << endl;
792 continue;
793 }
794
795 Info<< endl;
796
797 forAll(regionNames, regioni)
798 {
799 const word& regionName = regionNames[regioni];
800 const word& regionDir = polyMesh::regionName(regionName);
801
802 if (!hasRegionMesh[regioni])
803 {
804 Info<< "region=" << regionName << " (no mesh)" << nl << endl;
805 continue;
806 }
807
808 if (regionNames.size() > 1)
809 {
810 Info<< "region=" << regionName << nl;
811 }
812
813
814 // Addressing from processor to reconstructed case
815 labelListList cellProcAddressing(nProcs);
816 labelListList faceProcAddressing(nProcs);
817 labelListList pointProcAddressing(nProcs);
818 labelListList boundProcAddressing(nProcs);
819
820
821 // Internal faces on the final reconstructed mesh
822 label masterInternalFaces;
823
824 // Owner addressing on the final reconstructed mesh
825 labelList masterOwner;
826
827 if (procMatch)
828 {
829 // Read point on individual processors to determine
830 // merge tolerance
831 // (otherwise single cell domains might give problems)
832
833 const boundBox bb = procBounds(databases, regionDir);
834 const scalar mergeDist = mergeTol*bb.mag();
835
836 Info<< "Overall mesh bounding box : " << bb << nl
837 << "Relative tolerance : " << mergeTol << nl
838 << "Absolute matching distance : " << mergeDist << nl
839 << endl;
840
841
842 // Construct empty mesh.
843 PtrList<fvMesh> masterMesh(nProcs);
844
845 for (label proci=0; proci<nProcs; proci++)
846 {
847 masterMesh.set
848 (
849 proci,
850 new fvMesh
851 (
853 (
855 runTime.timeName(),
856 runTime,
857 IOobject::NO_READ
858 ),
859 Zero
860 )
861 );
862
863 fvMesh meshToAdd
864 (
866 (
868 databases[proci].timeName(),
869 databases[proci]
870 )
871 );
872
873 // Initialize its addressing
874 cellProcAddressing[proci] = identity(meshToAdd.nCells());
875 faceProcAddressing[proci] = identity(meshToAdd.nFaces());
876 pointProcAddressing[proci] = identity(meshToAdd.nPoints());
877 boundProcAddressing[proci] =
878 identity(meshToAdd.boundaryMesh().size());
879
880 // Find geometrically shared points/faces.
881 autoPtr<faceCoupleInfo> couples = determineCoupledFaces
882 (
883 fullMatch,
884 proci,
885 proci,
886 masterMesh[proci],
887 proci,
888 proci,
889 meshToAdd,
890 mergeDist
891 );
892
893 // Add elements to mesh
895 (
896 masterMesh[proci],
897 meshToAdd,
898 couples()
899 );
900
901 // Added processor
902 renumber(map().addedCellMap(), cellProcAddressing[proci]);
903 renumber(map().addedFaceMap(), faceProcAddressing[proci]);
904 renumber(map().addedPointMap(), pointProcAddressing[proci]);
905 renumber(map().addedPatchMap(), boundProcAddressing[proci]);
906 }
907 for (label step=2; step<nProcs*2; step*=2)
908 {
909 for (label proci=0; proci<nProcs; proci+=step)
910 {
911 label next = proci + step/2;
912 if(next >= nProcs)
913 {
914 continue;
915 }
916
917 Info<< "Merging mesh " << proci << " with "
918 << next << endl;
919
920 // Find geometrically shared points/faces.
921 autoPtr<faceCoupleInfo> couples = determineCoupledFaces
922 (
923 fullMatch,
924 proci,
925 next,
926 masterMesh[proci],
927 next,
928 proci+step,
929 masterMesh[next],
930 mergeDist
931 );
932
933 // Add elements to mesh
935 (
936 masterMesh[proci],
937 masterMesh[next],
938 couples()
939 );
940
941 // Processors that were already in masterMesh
942 for (label mergedI=proci; mergedI<next; mergedI++)
943 {
945 (
946 map().oldCellMap(),
947 cellProcAddressing[mergedI]
948 );
949
951 (
952 map().oldFaceMap(),
953 faceProcAddressing[mergedI]
954 );
955
957 (
958 map().oldPointMap(),
959 pointProcAddressing[mergedI]
960 );
961
962 // Note: boundary is special since can contain -1.
964 (
965 map().oldPatchMap(),
966 boundProcAddressing[mergedI]
967 );
968 }
969
970 // Added processor
971 for
972 (
973 label addedI=next;
974 addedI<min(proci+step, nProcs);
975 addedI++
976 )
977 {
979 (
980 map().addedCellMap(),
981 cellProcAddressing[addedI]
982 );
983
985 (
986 map().addedFaceMap(),
987 faceProcAddressing[addedI]
988 );
989
991 (
992 map().addedPointMap(),
993 pointProcAddressing[addedI]
994 );
995
997 (
998 map().addedPatchMap(),
999 boundProcAddressing[addedI]
1000 );
1001 }
1002
1003 masterMesh.set(next, nullptr);
1004 }
1005 }
1006
1007 for (label proci=0; proci<nProcs; proci++)
1008 {
1009 Info<< "Reading mesh to add from "
1010 << databases[proci].caseName()
1011 << " for time = " << databases[proci].timeName()
1012 << nl << nl << endl;
1013 }
1014
1015 // See if any points on the mastermesh have become connected
1016 // because of connections through processor meshes.
1017 mergeSharedPoints(mergeDist,masterMesh[0],pointProcAddressing);
1018
1019 // Save some properties on the reconstructed mesh
1020 masterInternalFaces = masterMesh[0].nInternalFaces();
1021 masterOwner = masterMesh[0].faceOwner();
1022
1023
1024 if (writeAddrOnly)
1025 {
1026 Info<< nl
1027 << "Disabled writing of merged mesh (-addressing-only)"
1028 << nl << nl;
1029 }
1030 else
1031 {
1032 // Write reconstructed mesh
1033 writeMesh(masterMesh[0], cellProcAddressing);
1034 if (writeCellDist)
1035 {
1036 writeDistribution
1037 (
1038 runTime,
1039 masterMesh[0],
1040 cellProcAddressing
1041 );
1042 }
1043 }
1044 }
1045 else
1046 {
1047 // Load all meshes
1048 PtrList<fvMesh> fvMeshes(nProcs);
1049 for (label proci=0; proci<nProcs; proci++)
1050 {
1051 fvMeshes.set
1052 (
1053 proci,
1054 new fvMesh
1055 (
1056 IOobject
1057 (
1058 regionName,
1059 databases[proci].timeName(),
1060 databases[proci]
1061 )
1062 )
1063 );
1064 }
1065
1066 // Construct pointers to meshes
1067 UPtrList<polyMesh> meshes(fvMeshes.size());
1068 forAll(fvMeshes, proci)
1069 {
1070 meshes.set(proci, &fvMeshes[proci]);
1071 }
1072
1073 // Get pairs of patches to stitch. These pairs have to
1074 // - have ordered, opposite faces (so one to one correspondence)
1075 List<DynamicList<label>> localPatch;
1076 List<DynamicList<label>> remoteProc;
1077 List<DynamicList<label>> remotePatch;
1078 const label nGlobalPatches = polyMeshAdder::procPatchPairs
1079 (
1080 meshes,
1081 localPatch,
1082 remoteProc,
1083 remotePatch
1084 );
1085
1086 // Collect matching boundary faces on patches-to-stitch
1087 labelListList localBoundaryFace;
1088 labelListList remoteFaceProc;
1089 labelListList remoteBoundaryFace;
1090 polyMeshAdder::patchFacePairs
1091 (
1092 meshes,
1093 localPatch,
1094 remoteProc,
1095 remotePatch,
1096 localBoundaryFace,
1097 remoteFaceProc,
1098 remoteBoundaryFace
1099 );
1100
1101 // All matched faces assumed to have vertex0 matched
1102 labelListList remoteFaceStart(meshes.size());
1103 forAll(meshes, proci)
1104 {
1105 const labelList& procFaces = localBoundaryFace[proci];
1106 remoteFaceStart[proci].setSize(procFaces.size(), 0);
1107 }
1108
1109
1110
1111 labelListList patchMap(meshes.size());
1112 labelListList pointZoneMap(meshes.size());
1113 labelListList faceZoneMap(meshes.size());
1114 labelListList cellZoneMap(meshes.size());
1115 forAll(meshes, proci)
1116 {
1117 const polyMesh& mesh = meshes[proci];
1118 patchMap[proci] = identity(mesh.boundaryMesh().size());
1119
1120 // Remove excess patches
1121 patchMap[proci].setSize(nGlobalPatches);
1122
1123 pointZoneMap[proci] = identity(mesh.pointZones().size());
1124 faceZoneMap[proci] = identity(mesh.faceZones().size());
1125 cellZoneMap[proci] = identity(mesh.cellZones().size());
1126 }
1127
1128
1129 refPtr<fvMesh> masterMeshPtr;
1130 {
1131 // Do in-place addition on proc0.
1132
1133 const labelList oldFaceOwner(fvMeshes[0].faceOwner());
1134
1136 (
1137 0, // index of mesh to modify (== mesh_)
1138 fvMeshes,
1139 oldFaceOwner,
1140
1141 // Coupling info
1142 localBoundaryFace,
1143 remoteFaceProc,
1144 remoteBoundaryFace,
1145
1146 boundProcAddressing,
1147 cellProcAddressing,
1148 faceProcAddressing,
1149 pointProcAddressing
1150 );
1151
1152 // Remove zero-faces processor patches
1153 const polyBoundaryMesh& pbm = fvMeshes[0].boundaryMesh();
1154 labelList oldToNew(pbm.size(), -1);
1155 label newi = 0;
1156 // Non processor patches first
1157 forAll(pbm, patchi)
1158 {
1159 const auto& pp = pbm[patchi];
1160 if (!isA<processorPolyPatch>(pp) || pp.size())
1161 {
1162 oldToNew[patchi] = newi++;
1163 }
1164 }
1165 const label nNonProcPatches = newi;
1166
1167 // Move all deletable patches to the end
1168 forAll(oldToNew, patchi)
1169 {
1170 if (oldToNew[patchi] == -1)
1171 {
1172 oldToNew[patchi] = newi++;
1173 }
1174 }
1175 fvMeshTools::reorderPatches
1176 (
1177 fvMeshes[0],
1178 oldToNew,
1179 nNonProcPatches,
1180 false
1181 );
1182
1183 masterMeshPtr = fvMeshes[0];
1184 }
1185
1186
1187 const fvMesh& masterMesh = masterMeshPtr();
1188
1189 // Number of internal faces on the final reconstructed mesh
1190 masterInternalFaces = masterMesh.nInternalFaces();
1191
1192 // Owner addressing on the final reconstructed mesh
1193 masterOwner = masterMesh.faceOwner();
1194
1195
1196 // Write reconstructed mesh
1197 // Override:
1198 // - caseName
1199 // - processorCase flag
1200 // so the resulting mesh goes to the correct location (even with
1201 // collated). The better way of solving this is to construct
1202 // (zero) mesh on the undecomposed runTime.
1203
1204 if (writeAddrOnly)
1205 {
1206 Info<< nl
1207 << "Disabled writing of merged mesh (-addressing-only)"
1208 << nl << nl;
1209 }
1210 else
1211 {
1212 Time& masterTime = const_cast<Time&>(masterMesh.time());
1213
1214 const word oldCaseName = masterTime.caseName();
1215 masterTime.caseName() = runTime.caseName();
1216 const bool oldProcCase(masterTime.processorCase(false));
1217
1218 // Write reconstructed mesh
1219 writeMesh(masterMesh, cellProcAddressing);
1220 if (writeCellDist)
1221 {
1222 writeDistribution
1223 (
1224 runTime,
1225 masterMesh,
1226 cellProcAddressing
1227 );
1228 }
1229 masterTime.caseName() = oldCaseName;
1230 masterTime.processorCase(oldProcCase);
1231 }
1232 }
1233
1234
1235 // Write the addressing
1236
1237 Info<< "Reconstructing addressing from processor meshes"
1238 << " to the newly reconstructed mesh" << nl << endl;
1239
1240 forAll(databases, proci)
1241 {
1242 Info<< "Processor " << proci << nl
1243 << "Read processor mesh: "
1244 << (databases[proci].caseName()/regionDir) << endl;
1245
1246 polyMesh procMesh
1247 (
1248 IOobject
1249 (
1250 regionName,
1251 databases[proci].timeName(),
1252 databases[proci]
1253 )
1254 );
1255
1256 writeMaps
1257 (
1258 masterInternalFaces,
1259 masterOwner,
1260 procMesh,
1261 cellProcAddressing[proci],
1262 faceProcAddressing[proci],
1263 pointProcAddressing[proci],
1264 boundProcAddressing[proci]
1265 );
1266 }
1267 }
1268 }
1269
1270 Info<< "End\n" << endl;
1271
1272 return 0;
1273}
1274
1275
1276// ************************************************************************* //
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
bool set(const Key &key, const T &obj)
Copy assign a new entry, overwriting existing entries.
Definition: HashTableI.H:202
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
A IOList wrapper for writing external data.
Definition: IOList.H:123
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 HashTable to objects of type <T> with a label key.
Definition: Map.H:60
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
const T * set(const label i) const
Definition: PtrList.H:138
bool processorCase() const noexcept
Return true if this is a processor case.
Definition: TimePathsI.H:36
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
const fileName & caseName() const
Return case name.
Definition: TimePathsI.H:62
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
A list of pointers to objects of type <T>, without allocation/deallocation management of the pointers...
Definition: UPtrList.H:71
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
const fileName & rootPath() const noexcept
Return root path.
Definition: argListI.H:63
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
fileName path() const
Return the full path to the (processor local) case.
Definition: argListI.H:81
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
const fileName & caseName() const noexcept
Return case name (parallel run) or global case (serial run)
Definition: argListI.H:69
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
void add(const boundBox &bb)
Extend to include the second box.
Definition: boundBoxI.H:191
scalar mag() const
The magnitude of the bounding box span.
Definition: boundBoxI.H:133
A class for handling file names.
Definition: fileName.H:76
virtual label nProcs(const fileName &dir, const fileName &local="") const
Get number of processor directories/results. Used for e.g.
Sums a given list of (at least two or more) fields and outputs the result into a new field,...
Definition: add.H:161
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
const Time & time() const noexcept
Return time registry.
const word & name() const noexcept
The patch name.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
static const word & regionName(const word &region)
The mesh region name or word::null if polyMesh::defaultRegion.
Definition: polyMesh.C:829
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:866
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:364
Direct mesh changes based on v1.3 polyTopoChange syntax.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
label nInternalFaces() const noexcept
Number of internal faces.
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
A class for managing references or pointers (no reference counting)
Definition: refPtr.H:58
virtual bool write(const bool valid=true) const
Write using setting from DB.
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
engineTime & runTime
Info<< "Creating field kinetic energy K\n"<< endl;volScalarField K("K", 0.5 *magSqr(U));if(U.nOldTimes()){ volVectorField *Uold=&U.oldTime();volScalarField *Kold=&K.oldTime();*Kold==0.5 *magSqr(*Uold);while(Uold->nOldTimes()) { Uold=&Uold-> oldTime()
Foam::word regionName(Foam::polyMesh::defaultRegion)
Foam::PtrList< Foam::fvMesh > meshes(regionNames.size())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const word & regionDir
wordList regionNames
const pointField & points
word timeName
Definition: getTimeIndex.H:3
Namespace for OpenFOAM.
const fileOperation & fileHandler()
Get current file handler.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
messageStream Info
Information stream (stdout output on master, null elsewhere)
IntListType renumber(const labelUList &oldToNew, const IntListType &input)
Renumber the values (not the indices) of a list.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label newPointi
Definition: readKivaGrid.H:496
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
static const char *const typeName
The type name used in ensight case files.