Optimisation of meshOctreeInsideOutside + bug fix when reading 0 length

LongList in binary format


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

meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C

@@ -379,8 +379,6 @@ cartesianMeshGenerator::cartesianMeshGenerator(const Time& time)
 
     surfacePtr_ = new triSurf(db_.path()/surfaceFile);
 
-    surfacePtr_->writeSurface("geom.fms");
-
     if( surfacePtr_->featureEdges().size() != 0 )
     {
         //- create surface patches based on the feature edges

meshLibrary/utilities/containers/LongList/LongList.C

@@ -160,48 +160,47 @@ Foam::Istream& Foam::operator>>
         // Set list length to that read
         DL.setSize(size);
 
-        if( size == 0 )
+        // Read list contents depending on data format
+        if( (is.format() == IOstream::ASCII) || !contiguous<T>() )
         {
+            // Read beginning of contents
             char listDelimiter = is.readBeginList("List");
-            if( listDelimiter != token::BEGIN_LIST )
-            {
-                WarningIn
-                (
-                    "template<class T, Foam::label Offset>"
 
-                    "Foam::Istream& Foam::operator>>"
-                    "("
-                        "Foam::Istream& ,"
-                        "Foam::LongList<T, Offset>& DL"
-                    ")"
-                ) << "Missing ( after 0" << endl;
+            if( size == 0 )
+            {
+                if( listDelimiter != token::BEGIN_LIST )
+                {
+                    WarningIn
+                    (
+                        "template<class T, Foam::label Offset>"
 
-                return is;
-            }
+                        "Foam::Istream& Foam::operator>>"
+                        "("
+                            "Foam::Istream& ,"
+                            "Foam::LongList<T, Offset>& DL"
+                        ")"
+                    ) << "Missing ( after 0" << endl;
 
-            listDelimiter = is.readEndList("List");
-            if( listDelimiter != token::END_LIST )
-            {
-                WarningIn
-                (
-                    "template<class T, Foam::label Offset>"
-
-                    "Foam::Istream& Foam::operator>>"
-                    "("
-                        "Foam::Istream& ,"
-                        "Foam::LongList<T, Offset>& DL"
-                    ")"
-                ) << "Missing ) after 0(" << endl;
-            }
+                    return is;
+                }
 
-            return is;
-        }
+                listDelimiter = is.readEndList("List");
+                if( listDelimiter != token::END_LIST )
+                {
+                    WarningIn
+                    (
+                        "template<class T, Foam::label Offset>"
+
+                        "Foam::Istream& Foam::operator>>"
+                        "("
+                            "Foam::Istream& ,"
+                            "Foam::LongList<T, Offset>& DL"
+                        ")"
+                    ) << "Missing ) after 0(" << endl;
+                }
 
-        // Read list contents depending on data format
-        if( (is.format() == IOstream::ASCII) || !contiguous<T>() )
-        {
-            // Read beginning of contents
-            char listDelimiter = is.readBeginList("List");
+                return is;
+            }
 
             if( listDelimiter == token::BEGIN_LIST )
             {

meshLibrary/utilities/octrees/meshOctree/meshOctreeInsideOutside/meshOctreeInsideOutside.C

@@ -139,7 +139,7 @@ void meshOctreeInsideOutside::frontalMarking()
     labelListPMG frontCubes;
     DynList<label> neighbours(24);
 
-    label nGroup(0), chunkI(0), nChunks, chunkSize;
+    label nGroup(0), nThreads(1);
 
     const LongList<meshOctreeCube*>& leaves = octreeModifier_.leavesAccess();
     const meshOctree& octree = octreeModifier_.octree();
@@ -147,150 +147,143 @@ void meshOctreeInsideOutside::frontalMarking()
     boolList commCubes(leaves.size(), false);
 
     # ifdef USE_OMP
-    # pragma omp parallel if( leaves.size() > 1000 ) \
+    if( leaves.size() > 1000 )
+        nThreads = 3 * omp_get_num_procs();
+
+    # pragma omp parallel num_threads(nThreads) \
     private(frontCubes, neighbours)
     # endif
     {
         LongList<std::pair<label, label> > threadCommPairs;
 
-        # ifdef USE_OMP
-        # pragma omp master
+        const label threadI = omp_get_thread_num();
 
-        {
-            nChunks = 3 * omp_get_num_threads();
-            chunkSize = leaves.size() / nChunks + 1;
-        }
+        const label chunkSize = leaves.size() / nThreads + 1;
 
-        # pragma omp barrier
-        # else
-        nChunks = 1;
-        chunkSize = leaves.size();
-        # endif
+        const label minLeaf = threadI * chunkSize;
 
-        while( chunkI < nChunks )
+        const label maxLeaf = Foam::min(leaves.size(), minLeaf + chunkSize);
+
+        for(label leafI=minLeaf;leafI<maxLeaf;++leafI)
         {
-            label minLeaf, maxLeaf;
+            if( leaves[leafI]->hasContainedElements() )
+                continue;
+            if( cubeGroup_[leafI] != -1 )
+                continue;
+
+            label groupI;
             # ifdef USE_OMP
             # pragma omp critical
             # endif
-            minLeaf = chunkI++ * chunkSize;
+            groupI = nGroup++;
 
-            if( minLeaf >= leaves.size() )
-                break;
+            direction cType(meshOctreeCubeBasic::UNKNOWN);
+            frontCubes.clear();
+            frontCubes.append(leafI);
+            cubeGroup_[leafI] = groupI;
 
-            maxLeaf = Foam::min(leaves.size(), minLeaf + chunkSize);
+            labelListPMG neiDATACubes;
 
-            for(label leafI=minLeaf;leafI<maxLeaf;++leafI)
+            while( frontCubes.size() )
             {
-                if( leaves[leafI]->hasContainedElements() )
-                    continue;
-                if( cubeGroup_[leafI] != -1 )
-                    continue;
-
-                label groupI;
-                # ifdef USE_OMP
-                # pragma omp critical
-                # endif
-                groupI = nGroup++;
+                const label fLabel = frontCubes.removeLastElement();
+                octree.findNeighboursForLeaf(fLabel, neighbours);
 
-                direction cType(meshOctreeCubeBasic::UNKNOWN);
-                frontCubes.clear();
-                frontCubes.append(leafI);
-                cubeGroup_[leafI] = groupI;
-
-                labelListPMG neiDATACubes;
-
-                while( frontCubes.size() )
+                forAll(neighbours, neiI)
                 {
-                    const label fLabel = frontCubes.removeLastElement();
-                    octree.findNeighboursForLeaf(fLabel, neighbours);
-
-                    forAll(neighbours, neiI)
+                    const label nei = neighbours[neiI];
+                    if( (nei >= minLeaf) && (nei < maxLeaf) )
                     {
-                        const label nei = neighbours[neiI];
-                        if( (nei >= minLeaf) && (nei < maxLeaf) )
-                        {
-                            if( cubeGroup_[nei] != -1 )
-                                continue;
+                        if( cubeGroup_[nei] != -1 )
+                            continue;
 
-                            if( leaves[nei]->hasContainedElements() )
-                            {
-                                neiDATACubes.append(nei);
-                            }
-                            else
-                            {
-                                frontCubes.append(nei);
-                                cubeGroup_[nei] = groupI;
-                            }
+                        if( leaves[nei]->hasContainedElements() )
+                        {
+                            neiDATACubes.append(nei);
                         }
-                        else if( nei == -1 )
+                        else
                         {
-                            cType = meshOctreeCubeBasic::OUTSIDE;
+                            frontCubes.append(nei);
+                            cubeGroup_[nei] = groupI;
                         }
-                        else if( nei == meshOctreeCubeBasic::OTHERPROC )
+                    }
+                    else if( nei == -1 )
+                    {
+                        cType = meshOctreeCubeBasic::OUTSIDE;
+                    }
+                    else if( nei == meshOctreeCubeBasic::OTHERPROC )
+                    {
+                        commCubes[fLabel] = true;
+                    }
+                    else
+                    {
+                        if( leaves[nei]->hasContainedElements() )
                         {
-                            commCubes[fLabel] = true;
+                            neiDATACubes.append(nei);
                         }
                         else
                         {
-                            if( leaves[nei]->hasContainedElements() )
-                            {
-                                neiDATACubes.append(nei);
-                            }
-                            else
-                            {
-                                threadCommPairs.append
-                                (
-                                    std::make_pair(fLabel, nei)
-                                );
-                            }
+                            threadCommPairs.append
+                            (
+                                std::make_pair(fLabel, nei)
+                            );
                         }
                     }
                 }
+            }
 
-                # ifdef USE_OMP
-                # pragma omp critical
-                # endif
-                {
-                    if( groupI >= boundaryDATACubes_.size() )
-                        boundaryDATACubes_.setSize(groupI+1);
+            # ifdef USE_OMP
+            # pragma omp critical
+            # endif
+            {
+                if( groupI >= boundaryDATACubes_.size() )
+                    boundaryDATACubes_.setSize(groupI+1);
 
-                    boundaryDATACubes_.setRow(groupI, neiDATACubes);
-                    groupType_[groupI] = cType;
-                }
+                boundaryDATACubes_.setRow(groupI, neiDATACubes);
+                groupType_[groupI] = cType;
             }
         }
 
-        # ifdef USE_OMP
-        # pragma omp barrier
-
-        # pragma omp master
-        # endif
-        neighbouringGroups_.setSize(nGroup);
-
         # ifdef USE_OMP
         # pragma omp barrier
         # endif
 
-        forAll(threadCommPairs, pairI)
+        //- find group to neighbouring groups addressing
+        List<DynList<label> > localNeiGroups(nGroup);
+        forAll(threadCommPairs, cfI)
         {
-            const std::pair<label, label>& cubesPair = threadCommPairs[pairI];
-            const label groupI = cubeGroup_[cubesPair.first];
-            const label neiGroup = cubeGroup_[cubesPair.second];
+            const std::pair<label, label>& lp = threadCommPairs[cfI];
+            const label groupI = cubeGroup_[lp.first];
+            const label neiGroup = cubeGroup_[lp.second];
 
-            if( neiGroup >= nGroup )
-                FatalError << "neiGroup " << neiGroup << " is >= than "
+            if( (neiGroup >= nGroup) || (groupI >= nGroup) )
+                FatalError << "neiGroup " << neiGroup
+                    << " groupI " << groupI << " are >= than "
                     << "nGroups " << nGroup << abort(FatalError);
 
-            if(
-                (neiGroup != -1) &&
-                !neighbouringGroups_.contains(groupI, neiGroup)
-            )
+            if( neiGroup != -1 )
             {
-                # ifdef USE_OMP
-                # pragma omp critical
-                # endif
-                neighbouringGroups_.append(groupI, neiGroup);
+                localNeiGroups[groupI].appendIfNotIn(neiGroup);
+                localNeiGroups[neiGroup].appendIfNotIn(groupI);
+            }
+        }
+
+        # ifdef USE_OMP
+        # pragma omp critical
+        # endif
+        {
+            neighbouringGroups_.setSize(nGroup);
+
+            forAll(localNeiGroups, groupI)
+            {
+                const DynList<label>& lGroups = localNeiGroups[groupI];
+
+                neighbouringGroups_.appendIfNotIn(groupI, groupI);
+
+                forAll(lGroups, i)
+                {
+                    neighbouringGroups_.append(groupI, lGroups[i]);
+                }
             }
         }
     }