Work in progress commit

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

meshLibrary/utilities/containers/DynList/DynList.H

@@ -101,6 +101,9 @@ public:
         //- Construct given size
         explicit inline DynList(const label);
 
+        //- Construct from given size and defualt value
+        explicit inline DynList(const label, const T&);
+
         //- Construct from UList. nextFree_ set to size().
         explicit inline DynList(const UList<T>&);
 

meshLibrary/utilities/containers/DynList/DynListI.H

@@ -94,7 +94,6 @@ inline Foam::DynList<T, staticSize>::DynList()
 {}
 
 
-//- Construct given size
 template<class T, Foam::label staticSize>
 inline Foam::DynList<T, staticSize>::DynList(const label s)
 :
@@ -104,7 +103,18 @@ inline Foam::DynList<T, staticSize>::DynList(const label s)
     setSize(s);
 }
 
-//- Construct given size
+template<class T, Foam::label staticSize>
+inline Foam::DynList<T, staticSize>::DynList(const label s, const T& val)
+:
+    UList<T>(staticData_, staticSize),
+    nextFree_(0)
+{
+    setSize(s);
+
+    for(label i=0;i<s;++i)
+        this->operator[](i) = val;
+}
+
 template<class T, Foam::label staticSize>
 inline Foam::DynList<T, staticSize>::DynList(const UList<T>& ul)
 :

meshLibrary/utilities/helperFunctions/helperFunctionsTopologyManipulation.H

@@ -119,7 +119,8 @@ namespace help
     inline label sharedVertex(const faceListType& fcs);
 
     //- check if selected elements are in one singly-connected chain
-    inline bool areElementsInChain(const boolList& sel);
+    template<class boolListType>
+    inline bool areElementsInChain(const boolListType& sel);
 
     //- sort the edge chain such that the end of one edge is the beginning
     //- of the next one (valid for singly-connected chains, only)

meshLibrary/utilities/helperFunctions/helperFunctionsTopologyManipulationI.H

@@ -423,9 +423,10 @@ inline label sharedVertex(const faceListType& fcs)
     return -1;
 }
 
-inline bool areElementsInChain(const boolList& sel)
+template<class boolListType>
+inline bool areElementsInChain(const boolListType& sel)
 {
-    boolList selInChain(sel.size(), false);
+    DynList<bool> selInChain(sel.size(), false);
 
     forAll(sel, eI)
     {

meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.H

@@ -89,7 +89,12 @@ class meshSurfaceOptimizer
         inline const VRWGraph& pointTriangles() const;
 
         //- find vertices which are inverted
-        label findInvertedVertices(boolList& smoothVertex) const;
+        //- mark additional layers of vertices near the inverted ones
+        label findInvertedVertices
+        (
+            boolList& smoothVertex,
+            const label nAdditionalLayers = 2
+        ) const;
 
         //- transform into a 2D space in plane
         inline bool transformIntoPlane
@@ -254,7 +259,7 @@ public:
 
     // Member Functions
         //- optimisation before the edges and corners are created
-        bool preOptimizeSurface();
+        bool preOptimizeSurface(const label nAdditionalLayers = 2);
 
         //- optimize boundary nodes after boundary regions are created
         void optimizeSurface(const label nIterations = 5);

meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizeSurface.C

@@ -56,7 +56,8 @@ namespace Foam
 
 label meshSurfaceOptimizer::findInvertedVertices
 (
-    boolList& smoothVertex
+    boolList& smoothVertex,
+    const label nAdditionalLayers
 ) const
 {
     const labelList& bPoints = surfaceEngine_.boundaryPoints();
@@ -93,7 +94,7 @@ label meshSurfaceOptimizer::findInvertedVertices
         return 0;
 
     //- add additional layers around inverted points
-    for(label i=0;i<2;++i)
+    for(label i=0;i<nAdditionalLayers;++i)
     {
         boolList originallySelected = smoothVertex;
         forAll(smoothVertex, bpI)
@@ -152,7 +153,7 @@ label meshSurfaceOptimizer::findInvertedVertices
     return nInvertedTria;
 }
 
-bool meshSurfaceOptimizer::preOptimizeSurface()
+bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
 {
     Info << "Optimizing positions of surface nodes" << endl;
 
@@ -193,7 +194,7 @@ bool meshSurfaceOptimizer::preOptimizeSurface()
 
         do
         {
-            nInvertedTria = findInvertedVertices(smoothVertex);
+            nInvertedTria = findInvertedVertices(smoothVertex, nLayers);
 
             if( nInvertedTria == 0 ) break;
 
@@ -466,6 +467,8 @@ void meshSurfaceOptimizer::optimizeSurface(const label nIterations)
     }
 
     Info << endl;
+
+    preOptimizeSurface(0);
 }
 
 void meshSurfaceOptimizer::optimizeSurface2D(const label nIterations)

meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C

@@ -1304,7 +1304,7 @@ bool edgeExtractor::checkFacePatches()
         //- in the current patch
         # ifdef USE_OMP
         # pragma omp parallel for schedule(dynamic, 40) \
-        if( candidates.size() > 1000 ) reduction(+ : nCorrected)
+        reduction(+ : nCorrected)
         # endif
         forAll(candidates, i)
         {
@@ -1314,7 +1314,7 @@ bool edgeExtractor::checkFacePatches()
             //- do not change patches of faces where all points are mapped
             //- onto the same patch
             bool allInSamePatch(true);
-            forAll(bFaces[bfI], pI)
+            forAll(bf, pI)
                 if( pointPatch_[bp[bf[pI]]] != newBoundaryPatches[bfI] )
                 {
                     allInSamePatch = false;
@@ -1551,10 +1551,8 @@ bool edgeExtractor::findCornerCandidates()
 
 void edgeExtractor::projectDeterminedFeatureVertices()
 {
-    VRWGraph pointPatches;
-    pointPatches.setSizeAndRowSize(pointValence_);
-    forAll(pointPatches, pI)
-        pointPatches.setRowSize(pI, 0);
+    List<DynList<label, 5> > pointPatches;
+    pointPatches.setSize(pointValence_.size());
 
     const meshSurfaceEngine& mse = surfaceEngine();
     const pointFieldPMG& points = mse.mesh().points();
@@ -1568,7 +1566,7 @@ void edgeExtractor::projectDeterminedFeatureVertices()
         const face& bf = bFaces[bfI];
 
         forAll(bf, pI)
-            pointPatches.appendIfNotIn(bp[bf[pI]], facePatch_[bfI]);
+            pointPatches[bp[bf[pI]]].appendIfNotIn(facePatch_[bfI]);
     }
 
     if( Pstream::parRun() )
@@ -1590,6 +1588,8 @@ void edgeExtractor::projectDeterminedFeatureVertices()
         {
             const label bpI = it();
 
+            const DynList<label, 5>& pPatches = pointPatches[bpI];
+
             forAllRow(bpAtProcs, bpI, i)
             {
                 const label neiProc = bpAtProcs(bpI, i);
@@ -1599,8 +1599,8 @@ void edgeExtractor::projectDeterminedFeatureVertices()
 
                 LongList<labelPair>& data = exchangeData[neiProc];
 
-                forAllRow(pointPatches, bpI, ppI)
-                    data.append(labelPair(it.key(), pointPatches(bpI, ppI)));
+                forAll(pPatches, ppI)
+                    data.append(labelPair(it.key(), pPatches[ppI]));
             }
         }
 
@@ -1613,7 +1613,7 @@ void edgeExtractor::projectDeterminedFeatureVertices()
         {
             const labelPair& lp = receivedData[i];
 
-            pointPatches.appendIfNotIn(globalToLocal[lp.first()], lp.second());
+            pointPatches[globalToLocal[lp.first()]].appendIfNotIn(lp.second());
         }
     }
 
@@ -1626,7 +1626,9 @@ void edgeExtractor::projectDeterminedFeatureVertices()
         # endif
         forAll(pointPatches, bpI)
         {
-            if( pointPatches.sizeOfRow(bpI) < 2 )
+            const DynList<label, 5>& pPatches = pointPatches[bpI];
+
+            if( pPatches.size() < 2 )
                 continue;
 
             const point& p = points[bPoints[bpI]];
@@ -1634,7 +1636,7 @@ void edgeExtractor::projectDeterminedFeatureVertices()
             point newP(vector::zero);
             label counter(0);
 
-            forAllRow(pointPatches, bpI, i)
+            forAll(pPatches, i)
             {
                 point pp;
                 scalar dSq;
@@ -1642,7 +1644,7 @@ void edgeExtractor::projectDeterminedFeatureVertices()
                 (
                     pp,
                     dSq,
-                    pointPatches(bpI, i),
+                    pPatches[i],
                     p
                 );
 
@@ -1673,8 +1675,8 @@ bool edgeExtractor::untangleSurface()
 
 void edgeExtractor::extractEdges()
 {
-//    bool changed;
-//    label nIter(0);
+    bool changed;
+    label nIter(0);
 
     distributeBoundaryFaces();
 
@@ -1694,24 +1696,10 @@ void edgeExtractor::extractEdges()
     deleteDemandDrivenData(sPtr);
     # endif
 
-//    do
-//    {
-//        changed = false;
-
-/*        Info << "Checking cells near concave edges" << endl;
-        if( checkConcaveEdgeCells() )
-        {
-            # ifdef DEBUGEdgeExtractor
-            Info << "Changes due to concave edge cells" << endl;
-            fileName sName("changedConcaveCells"+help::scalarToText(nIter)+".stl");
-            sPtr = surfaceWithPatches();
-            sPtr->writeSurface(sName);
-            deleteDemandDrivenData(sPtr);
-            # endif
+    do
+    {
+        changed = false;
 
-            //changed = true;
-        }
-*/
         Info << "Checking patch in the neighbourhood of each face" << endl;
         if( checkFacePatches() )
         {
@@ -1723,23 +1711,23 @@ void edgeExtractor::extractEdges()
             deleteDemandDrivenData(sPtr);
             # endif
 
-            //changed = true;
+            changed = true;
         }
 
         //findEdgeCandidates();
-/*
-        if( findCornerCandidates() )
-        {
-            # ifdef DEBUGEdgeExtractor
-            Info << "Changes due to corner candidates" << endl;
-            fileName sName("findCornerCandidates"+help::scalarToText(nIter)+".stl");
-            sPtr = surfaceWithPatches();
-            sPtr->writeSurface(sName);
-            deleteDemandDrivenData(sPtr);
-            # endif
 
-            changed = true;
-        }
+//        if( findCornerCandidates() )
+//        {
+//            # ifdef DEBUGEdgeExtractor
+//            Info << "Changes due to corner candidates" << endl;
+//            fileName sName("findCornerCandidates"+help::scalarToText(nIter)+".stl");
+//            sPtr = surfaceWithPatches();
+//            sPtr->writeSurface(sName);
+//            deleteDemandDrivenData(sPtr);
+//            # endif
+
+//            changed = true;
+//        }
 
         projectDeterminedFeatureVertices();
 
@@ -1759,9 +1747,8 @@ void edgeExtractor::extractEdges()
         {
             break;
         }
-        */
 
-//    } while( changed && ++nIter < 3 );
+    } while( changed && ++nIter < 3 );
 
     # ifdef DEBUGEdgeExtractor
     const triSurf* sPtr = surfaceWithPatches();

meshLibrary/utilities/surfaceTools/meshSurfaceCheckEdgeTypes/meshSurfaceCheckEdgeTypes.C

@@ -36,7 +36,9 @@ Description
 #include "labelledPoint.H"
 
 #include <map>
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -140,10 +142,7 @@ void meshSurfaceCheckEdgeTypes::classifyEdges()
                 const label f1 = edgeFaces(edgeI, 1);
 
                 if( facePatch[f0] == facePatch[f1] )
-                {
                     edgeTypes_[edgeI] |= PATCHEDGE;
-                    continue;
-                }
 
                 const edge e = edges[edgeI];
 

meshLibrary/utilities/surfaceTools/meshSurfaceCheckInvertedVertices/meshSurfaceCheckInvertedVertices.C

@@ -33,7 +33,10 @@ Description
 #include "helperFunctionsPar.H"
 
 #include <map>
+
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,7 +49,6 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
 {
     const meshSurfaceEngine& mse = *surfaceEnginePtr_;
     const pointFieldPMG& points = mse.points();
-    const labelList& bp = mse.bp();
     const VRWGraph& pointFaces = mse.pointFaces();
     const VRWGraph& pointInFaces = mse.pointInFaces();
     const faceList::subList& bFaces = mse.boundaryFaces();
@@ -64,6 +66,8 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
         if( activePointsPtr_ && !activePointsPtr_->operator[](bpI) )
             continue;
 
+        const vector& pNormal = pNormals[bpI];
+
         forAllRow(pointFaces, bpI, pfI)
         {
             const label pI = pointInFaces(bpI, pfI);
@@ -79,20 +83,20 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
                 fCentres[bfI]
             );
 
-            vector n = triNext.normal();
-            scalar m = mag(n);
-            if( m < VSMALL )
+            vector nNext = triNext.normal();
+            scalar mNext = mag(nNext);
+            if( mNext < VSMALL )
             {
                 # ifdef USE_OMP
                 # pragma omp critical
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
             else
             {
-                n /= m;
+                nNext /= mNext;
             }
 
             if( magSqr(triNext.a() - triNext.b()) < VSMALL )
@@ -102,7 +106,7 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
             if( magSqr(triNext.c() - triNext.a()) < VSMALL )
             {
@@ -111,17 +115,17 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
 
-            if( (n & pNormals[bp[bf[pI]]]) < 0.0 )
+            if( (nNext & pNormal) < 0.0 )
             {
                 # ifdef USE_OMP
                 # pragma omp critical
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
 
             //- check the second triangle (with previous node)
@@ -132,20 +136,20 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
                 points[bf.prevLabel(pI)]
             );
 
-            n = triPrev.normal();
-            m = mag(n);
-            if( m < VSMALL )
+            vector nPrev = triPrev.normal();
+            scalar mPrev = mag(nPrev);
+            if( mPrev < VSMALL )
             {
                 # ifdef USE_OMP
                 # pragma omp critical
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
             else
             {
-                n /= m;
+                nPrev /= mPrev;
             }
 
             if( magSqr(triPrev.a() - triPrev.b()) < VSMALL )
@@ -155,7 +159,7 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
             if( magSqr(triPrev.c() - triPrev.a()) < VSMALL )
             {
@@ -164,27 +168,29 @@ void meshSurfaceCheckInvertedVertices::checkVertices()
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
 
-            if( (n & pNormals[bp[bf[pI]]]) < 0.0 )
+            if( (nPrev & pNormal) < 0.0 )
             {
                 # ifdef USE_OMP
                 # pragma omp critical
                 # endif
                 invertedVertices_.insert(bf[pI]);
 
-                continue;
+                break;
             }
 
             //- check whether the normals of both triangles
             //- point in the same direction
-            if( (triNext.normal() & triPrev.normal()) < 0.0 )
+            if( (nNext & nPrev) < 0.0 )
             {
                 # ifdef USE_OMP
                 # pragma omp critical
                 # endif
                 invertedVertices_.insert(bf[pI]);
+
+                break;
             }
         }
     }

meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractorNonTopo/meshSurfaceEdgeExtractorNonTopo.C

@@ -31,6 +31,7 @@ Description
 
 #include "correctEdgesBetweenPatches.H"
 #include "meshSurfaceOptimizer.H"
+#include "meshSurfaceCheckInvertedVertices.H"
 
 // #define DEBUGSearch
 
@@ -53,9 +54,9 @@ meshSurfaceEdgeExtractorNonTopo::meshSurfaceEdgeExtractorNonTopo
 {
     distributeBoundaryFaces();
 
-    correctEdgesBetweenPatches featureEdges(mesh);
-
     remapBoundaryPoints();
+
+    correctEdgesBetweenPatches featureEdges(mesh);
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

meshLibrary/utilities/surfaceTools/meshSurfaceEngine/meshSurfaceEngineCalculateBoundaryNodesAndFaces.C

@@ -745,7 +745,7 @@ void meshSurfaceEngine::calculatePointNormals() const
     pointNormalsPtr_ = new vectorField(pFaces.size());
 
     # ifdef USE_OMP
-    # pragma omp parallel for if( pFaces.size() > 1000 )
+    # pragma omp parallel for if( pFaces.size() > 1000 ) schedule(dynamic, 50)
     # endif
     forAll(pFaces, pI)
     {