Bug fix in boundary laye refinement + disallow calculation of

polyMeshGenAddressing inside a parallel region.


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

meshLibrary/utilities/boundaryLayers/refineBoundaryLayersFunctions.C

@@ -276,6 +276,9 @@ public:
             return true;
         }
 
+        if( (nQuads + 2) != c.size() )
+            return false;
+
         if( nBndFaces != 1 )
             return false;
 
@@ -288,7 +291,7 @@ public:
             const bool sEdge =
                 help::shareAnEdge(faces[c[baseFace]], faces[c[fI]]);
 
-            if( faces[c[fI]].size() == 4 && sEdge )
+            if( (faces[c[fI]].size() == 4) && sEdge )
             {
                 ++nQuadsAttachedToBaseFace;
             }
@@ -328,6 +331,8 @@ void refineBoundaryLayers::analyseLayers()
     mse.faceOwners();
     mse.faceEdges();
     mse.edgeFaces();
+    mse.edges();
+    mse.boundaryPointEdges();
 
     //- find layers in patch
     labelLongList bndFaceInLayer;
@@ -380,6 +385,10 @@ void refineBoundaryLayers::analyseLayers()
         }
     }
 
+    # ifdef DEBUGLayer
+    Info << "Layer at patch " << layerAtPatch_ << endl;
+    # endif
+
     //- set the information which patches are a single boundary layer face
     patchesInLayer_.setSize(nValidLayers);
     forAll(layerAtPatch_, patchI)
@@ -394,6 +403,8 @@ void refineBoundaryLayers::analyseLayers()
     }
 
     # ifdef DEBUGLayer
+    Info << "Patches in layer " << patchesInLayer_ << endl;
+
     //- write layers to a subset
     std::map<label, label> layerId;
     for(label i=0;i<nValidLayers;++i)
@@ -417,18 +428,17 @@ void refineBoundaryLayers::analyseLayers()
     labelList nLayersAtPatch(layerAtPatch_.size(), -1);
     boolList protectedValue(layerAtPatch_.size(), false);
 
-    forAll(boundaries, patchI)
+    forAll(patchesInLayer_, layerI)
     {
-        const label layerI = layerAtPatch_[patchI];
-
-        if( layerI < 0 )
-            continue;
+        const DynList<word>& layerPatches = patchesInLayer_[layerI];
 
         label maxNumLayers(1);
+        bool hasLocalValue(false);
 
-        forAll(patchesInLayer_[layerI], lpI)
+        //- find the maximum requested number of layers over the layer
+        forAll(layerPatches, lpI)
         {
-            const word pName = patchesInLayer_[layerI][lpI];
+            const word pName = layerPatches[lpI];
 
             std::map<word, label>::const_iterator it =
                 numLayersForPatch_.find(pName);
@@ -441,28 +451,28 @@ void refineBoundaryLayers::analyseLayers()
                     discontinuousLayersForPatch_.end()
                 )
                 {
-                    //- set the numbe of layers and lock this location
+                    //- set the number of layers and lock this location
                     nLayersAtPatch[patchNameToIndex[pName]] = it->second;
                     protectedValue[patchNameToIndex[pName]] = true;
+                    hasLocalValue = true;
                 }
                 else
                 {
                     //- take the maximum number of layers
                     maxNumLayers = Foam::max(maxNumLayers, it->second);
+                    hasLocalValue = true;
                 }
             }
         }
 
-        if( maxNumLayers == 1 )
-        {
-            //- apply global settings to the patches which are not overriden
+        //- apply the global value if no local values exist
+        if( !hasLocalValue )
             maxNumLayers = globalNumLayers_;
-        }
 
-        //- set the number of layer to all patches which are not protected
-        forAll(patchesInLayer_[layerI], lpI)
+        //- apply the maximum number of ayer of all unprotected patches
+        forAll(layerPatches, lpI)
         {
-            const label ptchI = patchNameToIndex[patchesInLayer_[layerI][lpI]];
+            const label ptchI = patchNameToIndex[layerPatches[lpI]];
 
             if( !protectedValue[ptchI] )
                 nLayersAtPatch[ptchI] = maxNumLayers;

meshLibrary/utilities/helperFunctions/helperFunctionsFrontalMarking.C

@@ -119,7 +119,9 @@ public:
         //- operator for finding neighbours
         inline void operator()(const label groupI, DynList<label>& ng) const
         {
-            ng = neiGroups_[groupI];
+            ng.setSize(neiGroups_.sizeOfRow(groupI));
+            forAllRow(neiGroups_, groupI, i)
+                ng[i] = neiGroups_(groupI, i);
         }
 };
 
@@ -167,14 +169,17 @@ label groupMarking
     label nThreads(1);
 
     # ifdef USE_OMP
-    nThreads = 3 * omp_get_num_procs();
+    //nThreads = 3 * omp_get_num_procs();
     # endif
 
+    DynList<label> nGroupsAtThread(nThreads, 0);
+
     # ifdef USE_OMP
     # pragma omp parallel num_threads(nThreads)
     # endif
     {
-        const label chunkSize = neighbourCalculator.size() / nThreads + 1;
+        const label chunkSize =
+            Foam::max(1, neighbourCalculator.size() / nThreads);
 
         # ifdef USE_OMP
         const label threadI = omp_get_thread_num();
@@ -186,12 +191,12 @@ label groupMarking
 
         const label minEl = threadI * chunkSize;
 
-        const label maxEl =
-            Foam::min
-            (
-                neighbourCalculator.size(),
-                minEl + chunkSize
-            );
+        label maxEl = minEl + chunkSize;
+        if( threadI == (nThreads - 1) )
+            maxEl = Foam::max(maxEl, neighbourCalculator.size());
+
+        label& groupI = nGroupsAtThread[threadI];
+        groupI = 0;
 
         for(label elI=minEl;elI<maxEl;++elI)
         {
@@ -200,12 +205,6 @@ label groupMarking
             if( !selector(elI) )
                 continue;
 
-            label groupI;
-            # ifdef USE_OMP
-            # pragma omp critical
-            # endif
-            groupI = nGroups++;
-
             elementInGroup[elI] = groupI;
             labelLongList front;
             front.append(elI);
@@ -238,7 +237,28 @@ label groupMarking
                     }
                 }
             }
+
+            ++groupI;
+        }
+
+        # ifdef USE_OMP
+        # pragma omp barrier
+
+        # pragma omp master
+        {
+            forAll(nGroupsAtThread, i)
+                nGroups += nGroupsAtThread[i];
         }
+        # else
+        nGroups = groupI;
+        # endif
+
+        label startGroup(0);
+        for(label i=0;i<threadI;++i)
+            startGroup += nGroupsAtThread[i];
+
+        for(label elI=minEl;elI<maxEl;++elI)
+            elementInGroup[elI] += startGroup;
 
         # ifdef USE_OMP
         # pragma omp barrier
@@ -282,6 +302,17 @@ label groupMarking
         }
     }
 
+    forAll(neighbouringGroups, i)
+    {
+        labelList helper(neighbouringGroups.sizeOfRow(i));
+        forAllRow(neighbouringGroups, i, j)
+            helper[j] = neighbouringGroups(i, j);
+
+        sort(helper);
+
+        neighbouringGroups[i] = helper;
+    }
+
     //- start processing connections between the group and merge the connected
     //- ones into a new group
     DynList<label> globalGroupLabel;

meshLibrary/utilities/meshes/polyMeshGen/polyMeshGen.C

@@ -31,14 +31,13 @@ Description
 
 namespace Foam
 {
-    
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 polyMeshGen::polyMeshGen(const Time& t)
 :
     polyMeshGenCells(t)
-{
-}
+{}
 
 //- Construct from components without the boundary
 polyMeshGen::polyMeshGen
@@ -50,8 +49,7 @@ polyMeshGen::polyMeshGen
 )
 :
     polyMeshGenCells(t, points, faces, cells)
-{
-}
+{}
 
 //- Construct from components with the boundary
 polyMeshGen::polyMeshGen
@@ -75,14 +73,12 @@ polyMeshGen::polyMeshGen
         patchStart,
         nFacesInPatch
     )
-{
-}
-        
+{}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Destructor
 polyMeshGen::~polyMeshGen()
-{
-}
+{}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

meshLibrary/utilities/meshes/polyMeshGen/polyMeshGenCells.C

@@ -233,7 +233,18 @@ void polyMeshGenCells::calculateOwnersAndNeighbours() const
 void polyMeshGenCells::calculateAddressingData() const
 {
     if( !ownerPtr_ || !neighbourPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "inline label polyMeshGenCells::calculateAddressingData() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calculateOwnersAndNeighbours();
+    }
 
     addressingDataPtr_ = new polyMeshGenAddressing(*this);
 }
@@ -314,7 +325,18 @@ polyMeshGenCells::~polyMeshGenCells()
 const polyMeshGenAddressing& polyMeshGenCells::addressingData() const
 {
     if( !addressingDataPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "inline label polyMeshGenCells::addressingData() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calculateAddressingData();
+    }
 
     return *addressingDataPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingCellCells.C

@@ -27,7 +27,9 @@ License
 #include "polyMeshGenAddressing.H"
 #include "VRWGraphSMPModifier.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -133,7 +135,18 @@ void polyMeshGenAddressing::calcCellCells() const
 const VRWGraph& polyMeshGenAddressing::cellCells() const
 {
     if( !ccPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const VRWGraph& polyMeshGenAddressing::cellCells() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcCellCells();
+    }
 
     return *ccPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingCellEdges.C

@@ -28,7 +28,9 @@ License
 #include "DynList.H"
 #include "VRWGraphSMPModifier.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -126,7 +128,18 @@ void polyMeshGenAddressing::calcCellEdges() const
 const VRWGraph& polyMeshGenAddressing::cellEdges() const
 {
     if( !cePtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const VRWGraph& polyMeshGenAddressing::cellEdges() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcCellEdges();
+    }
 
     return *cePtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingCellPoints.C

@@ -27,7 +27,9 @@ License
 #include "polyMeshGenAddressing.H"
 #include "VRWGraphSMPModifier.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -121,7 +123,18 @@ void polyMeshGenAddressing::calcCellPoints() const
 const VRWGraph& polyMeshGenAddressing::cellPoints() const
 {
     if( !cpPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const VRWGraph& polyMeshGenAddressing::cellPoints() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcCellPoints();
+    }
 
     return *cpPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingCentresAndAreas.C

@@ -124,7 +124,18 @@ void polyMeshGenAddressing::makeFaceCentresAndAreas
 const vectorField& polyMeshGenAddressing::faceCentres() const
 {
     if( !faceCentresPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const vectorField& polyMeshGenAddressing::faceCentres() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcFaceCentresAndAreas();
+    }
 
     return *faceCentresPtr_;
 }
@@ -132,7 +143,18 @@ const vectorField& polyMeshGenAddressing::faceCentres() const
 const vectorField& polyMeshGenAddressing::faceAreas() const
 {
     if( !faceAreasPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const vectorField& polyMeshGenAddressing::faceAreas() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcFaceCentresAndAreas();
+    }
 
     return *faceAreasPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingCentresAndVols.C

@@ -33,7 +33,9 @@ Description
 
 #include "polyMeshGenAddressing.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -125,7 +127,18 @@ void polyMeshGenAddressing::makeCellCentresAndVols
 const vectorField& polyMeshGenAddressing::cellCentres() const
 {
     if( !cellCentresPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const vectorField& polyMeshGenAddressing::cellCentres() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcCellCentresAndVols();
+    }
 
     return *cellCentresPtr_;
 }
@@ -133,7 +146,18 @@ const vectorField& polyMeshGenAddressing::cellCentres() const
 const scalarField& polyMeshGenAddressing::cellVolumes() const
 {
     if( !cellVolumesPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const scalarField& polyMeshGenAddressing::cellVolumes() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcCellCentresAndVols();
+    }
 
     return *cellVolumesPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingEdgeCells.C

@@ -59,7 +59,18 @@ void polyMeshGenAddressing::calcEdgeCells() const
 const VRWGraph& polyMeshGenAddressing::edgeCells() const
 {
     if( !ecPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const VRWGraph& polyMeshGenAddressing::edgeCells() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcEdgeCells();
+    }
 
     return *ecPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingEdgeFaces.C

@@ -137,7 +137,18 @@ void polyMeshGenAddressing::calcEdgeFaces() const
 const VRWGraph& polyMeshGenAddressing::edgeFaces() const
 {
     if( !efPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const VRWGraph& polyMeshGenAddressing::edgeFaces() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcEdgeFaces();
+    }
 
     return *efPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingEdges.C

@@ -29,7 +29,9 @@ License
 #include "helperFunctions.H"
 #include "demandDrivenData.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -156,7 +158,18 @@ void polyMeshGenAddressing::calcEdges() const
 const edgeList& polyMeshGenAddressing::edges() const
 {
     if( !edgesPtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const edgeList& polyMeshGenAddressing::edges() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcEdges();
+    }
 
     return *edgesPtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingFaceEdges.C

@@ -27,7 +27,9 @@ License
 #include "polyMeshGenAddressing.H"
 #include "VRWGraphSMPModifier.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -111,7 +113,18 @@ void polyMeshGenAddressing::calcFaceEdges() const
 const VRWGraph& polyMeshGenAddressing::faceEdges() const
 {
     if( !fePtr_ )
+    {
+        # ifdef USE_OMP
+        if( omp_in_parallel() )
+            FatalErrorIn
+            (
+                "const VRWGraph& polyMeshGenAddressing::faceEdges() const"
+            ) << "Calculating addressing inside a parallel region."
+                << " This is not thread safe" << exit(FatalError);
+        # endif
+
         calcFaceEdges();
+    }
 
     return *fePtr_;
 }

meshLibrary/utilities/meshes/polyMeshGenAddressing/polyMeshGenAddressingPointCells.C

@@ -45,11 +45,11 @@ void polyMeshGenAddressing::calcPointCells() const
     else
     {
         const VRWGraph& cellPoints = this->cellPoints();
-        
+
         //- create the storage
         pcPtr_ = new VRWGraph();
         VRWGraph& pointCellAddr = *pcPtr_;
-    
+
         VRWGraphSMPModifier(pointCellAddr).reverseAddressing(cellPoints);
         pointCellAddr.setSize(mesh_.points().size());
     }
@@ -60,7 +60,18 @@ void polyMeshGenAddressing::calcPointCells() const
 const VRWGraph& polyMeshGenAddressing::pointCells() const
 {
     if( !pcPtr_ )
+    {