Fixed bug in shared points for reconstruction

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C

@@ -290,7 +290,23 @@ Foam::processorMeshesReconstructor::neighbourProcPatch
             // Check neighbour processor index
             if (masterProcPatch.neighbProcNo() == procPatch.myProcNo())
             {
-                // Found matching patch
+                // Found matching patch.  Check sizes
+                if (masterProcPatch.size() != procPatch.size())
+                {
+                    FatalErrorIn
+                    (
+                        "const processorPolyPatch&\n"
+                        "processorMeshesReconstructor::neighbourProcPatch\n"
+                        "(\n"
+                        "    const processorPolyPatch& procPatch\n"
+                        ") const"
+                    )   << "Processor patch pair ("
+                        << procPatch.myProcNo() << " "
+                        << procPatch.neighbProcNo() << ") sizes do not match: "
+                        << masterProcPatch.size() << " and " << procPatch.size()
+                        << abort(FatalError);
+                }
+
                 return masterProcPatch;
             }
         }
@@ -528,7 +544,6 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     // Note: for easier debugging, set mapping, owner and neighbour to -1
     pointField reconPoints(nReconPoints);
-    labelList globalPointMapping(nReconPoints, -1);
     faceList reconFaces(nReconFaces);
     labelList cellOffset(meshes_.size(), 0);
     labelList reconOwner(nReconFaces, -1);
@@ -642,6 +657,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     // on the same processor
     sharedPoints sharedData(meshes_);
 
+    // Record global point index for shared points
+    labelList globalPointMapping(sharedData.nGlobalPoints(), -1);
+
     // Before assembling the meshes, create unique ordering for all passive
     // processor patches that will be merged.  HJ, 6/May/2018
 
@@ -655,7 +673,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     // Memory management
     {
-        // This list provides inser offset for the patch to unpack the list
+        // This list provides insert offset for the patch to unpack the list
         labelListList passivePatchInsertOffset(meshes_.size());
 
         // This list records all global indices from passive faces
@@ -789,12 +807,13 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
         }
 
         // Collect globally shared point labels
+        const labelList& curSpAddr = sharedData.sharedPointAddr()[fvmId];
         const labelList& curSpl = sharedData.sharedPointLabels()[fvmId];
 
-        forAll (curSpl, splI)
+        forAll (curSpAddr, splI)
         {
             // From processor 0, mark points without checking
-            globalPointMapping[curSpl[splI]] = ppAddr[curSpl[splI]];
+            globalPointMapping[curSpAddr[splI]] = ppAddr[curSpl[splI]];
         }
 
         // Dump all internal faces into the list
@@ -962,7 +981,8 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
         // Set cell offset for the next mesh
         if (cellOffset.size() > firstValidMesh() + 1)
         {
-            cellOffset[firstValidMesh() + 1] = meshes_[firstValidMesh()].nCells();
+            cellOffset[firstValidMesh() + 1] =
+                meshes_[firstValidMesh()].nCells();
         }
     }
 
@@ -973,7 +993,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     {
         if (meshes_.set(procI))
         {
-            Pout<< "Dump mesh " << procI << " cell offset: " << cellOffset[procI] << endl;
+            Pout<< "Dump mesh " << procI
+                << " cell offset: " << cellOffset[procI]
+                << endl;
 
             const polyMesh& curMesh = meshes_[procI];
             const polyBoundaryMesh& procPatches = curMesh.boundaryMesh();
@@ -992,10 +1014,29 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
             ppAddr.setSize(curPoints.size());
             ppAddr = -1;
 
+            // Collect globally shared point labels
+            const labelList& curSpAddr = sharedData.sharedPointAddr()[procI];
+            const labelList& curSpl = sharedData.sharedPointLabels()[procI];
+
+            // Go through globally shared points.  If already marked, mark
+            // their point index.  If not, leave for later
+
+            forAll (curSpAddr, splI)
+            {
+                // From other processors, check if point is already marked
+                // If yes, insert the point label into ppAddr
+                if (globalPointMapping[curSpAddr[splI]] > -1)
+                {
+                    // Point is marked.  Record it locally
+                    ppAddr[curSpl[splI]] = globalPointMapping[curSpAddr[splI]];
+                }
+            }
+
             // Find points already added via processor patches and mark them
             // in ppAddr
 
-            // Collect point-processor addressing for points on processor patches
+            // Collect point-processor addressing for points on
+            // processor patches
 
             // Go through all processor patches.  For neighbour patches, access
             // owner addressing and dump into ppAddr
@@ -1053,11 +1094,39 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                             const face& curMF = masterProcPatch[faceI];
 
                             forAll (curRF, pointI)
+                            {
+                                if (ppAddr[curRF[pointI]] == -1)
                                 {
                                     // Mapping is established
                                     ppAddr[curRF[pointI]] =
                                         masterPpAddr[curMF[pointI]];
                                 }
+                                else
+                                {
+                                    // Mapping already exists.  Check it
+                                    if
+                                    (
+                                        ppAddr[curRF[pointI]]
+                                     != masterPpAddr[curMF[pointI]]
+                                    )
+                                    {
+                                        WarningIn
+                                        (
+                                            "autoPtr<fvMesh> "
+                                            "processorMeshesReconstructor::"
+                                            "reconstructMesh(const Time& db)"
+                                        )   << "Loss of proc sync: "
+                                            << ppAddr[curRF[pointI]] << " and "
+                                            << masterPpAddr[curMF[pointI]]
+                                            << " for point "
+                                            << curRF[pointI]
+                                            << " at "
+                                            << curMesh.points()[curRF[pointI]]
+                                            << nl
+                                            << abort(FatalError);
+                                    }
+                                }
+                            }
                         }
                     } // End of "is neighbour"
                 } // End of "is processor"
@@ -1087,24 +1156,20 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
             const labelList& curNeighbour = curMesh.faceNeighbour();
             fpAddr.setSize(curFaces.size());
 
-            // Collect globally shared point labels
+            // Collect newly marked globally shared point labels
-            const labelList& curSpl = sharedData.sharedPointLabels()[procI];
+            forAll (curSpAddr, splI)
-
-            forAll (curSpl, splI)
             {
-                // From other processors, check if point is already marked
+                if (globalPointMapping[curSpAddr[splI]] < 0)
-                // If not, mark it; otherwise compare (and correct?) with local
-                // mark
-                if (globalPointMapping[curSpl[splI]] < 0)
                 {
-                    globalPointMapping[curSpl[splI]] = ppAddr[curSpl[splI]];
+                    globalPointMapping[curSpAddr[splI]] = ppAddr[curSpl[splI]];
                 }
                 else
                 {
-                    // Compare.  Is this needed - should always be OK.
+                    // If point is already set, compare sync
+                    // to detect merge error.  Debug
                     if
                     (
-                        globalPointMapping[curSpl[splI]]
+                        globalPointMapping[curSpAddr[splI]]
                      != ppAddr[curSpl[splI]]
                     )
                     {
@@ -1113,7 +1178,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
                             "autoPtr<fvMesh> "
                             "processorMeshesReconstructor::"
                             "reconstructMesh(const Time& db)"
-                        )   << "Loss of sync???"
+                        )   << "Loss of global sync: "
+                            << globalPointMapping[curSpAddr[splI]] << " and "
+                            << ppAddr[curSpl[splI]] << nl
                             << abort(FatalError);
                     }
                 }
@@ -1542,20 +1609,6 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
 
     // TODO: point, face and cell zones
 
-    Info<< "Reconstructed addressing: " << nl;
-    forAll (meshes_, procI)
-    {
-        if (meshes_.set(procI))
-        {
-            Info<< "Proc " << procI
-                << " point addr: " << pointProcAddressing_[procI].size()
-                << " face addr: " << faceProcAddressing_[procI].size()
-                << " cell addr: " << cellProcAddressing_[procI].size()
-                << " boundary addr: " << boundaryProcAddressing_[procI].size()
-                << endl;
-        }
-    }
-
     return globalMeshPtr;
 }
 

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C

@@ -194,19 +194,18 @@ void Foam::sharedPoints::calcSharedPoints()
     labelListList patchPairs = procPatchPairs();
 
     // Communicate and count global points
-    labelList nGlobalPointsPerProc(meshes_.size(), 0);
 
     // Identify, count and communicate points across processor boundaries
     // Repeat until the number of points per processor stabilises,
     // ie. no further points are found through communication
-    label oldNTotalPoints, newNTotalPoints;
+
+    label nSynced;
+
     do
     {
-        oldNTotalPoints = sum(nGlobalPointsPerProc);
+        nSynced = 0;
-
-        // Reset the list
-        nGlobalPointsPerProc = 0;
 
+        // Sync mark across processor boundaries
         forAll (meshes_, meshI)
         {
             if (meshes_.set(meshI))
@@ -223,13 +222,16 @@ void Foam::sharedPoints::calcSharedPoints()
 
                     if (isA<processorPolyPatch>(pp))
                     {
-                        // Found processor patch
-                        const labelList& patchMeshPoints = pp.meshPoints();
-
                         // My processor patch
                         const processorPolyPatch& myProcPatch =
                             refCast<const processorPolyPatch>(pp);
 
+                        // Get local mesh points
+                        const labelList& patchMeshPoints = pp.meshPoints();
+
+                        // Get my local faces
+                        const faceList& patchLocalFaces = pp.localFaces();
+
                         // Get neighbour processor ID
                         const int nbrProcID = myProcPatch.neighbProcNo();
 
@@ -238,65 +240,115 @@ void Foam::sharedPoints::calcSharedPoints()
                             meshes_[nbrProcID].boundaryMesh()
                             [patchPairs[meshI][patchI]];
 
+                        // Get neighbour mesh points
                         const labelList& nbrMeshPoints = nbrPatch.meshPoints();
                         
-                        forAll (patchMeshPoints, mpI)
+                        // Get neighbour local faces
+                        const faceList& nbrLocalFaces = nbrPatch.localFaces();
+                        
+                        labelList& nbrMarkedPoints = markedPoints[nbrProcID];
+
+                        // Note:
+                        // Cannot loop over mesh points because they are sorted
+                        // Use face loops as they are synchronised
+                        // HJ, 21/May/2018
+                        forAll (patchLocalFaces, faceI)
                         {
-                            if (curMarkedPoints[patchMeshPoints[mpI]] > 1)
+                            const face& curFace = patchLocalFaces[faceI];
+
+                            // Reverse neighbour face (copy)
+                            const face nbrFace =
+                                nbrLocalFaces[faceI].reverseFace();
+                            
+                            forAll (curFace, fpI)
                             {
-                                // Mark the point on the other processor/side
+                                const label patchMpI =
-                                markedPoints[nbrProcID][nbrMeshPoints[mpI]] =
+                                    patchMeshPoints[curFace[fpI]];
-                                    curMarkedPoints[patchMeshPoints[mpI]];
-                            }
-                        }
-                    }
-                }
 
-                // Count number of shared points per processor
+                                const label nbrMpI =
-                forAll (curMarkedPoints, cpI)
+                                    nbrMeshPoints[nbrFace[fpI]];
+
+                                if
+                                (
+                                    curMarkedPoints[patchMpI] > 1
+                                 && nbrMarkedPoints[nbrMpI] <= 1
+                                )
                                 {
-                    if (curMarkedPoints[cpI] > 1)
+                                    nbrMarkedPoints[nbrMpI] =
+                                        curMarkedPoints[patchMpI];
+                                         
+                                    nSynced++;
+                                }
+
+                                if
+                                (
+                                    curMarkedPoints[patchMpI] <= 1
+                                 && nbrMarkedPoints[nbrMpI] > 1
+                                )
                                 {
-                        nGlobalPointsPerProc[meshI]++;
+                                    curMarkedPoints[patchMpI] =
+                                        nbrMarkedPoints[nbrMpI];
+
+                                    nSynced++;
+                                }
+                            }
+                        }
                     }
                 }
             }
         }
 
-        newNTotalPoints = sum(nGlobalPointsPerProc);
+        Info<< "Proc merge pass.  nSynced = " << nSynced << endl;
+    } while (nSynced > 0);
 
-        Info<< "Proc merge pass: " << oldNTotalPoints << " "
+    // Grab marked points
-            << newNTotalPoints << endl;
-    } while (oldNTotalPoints != newNTotalPoints);
-
-    Info<< "Number of shared points per processor: " << nGlobalPointsPerProc
-        << endl;
-
-    // Collect points for every processor, in order to re-use the markedPoints
-    // list.  Note: the list of global labels of shared points
-    // will be collected later
     forAll (meshes_, meshI)
     {
         if (meshes_.set(meshI))
         {
-            labelList& curSharedPoints = sharedPointLabels_[meshI];
+            // Count the points and then collect
-
-            curSharedPoints.setSize(nGlobalPointsPerProc[meshI]);
-
-            // Count inserted points
             label nShared = 0;
 
-            // Get point marking
             const labelList& curMarkedPoints = markedPoints[meshI];
 
             forAll (curMarkedPoints, pointI)
             {
                 if (curMarkedPoints[pointI] > 1)
                 {
-                    curSharedPoints[nShared] = pointI;
                     nShared++;
                 }
             }
+
+            labelList& curShared = sharedPointLabels_[meshI];
+            curShared.setSize(nShared);
+
+            // Re-use the counter for insertion
+            nShared = 0;
+
+            forAll (curMarkedPoints, pointI)
+            {
+                if (curMarkedPoints[pointI] > 1)
+                {
+                    curShared[nShared] = pointI;
+                    nShared++;
+                }
+            }
+
+            const pointField& P = meshes_[meshI].points();
+            OFstream ppp("points" + name(meshI) + ".vtk");
+            ppp << "# vtk DataFile Version 2.0" << nl
+                << "points" + name(meshI) << ".vtk" << nl
+                << "ASCII" << nl
+                << "DATASET POLYDATA" << nl
+                << "POINTS " << curShared.size() << " float" << nl;
+
+            forAll (curShared, i)
+            {
+                ppp << float(P[curShared[i]].x()) << ' '
+                    << float(P[curShared[i]].y()) << ' '
+                    << float(P[curShared[i]].z())
+                    << nl;
+            }
         }
     }
 
@@ -324,6 +376,7 @@ void Foam::sharedPoints::calcSharedPoints()
 
             labelList& curMarkedPoints = markedPoints[meshI];
 
+            // Collect shared point addressing
             forAll (curSharedPoints, spI)
             {
                 if (curMarkedPoints[curSharedPoints[spI]] == -1)
@@ -356,14 +409,19 @@ void Foam::sharedPoints::calcSharedPoints()
 
                 if (isA<processorPolyPatch>(pp))
                 {
-                    // Found processor patch
-
                     // My processor patch
                     const processorPolyPatch& myProcPatch =
                         refCast<const processorPolyPatch>(pp);
 
+                    if (myProcPatch.master())
+                    {
+                        // Get local mesh points
                         const labelList& patchMeshPoints = pp.meshPoints();
 
+                        // Get my local faces
+                        const faceList& patchLocalFaces = pp.localFaces();
+
+
                         // Get neighbour processor ID
                         const int nbrProcID = myProcPatch.neighbProcNo();
 
@@ -372,15 +430,39 @@ void Foam::sharedPoints::calcSharedPoints()
                             meshes_[nbrProcID].boundaryMesh()
                             [patchPairs[meshI][patchI]];
 
+                        // Get neighbour mesh points
                         const labelList& nbrMeshPoints = nbrPatch.meshPoints();
 
-                    forAll (patchMeshPoints, mpI)
+                        // Get neighbour local faces
+                        const faceList& nbrLocalFaces = nbrPatch.localFaces();
+                        
+                        // Note:
+                        // Cannot loop over mesh points because they are sorted
+                        // Use face loops as they are synchronised
+                        // HJ, 21/May/2018
+                        forAll (patchLocalFaces, faceI)
                         {
-                        if (curMarkedPoints[patchMeshPoints[mpI]] > -1)
+                            const face& curFace = patchLocalFaces[faceI];
+
+                            // Reverse neighbour face (copy)
+                            const face nbrFace =
+                                nbrLocalFaces[faceI].reverseFace();
+
+                            forAll (curFace, fpI)
+                            {
+                                const label patchMpI =
+                                    patchMeshPoints[curFace[fpI]];
+
+                                const label nbrMpI =
+                                    nbrMeshPoints[nbrFace[fpI]];
+
+                                if (curMarkedPoints[patchMpI] > -1)
                                 {
                                     // Mark opposite side
-                            markedPoints[nbrProcID][nbrMeshPoints[mpI]] =
+                                    markedPoints[nbrProcID][nbrMpI] =
-                                curMarkedPoints[patchMeshPoints[mpI]];
+                                        curMarkedPoints[patchMpI];
+                                }
+                            }
                         }
                     }
                 }

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H

@@ -35,18 +35,12 @@ Description
       - patch face ordering to be ok
       - f[0] ordering on patch faces to be ok.
 
-    Works by constructing equivalence lists for all the points on processor
+    The class will create a unique list of multiply shared points
-    patches. These list are procPointList and give processor and meshPoint
+    (on more than 2 mesh pieces).  On each mesh piece, it will provide a
-    label on that processor.
+    a list of local point labels that are globally shared (sharedPointLabels)
-    E.g.
+    and for each of such points the index in the list of globally shared
-    @verbatim
+    points (sharedPointAddr).  The class also provides the number of globally
-          ((7 93)(4 731)(3 114))
+    shared points across all mesh pieces (nGlobalPoints)
-    @endverbatim
-
-    means point 93 on proc7 is connected to point 731 on proc4 and 114 on proc3.
-    It then gets the lowest numbered processor (the 'master') to request a
-    sharedPoint label from processor0 and it redistributes this label back to
-    the other processors in the equivalence list.
 
 Note
     Currently operating with a PtrList<fvMesh>, whereas the operation actually
@@ -118,9 +112,6 @@ class sharedPoints
 
 public:
 
-    // Static data members
-
-
     // Constructors
 
         //- Construct from the list of meshes