From 5f840ba74fd4bc5e085c8c985199033e155369a7 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 21 May 2018 13:00:34 +0100 Subject: [PATCH] Fixed bug in shared points for reconstruction --- .../finiteVolume/processorMeshesRebuild.C | 127 +++++++---- .../finiteVolume/sharedPoints.C | 202 ++++++++++++------ .../finiteVolume/sharedPoints.H | 21 +- 3 files changed, 238 insertions(+), 112 deletions(-) diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C index a5ed04ee3..3559ab3f9 100644 --- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C +++ b/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 @@ -1054,9 +1095,37 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db) forAll (curRF, pointI) { - // Mapping is established - ppAddr[curRF[pointI]] = - masterPpAddr[curMF[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 " + "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" @@ -1087,24 +1156,20 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db) const labelList& curNeighbour = curMesh.faceNeighbour(); fpAddr.setSize(curFaces.size()); - // Collect globally shared point labels - const labelList& curSpl = sharedData.sharedPointLabels()[procI]; - - forAll (curSpl, splI) + // Collect newly marked globally shared point labels + forAll (curSpAddr, splI) { - // From other processors, check if point is already marked - // If not, mark it; otherwise compare (and correct?) with local - // mark - if (globalPointMapping[curSpl[splI]] < 0) + if (globalPointMapping[curSpAddr[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 " "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; } diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C index af30fb647..68e223a4e 100644 --- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C +++ b/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); - - // Reset the list - nGlobalPointsPerProc = 0; + nSynced = 0; + // Sync mark across processor boundaries forAll (meshes_, meshI) { if (meshes_.set(meshI)) @@ -223,13 +222,16 @@ void Foam::sharedPoints::calcSharedPoints() if (isA(pp)) { - // Found processor patch - const labelList& patchMeshPoints = pp.meshPoints(); - // My processor patch const processorPolyPatch& myProcPatch = refCast(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(); + + // Get neighbour local faces + const faceList& nbrLocalFaces = nbrPatch.localFaces(); + + labelList& nbrMarkedPoints = markedPoints[nbrProcID]; - forAll (patchMeshPoints, mpI) + // 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 - markedPoints[nbrProcID][nbrMeshPoints[mpI]] = - curMarkedPoints[patchMeshPoints[mpI]]; + const label patchMpI = + patchMeshPoints[curFace[fpI]]; + + const label nbrMpI = + nbrMeshPoints[nbrFace[fpI]]; + + if + ( + curMarkedPoints[patchMpI] > 1 + && nbrMarkedPoints[nbrMpI] <= 1 + ) + { + nbrMarkedPoints[nbrMpI] = + curMarkedPoints[patchMpI]; + + nSynced++; + } + + if + ( + curMarkedPoints[patchMpI] <= 1 + && nbrMarkedPoints[nbrMpI] > 1 + ) + { + curMarkedPoints[patchMpI] = + nbrMarkedPoints[nbrMpI]; + + nSynced++; + } } } } } - - // Count number of shared points per processor - forAll (curMarkedPoints, cpI) - { - if (curMarkedPoints[cpI] > 1) - { - nGlobalPointsPerProc[meshI]++; - } - } } } - newNTotalPoints = sum(nGlobalPointsPerProc); + Info<< "Proc merge pass. nSynced = " << nSynced << endl; + } while (nSynced > 0); - Info<< "Proc merge pass: " << oldNTotalPoints << " " - << 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 + // Grab marked points forAll (meshes_, meshI) { if (meshes_.set(meshI)) { - labelList& curSharedPoints = sharedPointLabels_[meshI]; - - curSharedPoints.setSize(nGlobalPointsPerProc[meshI]); - - // Count inserted points + // Count the points and then collect 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,31 +409,60 @@ void Foam::sharedPoints::calcSharedPoints() if (isA(pp)) { - // Found processor patch - // My processor patch const processorPolyPatch& myProcPatch = refCast(pp); - const labelList& patchMeshPoints = pp.meshPoints(); - - // Get neighbour processor ID - const int nbrProcID = myProcPatch.neighbProcNo(); - - // Neighbour patch - const polyPatch& nbrPatch = - meshes_[nbrProcID].boundaryMesh() - [patchPairs[meshI][patchI]]; - - const labelList& nbrMeshPoints = nbrPatch.meshPoints(); - - forAll (patchMeshPoints, mpI) + if (myProcPatch.master()) { - if (curMarkedPoints[patchMeshPoints[mpI]] > -1) + // 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(); + + // Neighbour patch + const polyPatch& nbrPatch = + meshes_[nbrProcID].boundaryMesh() + [patchPairs[meshI][patchI]]; + + // Get neighbour mesh points + const labelList& nbrMeshPoints = nbrPatch.meshPoints(); + + // 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) { - // Mark opposite side - markedPoints[nbrProcID][nbrMeshPoints[mpI]] = - curMarkedPoints[patchMeshPoints[mpI]]; + 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][nbrMpI] = + curMarkedPoints[patchMpI]; + } + } } } } diff --git a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H b/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H index 16c595da6..b39109b32 100644 --- a/src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H +++ b/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 - patches. These list are procPointList and give processor and meshPoint - label on that processor. - E.g. - @verbatim - ((7 93)(4 731)(3 114)) - @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. + The class will create a unique list of multiply shared points + (on more than 2 mesh pieces). On each mesh piece, it will provide a + a list of local point labels that are globally shared (sharedPointLabels) + and for each of such points the index in the list of globally shared + points (sharedPointAddr). The class also provides the number of globally + shared points across all mesh pieces (nGlobalPoints) Note Currently operating with a PtrList, whereas the operation actually @@ -118,9 +112,6 @@ class sharedPoints public: - // Static data members - - // Constructors //- Construct from the list of meshes