Frontal marking SMP implementation

git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@28 fdcce57e-7e00-11e2-b579-49867b4cea03

meshLibrary/utilities/helperFunctions/helperFunctionsFrontalMarking.C

@@ -109,41 +109,212 @@ label groupMarking
     elementInGroup.setSize(neighbourCalculator.size());
     elementInGroup = -1;
 
-    forAll(neighbourCalculator, elI)
+    VRWGraph neighbouringGroups;
+
+    label nThreads(1);
+
+    # ifdef USE_OMP
+    nThreads = 3 * omp_get_num_procs();
+    # endif
+
+    # ifdef USE_OMP
+    # pragma omp parallel if( neighbourCalculator.size() > 1000 ) \
+    num_threads(nThreads)
+    # endif
     {
-        if( elementInGroup[elI] != -1 )
-            continue;
-        if( !selector(elI) )
-            continue;
+        const label chunkSize = neighbourCalculator.size() / nThreads + 1;
+
+        const label threadI = omp_get_thread_num();
+
+        LongList<std::pair<label, label> > threadCommPairs;
+
+        const label minEl = threadI * chunkSize;
+
+        const label maxEl =
+            Foam::min
+            (
+                neighbourCalculator.size(),
+                minEl + chunkSize
+            );
+
+        for(label elI=minEl;elI<maxEl;++elI)
+        {
+            if( elementInGroup[elI] != -1 )
+                continue;
+            if( !selector(elI) )
+                continue;
+
+            label groupI;
+            # ifdef USE_OMP
+            # pragma omp critical
+            # endif
+            groupI = nGroups++;
+
+            elementInGroup[elI] = groupI;
+            labelListPMG front;
+            front.append(elI);
 
-        elementInGroup[elI] = nGroups;
-        labelListPMG front;
-        front.append(elI);
+            while( front.size() )
+            {
+                const label eLabel = front.removeLastElement();
+
+                DynList<label> neighbours;
+                neighbourCalculator(eLabel, neighbours);
+
+                forAll(neighbours, neiI)
+                {
+                    const label nei = neighbours[neiI];
+
+                    if( (nei < 0) || (elementInGroup[nei] != -1) )
+                        continue;
+
+                    if( (nei < minEl) || (nei >= maxEl) )
+                    {
+                        //- this is a communication interface between
+                        //- two threads
+                        threadCommPairs.append(std::make_pair(elI, nei));
+                    }
+                    else if( selector(nei) )
+                    {
+                        //- this element is part of the same group
+                        elementInGroup[nei] = groupI;
+                        front.append(nei);
+                    }
+                }
+            }
+        }
+
+        # ifdef USE_OMP
+        # pragma omp barrier
+        # endif
 
-        while( front.size() )
+        //- find group to neighbouring groups addressing
+        List<DynList<label> > localNeiGroups(nGroups);
+        forAll(threadCommPairs, cfI)
         {
-            const label eLabel = front.removeLastElement();
+            const std::pair<label, label>& lp = threadCommPairs[cfI];
+            const label groupI = elementInGroup[lp.first];
+            const label neiGroup = elementInGroup[lp.second];
 
-            DynList<label> neighbours;
-            neighbourCalculator(eLabel, neighbours);
+            if( (neiGroup >= nGroups) || (groupI >= nGroups) )
+                FatalError << "neiGroup " << neiGroup
+                    << " groupI " << groupI << " are >= than "
+                    << "nGroups " << nGroups << abort(FatalError);
 
-            forAll(neighbours, neiI)
+            if( neiGroup != -1 )
             {
-                const label nei = neighbours[neiI];
+                localNeiGroups[groupI].appendIfNotIn(neiGroup);
+                localNeiGroups[neiGroup].appendIfNotIn(groupI);
+            }
+        }
+
+        # ifdef USE_OMP
+        # pragma omp critical
+        # endif
+        {
+            neighbouringGroups.setSize(nGroups);
+
+            forAll(localNeiGroups, groupI)
+            {
+                const DynList<label>& lGroups = localNeiGroups[groupI];
+
+                neighbouringGroups.appendIfNotIn(groupI, groupI);
 
-                if( (nei < 0) || (elementInGroup[nei] != -1) )
+                forAll(lGroups, i)
+                    neighbouringGroups.appendIfNotIn(groupI, lGroups[i]);
+            }
+        }
+    }
+
+    DynList<label> globalGroupLabel;
+    globalGroupLabel.setSize(nGroups);
+
+    bool changed;
+
+    do
+    {
+        changed = false;
+
+        DynList<label> newGroupLabel;
+        newGroupLabel.setSize(globalGroupLabel.size());
+        newGroupLabel = -1;
+
+        //- reduce the information about the groups
+        forAllReverse(neighbouringGroups, groupI)
+        {
+            forAllRow(neighbouringGroups, groupI, ngI)
+            {
+                const label neiGroup = neighbouringGroups(groupI, ngI);
+
+                if( neiGroup >= groupI )
                     continue;
 
-                if( selector(nei) )
+                forAllRow(neighbouringGroups, groupI, i)
                 {
-                    elementInGroup[nei] = nGroups;
-                    front.append(nei);
+                    const label grI = neighbouringGroups(groupI, i);
+                    neighbouringGroups.appendIfNotIn(neiGroup, grI);
                 }
             }
         }
 
-        //- increment the number of groups
-        ++nGroups;
+        label counter(0);
+        forAll(neighbouringGroups, groupI)
+        {
+            if( newGroupLabel[groupI] != -1 )
+                continue;
+
+            forAllRow(neighbouringGroups, groupI, ngI)
+                newGroupLabel[neighbouringGroups(groupI, ngI)] = counter;
+
+            ++counter;
+        }
+
+        if( counter < nGroups )
+        {
+            changed = true;
+            nGroups = counter;
+        }
+
+        globalGroupLabel = newGroupLabel;
+
+        if( Pstream::parRun() )
+        {
+            //- reduce the groups over processors
+            labelList nGroupsAtProc(Pstream::nProcs());
+            nGroupsAtProc[Pstream::myProcNo()] = nGroups;
+
+            Pstream::gatherList(nGroupsAtProc);
+            Pstream::scatterList(nGroupsAtProc);
+
+            label startGroup(0);
+            for(label procI=0;procI<Pstream::myProcNo();++procI)
+                startGroup += nGroupsAtProc[procI];
+
+            //- find the neighbouring groups
+            std::map<label, DynList<label> > neiGroups;
+
+            //- collect groups on other processors
+            //- this operator implements the algorithm for exchanging data
+            //- over processors and collects information which groups
+            //- are connected ovr inter-processor boundaries
+            neighbourCalculator.collectGroups
+            (
+                neiGroups,
+                elementInGroup,
+                newGroupLabel,
+                startGroup
+            );
+        }
+
+    } while( changed );
+
+    //- set the global group label
+    forAll(elementInGroup, elI)
+    {
+        if( elementInGroup[elI] < 0 )
+            continue;
+
+        elementInGroup[elI] = globalGroupLabel[elementInGroup[elI]];
     }
 
     return nGroups;

testingInterfaces/testFrontalMarking/Make/options

-EXE_INC = -I$(LIB_SRC)/triSurface/lnInclude -I$(LIB_SRC)/meshTools/lnInclude -I../../meshLibrary/lnInclude
+EXE_INC = -DUSE_OMP -I$(LIB_SRC)/triSurface/lnInclude -I$(LIB_SRC)/meshTools/lnInclude -I../../meshLibrary/lnInclude
 EXE_LIBS = -lMeshLibrary -ltriSurface -lmeshTools

testingInterfaces/testFrontalMarking/testFrontalMarking.C

@@ -126,6 +126,83 @@ public:
                 neighbourCells.append(nei);
         }
     }
+
+    template<class labelListType>
+    void collectGroups
+    (
+        std::map<label, DynList<label> >& neiGroups,
+        const labelListType& elementInGroup,
+        const DynList<label>& localGroupLabel,
+        const label startGroupAtProc
+    ) const
+    {
+        const PtrList<writeProcessorPatch>& procBoundaries =
+            mesh_.procBoundaries();
+        const labelList& owner = mesh_.owner();
+
+        //- send the data to other processors
+        forAll(procBoundaries, patchI)
+        {
+            const label start = procBoundaries[patchI].patchStart();
+            const label size = procBoundaries[patchI].patchSize();
+
+            labelList groupOwner(procBoundaries[patchI].patchSize());
+            for(label faceI=0;faceI<size;++faceI)
+            {
+                const label groupI = elementInGroup[owner[start+faceI]];
+
+                if( groupI < 0 )
+                {
+                    groupOwner[faceI] = groupI;
+                    continue;
+                }
+
+                groupOwner[faceI] = startGroupAtProc + localGroupLabel[groupI];
+            }
+
+            OPstream toOtherProc
+            (
+                Pstream::blocking,
+                procBoundaries[patchI].neiProcNo(),
+                groupOwner.byteSize()
+            );
+
+            toOtherProc << groupOwner;
+        }
+
+        //- receive data from other processors
+        forAll(procBoundaries, patchI)
+        {
+            const label start = procBoundaries[patchI].patchStart();
+
+            labelList receivedData;
+
+            IPstream fromOtherProc
+            (
+                Pstream::blocking,
+                procBoundaries[patchI].neiProcNo()
+            );
+
+            fromOtherProc >> receivedData;
+
+            forAll(receivedData, faceI)
+            {
+                if( receivedData[faceI] < 0 )
+                    continue;
+
+                const label groupI = elementInGroup[owner[start+faceI]];
+
+                if( groupI < 0 )
+                    continue;
+
+                DynList<label>& ng =
+                    neiGroups[startGroupAtProc + localGroupLabel[groupI]];
+
+                //- store the connection fo the inter-processor boundary
+                ng.appendIfNotIn(receivedData[faceI]);
+            }
+        }
+    }
 };
 
 class meshSelectorOperator
@@ -187,7 +264,8 @@ int main(int argc, char *argv[])
 
     help::frontalMarking(result, start, on, oso);
 
-    Info << "Number of octree elements meeting the criteria " << result << endl;
+    Info << "Number of octree elements meeting the criteria "
+         << result.size() << endl;
 
     Info<< "Execution time for octree creation = "
         << runTime.elapsedCpuTime()