From 76c16434abf2c47d541b9f1125fec3c5b8e86175 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 5 Oct 2021 21:30:17 +0200
Subject: [PATCH] BUG: incorrect finite-area faceProcAddressing (fixes #2155)

- when creating a finite-area mesh in parallel, need to determine
  the equivalent ProcAddressing for the faMesh.

  In the faceProcAddressing the collected and sorted order was being
  scattered directly back to the individual processors instead of only
  the sections relevant to each particular processor.

  This caused the observed jumbled order for reconstructed fields.
---
 .../faReconstruct/faMeshReconstructor.C       | 92 +++++++++++++++++--
 1 file changed, 84 insertions(+), 8 deletions(-)

diff --git a/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C b/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C
index 0f85a81d5e2..67e41128620 100644
--- a/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C
+++ b/src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C
@@ -88,11 +88,85 @@ void Foam::faMeshReconstructor::calcAddressing
     {
         globalFaceNum.gather(faFaceProcAddr_, singlePatchFaceLabels_);
 
-        labelList order(Foam::sortedOrder(singlePatchFaceLabels_));
+        const labelList globalOrder(Foam::sortedOrder(singlePatchFaceLabels_));
 
-        singlePatchFaceLabels_ = labelList(singlePatchFaceLabels_, order);
+        singlePatchFaceLabels_ =
+            labelList(singlePatchFaceLabels_, globalOrder);
 
-        globalFaceNum.scatter(order, faFaceProcAddr_);
+        // Set first face to be zero relative to the finiteArea patch
+        // ie, local-face numbering with the first being 0 on any given patch
+        {
+            label patchFirstMeshfacei
+            (
+                singlePatchFaceLabels_.empty()
+              ? 0
+              : singlePatchFaceLabels_.first()
+            );
+            Pstream::scatter(patchFirstMeshfacei);
+
+            for (label& facei : faFaceProcAddr_)
+            {
+                facei -= patchFirstMeshfacei;
+            }
+        }
+
+        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+        if (Pstream::master())
+        {
+            // Determine the respective local portions of the global ordering
+
+            labelList procTargets(globalFaceNum.size());
+
+            for (const label proci : Pstream::allProcs())
+            {
+                labelList::subList
+                (
+                    globalFaceNum.range(proci),
+                    procTargets
+                ) = proci;
+            }
+
+            labelList procStarts(globalFaceNum.offsets());
+            labelList procOrders(globalFaceNum.size());
+
+            for (const label globali : globalOrder)
+            {
+                const label proci = procTargets[globali];
+
+                procOrders[procStarts[proci]++] =
+                    (globali - globalFaceNum.localStart(proci));
+            }
+
+            // Send the local portions
+            for (const int proci : Pstream::subProcs())
+            {
+                SubList<label> localOrder
+                (
+                    procOrders,
+                    globalFaceNum.range(proci)
+                );
+
+                UOPstream toProc(proci, pBufs);
+                toProc << localOrder;
+            }
+
+            SubList<label> localOrder(procOrders, globalFaceNum.range(0));
+
+            faFaceProcAddr_ = labelList(faFaceProcAddr_, localOrder);
+        }
+
+        pBufs.finishedSends();
+
+        if (!Pstream::master())
+        {
+            labelList localOrder;
+
+            UIPstream fromProc(Pstream::master(), pBufs);
+            fromProc >> localOrder;
+
+            faFaceProcAddr_ = labelList(faFaceProcAddr_, localOrder);
+        }
     }
 
     // Broadcast the same information everywhere
@@ -201,12 +275,14 @@ void Foam::faMeshReconstructor::calcAddressing
         tmpFaces = initialPatch.localFaces();
         pointField tmpPoints(initialPatch.localPoints());
 
-        // The meshPointMap is contiguous, so flatten as linear list
-        /// Map<label> mpm(initialPatch.meshPointMap());
-        labelList mpm(initialPatch.nPoints());
-        forAllConstIters(initialPatch.meshPointMap(), iter)
+        // Equivalent to a flattened meshPointMap
+        labelList mpm(initialPatch.points().size(), -1);
         {
-            mpm[iter.key()] = iter.val();
+            const labelList& mp = initialPatch.meshPoints();
+            forAll(mp, i)
+            {
+                mpm[mp[i]] = i;
+            }
         }
         Pstream::scatter(mpm);
 
-- 
GitLab