ENH: multi-region support for reconstructParMesh (#2072)

applications/test/checkDecomposePar/Test-checkDecomposePar.C

@@ -71,7 +71,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
 
     const auto decompFile = args.get<fileName>(1);
-    const bool region     = args.found("region");
     const bool verbose    = args.found("verbose");
 
     // Set time from database
@@ -87,23 +86,17 @@ int main(int argc, char *argv[])
     // Get region names
     #include "getAllRegionOptions.H"
 
-    wordList regionDirs(regionNames);
-    if (regionDirs.size() == 1 && regionDirs[0] == polyMesh::defaultRegion)
-    {
-        regionDirs[0].clear();
-    }
-    else
-    {
-        Info<< "Decomposing regions: "
-            << flatOutput(regionNames) << nl << endl;
-    }
-
     labelList cellToProc;
 
     forAll(regionNames, regioni)
     {
         const word& regionName = regionNames[regioni];
-        const word& regionDir = regionDirs[regioni];
+        // const word& regionDir =
+        // (
+        //     regionName != polyMesh::defaultRegion
+        //   ? regionName
+        //   : word::null
+        // );
 
         Info<< "\n\nDecomposing mesh " << regionName << nl << endl;
         Info<< "Create mesh..." << flush;

applications/utilities/parallelProcessing/decomposePar/domainDecompositionWrite.C

@@ -52,7 +52,7 @@ void Foam::domainDecomposition::writeVolField
             false
         ),
         this->mesh(),
-        dimensionedScalar("cellDist", dimless, -1),
+        dimensionedScalar(dimless, Zero),
         zeroGradientFvPatchScalarField::typeName
     );
 

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -172,23 +172,29 @@ int main(int argc, char *argv[])
     // Get region names
     #include "getAllRegionOptions.H"
 
-    wordList regionDirs(regionNames);
+    // Determine the processor count
+    label nProcs{0};
 
-    if (regionNames.size() == 1)
+    if (regionNames.empty())
     {
-        if (regionNames[0] == polyMesh::defaultRegion)
+        FatalErrorInFunction
+            << "No regions specified or detected."
+            << exit(FatalError);
+    }
+    else if (regionNames[0] == polyMesh::defaultRegion)
     {
-            regionDirs[0].clear();
+        nProcs = fileHandler().nProcs(args.path());
     }
     else
+    {
+        nProcs = fileHandler().nProcs(args.path(), regionNames[0]);
+
+        if (regionNames.size() == 1)
         {
             Info<< "Using region: " << regionNames[0] << nl << endl;
         }
     }
 
-    // Determine the processor count
-    label nProcs = fileHandler().nProcs(args.path(), regionDirs[0]);
-
     if (!nProcs)
     {
         FatalErrorInFunction
@@ -261,16 +267,21 @@ int main(int argc, char *argv[])
     forAll(regionNames, regioni)
     {
         const word& regionName = regionNames[regioni];
-        const word& regionDir = regionDirs[regioni];
+        const word& regionDir =
+        (
+            regionName != polyMesh::defaultRegion
+          ? regionName
+          : word::null
+        );
 
-        Info<< "\n\nReconstructing fields for mesh " << regionName << nl
+        Info<< "\n\nReconstructing fields" << nl
-            << endl;
+            << "region=" << regionName << nl << endl;
 
         if
         (
             newTimes
          && regionNames.size() == 1
-         && regionDirs[0].empty()
+         && regionDir.empty()
          && haveAllTimes(masterTimeDirSet, timeDirs)
         )
         {

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -71,6 +71,7 @@ Usage
 #include "polyTopoChange.H"
 #include "extrapolatedCalculatedFvPatchFields.H"
 #include "topoSet.H"
+#include "regionProperties.H"
 #include "fvMeshTools.H"
 
 using namespace Foam;
@@ -304,31 +305,29 @@ autoPtr<mapPolyMesh> mergeSharedPoints
 
 boundBox procBounds
 (
-    const argList& args,
     const PtrList<Time>& databases,
     const word& regionDir
 )
 {
-    boundBox bb = boundBox::invertedBox;
+    boundBox bb;
 
-    forAll(databases, proci)
+    for (const Time& procDb : databases)
     {
         fileName pointsInstance
         (
-            databases[proci].findInstance
+            procDb.findInstance
             (
                 regionDir/polyMesh::meshSubDir,
                 "points"
             )
         );
 
-        if (pointsInstance != databases[proci].timeName())
+        if (pointsInstance != procDb.timeName())
         {
             FatalErrorInFunction
-                << "Your time was specified as " << databases[proci].timeName()
+                << "Your time was specified as " << procDb.timeName()
-                << " but there is no polyMesh/points in that time." << endl
+                << " but there is no polyMesh/points in that time." << nl
-                << "(there is a points file in " << pointsInstance
+                << "(points file at " << pointsInstance << ')' << nl
-                << ")" << endl
                 << "Please rerun with the correct time specified"
                 << " (through the -constant, -time or -latestTime "
                 << "(at your option)."
@@ -336,8 +335,8 @@ boundBox procBounds
         }
 
         Info<< "Reading points from "
-            << databases[proci].caseName()
+            << procDb.caseName()
-            << " for time = " << databases[proci].timeName()
+            << " for time = " << procDb.timeName()
             << nl << endl;
 
         pointIOField points
@@ -345,13 +344,13 @@ boundBox procBounds
             IOobject
             (
                 "points",
-                databases[proci].findInstance
+                procDb.findInstance
                 (
                     regionDir/polyMesh::meshSubDir,
                     "points"
                 ),
                 regionDir/polyMesh::meshSubDir,
-                databases[proci],
+                procDb,
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE,
                 false
@@ -397,10 +396,9 @@ void writeDistribution
     cellDecomposition.write();
 
     Info<< nl << "Wrote decomposition to "
-        << cellDecomposition.objectPath()
+        << runTime.relativePath(cellDecomposition.objectPath())
         << " for use in manual decomposition." << endl;
 
-
     // Write as volScalarField for postprocessing. Change time to 0
     // if was 'constant'
     {
@@ -419,7 +417,8 @@ void writeDistribution
                 runTime.timeName(),
                 masterMesh,
                 IOobject::NO_READ,
-                IOobject::AUTO_WRITE
+                IOobject::NO_WRITE,
+                false
             ),
             masterMesh,
             dimensionedScalar(dimless, Zero),
@@ -430,12 +429,13 @@ void writeDistribution
         {
             cellDist[celli] = cellDecomposition[celli];
         }
-        cellDist.correctBoundaryConditions();
 
+        cellDist.correctBoundaryConditions();
         cellDist.write();
 
-        Info<< nl << "Wrote decomposition as volScalarField to "
+        Info<< nl << "Wrote decomposition to "
-            << cellDist.name() << " for use in postprocessing."
+            << runTime.relativePath(cellDist.objectPath())
+            << " (volScalarField) for visualization."
             << endl;
 
         // Restore time
@@ -450,9 +450,23 @@ void writeMesh
     const labelListList& cellProcAddressing
 )
 {
-    Info<< "\nWriting merged mesh to "
+    const word& regionDir =
-        << mesh.time().path()/mesh.time().timeName()
+    (
-        << nl << endl;
+        mesh.name() == polyMesh::defaultRegion
+      ? word::null
+      : mesh.name()
+    );
+
+    const fileName outputDir
+    (
+        mesh.time().path()
+      / mesh.time().timeName()
+      / regionDir
+      / polyMesh::meshSubDir
+    );
+
+    Info<< nl << "Writing merged mesh to "
+        << mesh.time().relativePath(outputDir) << nl << endl;
 
     if (!mesh.write())
     {
@@ -472,13 +486,21 @@ void writeMaps
     const labelUList& cellProcAddressing,
     const labelUList& faceProcAddressing,
     const labelUList& pointProcAddressing,
-    const labelUList& boundaryProcAddressing
+    const labelUList& boundProcAddressing
 )
 {
+    const word& regionDir =
+    (
+        procMesh.name() == polyMesh::defaultRegion
+      ? word::null
+      : procMesh.name()
+    );
+
     const fileName outputDir
     (
         procMesh.time().caseName()
       / procMesh.facesInstance()
+      / regionDir
       / polyMesh::meshSubDir
     );
 
@@ -494,16 +516,18 @@ void writeMaps
     );
 
 
+    Info<< "Writing addressing : " << outputDir << nl;
+
     // From processor point to reconstructed mesh point
 
-    Info<< "Writing pointProcAddressing to " << outputDir << endl;
+    Info<< "    pointProcAddressing" << endl;
     ioAddr.rename("pointProcAddressing");
     labelIOList(ioAddr, pointProcAddressing).write();
 
 
     // From processor face to reconstructed mesh face
 
-    Info<< "Writing faceProcAddressing to " << outputDir << endl;
+    Info<< "    faceProcAddressing" << endl;
     ioAddr.rename("faceProcAddressing");
 
     labelIOList faceProcAddr(ioAddr, faceProcAddressing);
@@ -549,26 +573,39 @@ void writeMaps
 
     // From processor cell to reconstructed mesh cell
 
-    Info<< "Writing cellProcAddressing to " << outputDir << endl;
+    Info<< "    cellProcAddressing" << endl;
     ioAddr.rename("cellProcAddressing");
     labelIOList(ioAddr, cellProcAddressing).write();
 
 
     // From processor patch to reconstructed mesh patch
 
-    Info<< "Writing boundaryProcAddressing to " << outputDir << endl;
+    Info<< "    boundaryProcAddressing" << endl;
     ioAddr.rename("boundaryProcAddressing");
-    labelIOList(ioAddr, boundaryProcAddressing).write();
+    labelIOList(ioAddr, boundProcAddressing).write();
 
     Info<< endl;
 }
 
 
+void printWarning()
+{
+    Info<<
+"Merge individual processor meshes back into one master mesh.\n"
+"Use if the original master mesh has been deleted or the processor meshes\n"
+"have been modified (topology change).\n"
+"This tool will write the resulting mesh to a new time step and construct\n"
+"xxxxProcAddressing files in the processor meshes so reconstructPar can be\n"
+"used to regenerate the fields on the master mesh.\n\n"
+"Not well tested & use at your own risk!\n\n";
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::addNote
     (
-        "Reconstruct a mesh using geometric information only"
+        "Reconstruct a mesh using geometric/topological information only"
     );
 
     // Enable -constant ... if someone really wants it
@@ -576,6 +613,12 @@ int main(int argc, char *argv[])
     timeSelector::addOptions(true, true); // constant(true), zero(true)
 
     argList::noParallel();
+
+    argList::addBoolOption
+    (
+        "addressing-only",
+        "Create procAddressing only without overwriting the mesh"
+    );
     argList::addOption
     (
         "mergeTol",
@@ -599,61 +642,40 @@ int main(int argc, char *argv[])
         "decomposition method or as a volScalarField for post-processing."
     );
 
-    #include "addRegionOption.H"
+    #include "addAllRegionOptions.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
 
-    Info<< "This is an experimental tool which tries to merge"
+    printWarning();
-        << " individual processor" << nl
-        << "meshes back into one master mesh. Use it if the original"
-        << " master mesh has" << nl
-        << "been deleted or if the processor meshes have been modified"
-        << " (topology change)." << nl
-        << "This tool will write the resulting mesh to a new time step"
-        << " and construct" << nl
-        << "xxxxProcAddressing files in the processor meshes so"
-        << " reconstructPar can be" << nl
-        << "used to regenerate the fields on the master mesh." << nl
-        << nl
-        << "Not well tested & use at your own risk!" << nl
-        << endl;
-
-
-    word regionName(polyMesh::defaultRegion);
-    word regionDir;
-
-    if
-    (
-        args.readIfPresent("region", regionName)
-     && regionName != polyMesh::defaultRegion
-    )
-    {
-        regionDir = regionName;
-        Info<< "Operating on region " << regionName << nl << endl;
-    }
 
     const bool fullMatch = args.found("fullMatch");
     const bool procMatch = args.found("procMatch");
     const bool writeCellDist = args.found("cellDist");
 
+    const bool writeAddrOnly = args.found("addressing-only");
+
+    const scalar mergeTol =
+        args.getOrDefault<scalar>("mergeTol", defaultMergeTol);
+
     if (fullMatch)
     {
-        Info<< "Doing geometric matching on all boundary faces." << nl << endl;
+        Info<< "Use geometric matching on all boundary faces." << nl << endl;
     }
     else if (procMatch)
     {
-        Info<< "Doing geometric matching on correct procBoundaries only."
+        Info<< "Use geometric matching on correct procBoundaries only." << nl
-            << nl << "This assumes a correct decomposition." << endl;
+            << "This assumes a correct decomposition." << endl;
     }
     else
     {
-        Info<< "Assuming correct, fully matched procBoundaries." << nl << endl;
+        Info<< "Merge assuming correct, fully matched procBoundaries." << nl
+            << endl;
     }
 
-    scalar mergeTol = args.getOrDefault<scalar>("mergeTol", defaultMergeTol);
     if (fullMatch || procMatch)
     {
-        scalar writeTol =
+        const scalar writeTol =
             Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
 
         Info<< "Merge tolerance : " << mergeTol << nl
@@ -673,16 +695,46 @@ int main(int argc, char *argv[])
         }
     }
 
-    label nProcs = fileHandler().nProcs(args.path());
+    // Get region names
+    #include "getAllRegionOptions.H"
 
-    Info<< "Found " << nProcs << " processor directories" << nl << endl;
+    // Determine the processor count
+    label nProcs{0};
+
+    if (regionNames.empty())
+    {
+        FatalErrorInFunction
+            << "No regions specified or detected."
+            << exit(FatalError);
+    }
+    else if (regionNames[0] == polyMesh::defaultRegion)
+    {
+        nProcs = fileHandler().nProcs(args.path());
+    }
+    else
+    {
+        nProcs = fileHandler().nProcs(args.path(), regionNames[0]);
+
+        if (regionNames.size() == 1)
+        {
+            Info<< "Using region: " << regionNames[0] << nl << endl;
+        }
+    }
+
+    if (!nProcs)
+    {
+        FatalErrorInFunction
+            << "No processor* directories found"
+            << exit(FatalError);
+    }
 
     // Read all time databases
     PtrList<Time> databases(nProcs);
 
+    Info<< "Found " << nProcs << " processor directories" << endl;
     forAll(databases, proci)
     {
-        Info<< "Reading database "
+        Info<< "    Reading database "
             << args.caseName()/("processor" + Foam::name(proci))
             << endl;
 
@@ -697,6 +749,7 @@ int main(int argc, char *argv[])
             )
         );
     }
+    Info<< endl;
 
     // Use the times list from the master processor
     // and select a subset based on the command-line options
@@ -707,19 +760,32 @@ int main(int argc, char *argv[])
     );
 
     // Loop over all times
-    forAll(timeDirs, timeI)
+    forAll(timeDirs, timei)
     {
         // Set time for global database
-        runTime.setTime(timeDirs[timeI], timeI);
+        runTime.setTime(timeDirs[timei], timei);
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // Set time for all databases
         forAll(databases, proci)
         {
-            databases[proci].setTime(timeDirs[timeI], timeI);
+            databases[proci].setTime(timeDirs[timei], timei);
         }
 
+        Info<< "Time = " << runTime.timeName() << endl;
+
+        // Check for any mesh changes
+        label nMeshChanged = 0;
+        boolList hasRegionMesh(regionNames.size(), false);
+        forAll(regionNames, regioni)
+        {
+            const word& regionName = regionNames[regioni];
+            const word& regionDir =
+            (
+                regionName == polyMesh::defaultRegion
+              ? word::null
+              : regionName
+            );
+
             IOobject facesIO
             (
                 "faces",
@@ -730,22 +796,55 @@ int main(int argc, char *argv[])
                 IOobject::NO_WRITE
             );
 
-
             // Problem: faceCompactIOList recognises both 'faceList' and
             //          'faceCompactList' so we should be lenient when doing
             //          typeHeaderOk
-        if (!facesIO.typeHeaderOk<faceCompactIOList>(false))
+
+            const bool ok = facesIO.typeHeaderOk<faceCompactIOList>(false);
+
+            if (ok)
             {
-            Info<< "No mesh." << nl << endl;
+                hasRegionMesh[regioni] = true;
+                ++nMeshChanged;
+            }
+        }
+
+        // Check for any mesh changes
+        if (!nMeshChanged)
+        {
+            Info<< "No meshes" << nl << endl;
             continue;
         }
 
+        Info<< endl;
+
+        forAll(regionNames, regioni)
+        {
+            const word& regionName = regionNames[regioni];
+            const word& regionDir =
+            (
+                regionName == polyMesh::defaultRegion
+              ? word::null
+              : regionName
+            );
+
+            if (!hasRegionMesh[regioni])
+            {
+                Info<< "region=" << regionName << " (no mesh)" << nl << endl;
+                continue;
+            }
+
+            if (regionNames.size() > 1)
+            {
+                Info<< "region=" << regionName << nl;
+            }
+
 
             // Addressing from processor to reconstructed case
             labelListList cellProcAddressing(nProcs);
             labelListList faceProcAddressing(nProcs);
             labelListList pointProcAddressing(nProcs);
-        labelListList boundaryProcAddressing(nProcs);
+            labelListList boundProcAddressing(nProcs);
 
 
             // Internal faces on the final reconstructed mesh
@@ -756,10 +855,11 @@ int main(int argc, char *argv[])
 
             if (procMatch)
             {
-            // Read point on individual processors to determine merge tolerance
+                // Read point on individual processors to determine
+                // merge tolerance
                 // (otherwise single cell domains might give problems)
 
-            const boundBox bb = procBounds(args, databases, regionDir);
+                const boundBox bb = procBounds(databases, regionDir);
                 const scalar mergeDist = mergeTol*bb.mag();
 
                 Info<< "Overall mesh bounding box  : " << bb << nl
@@ -769,7 +869,6 @@ int main(int argc, char *argv[])
 
 
                 // Construct empty mesh.
-            // fvMesh** masterMesh = new fvMesh*[nProcs];
                 PtrList<fvMesh> masterMesh(nProcs);
 
                 for (label proci=0; proci<nProcs; proci++)
@@ -804,7 +903,7 @@ int main(int argc, char *argv[])
                     cellProcAddressing[proci] = identity(meshToAdd.nCells());
                     faceProcAddressing[proci] = identity(meshToAdd.nFaces());
                     pointProcAddressing[proci] = identity(meshToAdd.nPoints());
-                boundaryProcAddressing[proci] =
+                    boundProcAddressing[proci] =
                         identity(meshToAdd.boundaryMesh().size());
 
                     // Find geometrically shared points/faces.
@@ -832,7 +931,7 @@ int main(int argc, char *argv[])
                     renumber(map().addedCellMap(), cellProcAddressing[proci]);
                     renumber(map().addedFaceMap(), faceProcAddressing[proci]);
                     renumber(map().addedPointMap(), pointProcAddressing[proci]);
-                renumber(map().addedPatchMap(), boundaryProcAddressing[proci]);
+                    renumber(map().addedPatchMap(), boundProcAddressing[proci]);
                 }
                 for (label step=2; step<nProcs*2; step*=2)
                 {
@@ -844,7 +943,8 @@ int main(int argc, char *argv[])
                             continue;
                         }
 
-                    Info<< "Merging mesh " << proci << " with " << next << endl;
+                        Info<< "Merging mesh " << proci << " with "
+                            << next << endl;
 
                         // Find geometrically shared points/faces.
                         autoPtr<faceCoupleInfo> couples = determineCoupledFaces
@@ -892,7 +992,7 @@ int main(int argc, char *argv[])
                             renumber
                             (
                                 map().oldPatchMap(),
-                            boundaryProcAddressing[mergedI]
+                                boundProcAddressing[mergedI]
                             );
                         }
 
@@ -925,7 +1025,7 @@ int main(int argc, char *argv[])
                             renumber
                             (
                                 map().addedPatchMap(),
-                            boundaryProcAddressing[addedI]
+                                boundProcAddressing[addedI]
                             );
                         }
 
@@ -943,17 +1043,32 @@ int main(int argc, char *argv[])
 
                 // See if any points on the mastermesh have become connected
                 // because of connections through processor meshes.
-            mergeSharedPoints(mergeDist, masterMesh[0], pointProcAddressing);
+                mergeSharedPoints(mergeDist,masterMesh[0],pointProcAddressing);
 
                 // Save some properties on the reconstructed mesh
                 masterInternalFaces = masterMesh[0].nInternalFaces();
                 masterOwner = masterMesh[0].faceOwner();
 
+
+                if (writeAddrOnly)
+                {
+                    Info<< nl
+                        << "Disabled writing of merged mesh (-addressing-only)"
+                        << nl << nl;
+                }
+                else
+                {
                     // Write reconstructed mesh
                     writeMesh(masterMesh[0], cellProcAddressing);
                     if (writeCellDist)
                     {
-                writeDistribution(runTime, masterMesh[0], cellProcAddressing);
+                        writeDistribution
+                        (
+                            runTime,
+                            masterMesh[0],
+                            cellProcAddressing
+                        );
+                    }
                 }
             }
             else
@@ -1057,7 +1172,7 @@ int main(int argc, char *argv[])
                         remoteFaceProc,
                         remoteBoundaryFace,
 
-                    boundaryProcAddressing,
+                        boundProcAddressing,
                         cellProcAddressing,
                         faceProcAddressing,
                         pointProcAddressing
@@ -1102,9 +1217,11 @@ int main(int argc, char *argv[])
 
                 // Number of internal faces on the final reconstructed mesh
                 masterInternalFaces = masterMesh.nInternalFaces();
+
                 // Owner addressing on the final reconstructed mesh
                 masterOwner = masterMesh.faceOwner();
 
+
                 // Write reconstructed mesh
                 // Override:
                 //  - caseName
@@ -1112,32 +1229,48 @@ int main(int argc, char *argv[])
                 // so the resulting mesh goes to the correct location (even with
                 // collated). The better way of solving this is to construct
                 // (zero) mesh on the undecomposed runTime.
+
+                if (writeAddrOnly)
+                {
+                    Info<< nl
+                        << "Disabled writing of merged mesh (-addressing-only)"
+                        << nl << nl;
+                }
+                else
+                {
                     Time& masterTime = const_cast<Time&>(masterMesh.time());
 
                     const word oldCaseName = masterTime.caseName();
                     masterTime.caseName() = runTime.caseName();
                     const bool oldProcCase(masterTime.processorCase(false));
 
+                    // Write reconstructed mesh
                     writeMesh(masterMesh, cellProcAddressing);
                     if (writeCellDist)
                     {
-                writeDistribution(runTime, masterMesh, cellProcAddressing);
+                        writeDistribution
+                        (
+                            runTime,
+                            masterMesh,
+                            cellProcAddressing
+                        );
                     }
                     masterTime.caseName() = oldCaseName;
                     masterTime.processorCase(oldProcCase);
                 }
-
+            }
 
 
             // Write the addressing
 
-        Info<< "Reconstructing the addressing from the processor meshes"
+            Info<< "Reconstructing addressing from processor meshes"
                 << " to the newly reconstructed mesh" << nl << endl;
 
             forAll(databases, proci)
             {
-            Info<< "Reading processor " << proci << " mesh from "
+                Info<< "Processor " << proci << nl
-                << databases[proci].caseName() << endl;
+                    << "Read processor mesh: "
+                    << (databases[proci].caseName()/regionDir) << endl;
 
                 polyMesh procMesh
                 (
@@ -1157,13 +1290,13 @@ int main(int argc, char *argv[])
                     cellProcAddressing[proci],
                     faceProcAddressing[proci],
                     pointProcAddressing[proci],
-                boundaryProcAddressing[proci]
+                    boundProcAddressing[proci]
                 );
             }
         }
+    }
 
-
+    Info<< "End\n" << endl;
-    Info<< "\nEnd\n" << endl;
 
     return 0;
 }