Improvements of surface smoothing, feature edge capturing + various

fixes


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

meshLibrary/Make/files

@@ -98,6 +98,7 @@ polyMeshGenModifier = utilities/meshes/polyMeshGenModifier
 polyMeshGenAddressing = utilities/meshes/polyMeshGenAddressing
 polyMeshGenChecks = utilities/meshes/polyMeshGenChecks
 partTetMesh = utilities/meshes/partTetMesh
+partTriMesh = utilities/meshes/partTriMesh
 primitiveMesh = utilities/meshes/primitiveMesh
 triSurf = utilities/meshes/triSurf
 cell = utilities/meshes/primitives/cell
@@ -195,6 +196,11 @@ $(partTetMesh)/partTetMeshAddressing.C
 $(partTetMesh)/partTetMeshParallelAddressing.C
 $(partTetMesh)/partTetMeshSimplex.C
 
+$(partTriMesh)/partTriMesh.C
+$(partTriMesh)/partTriMeshAddressing.C
+$(partTriMesh)/partTriMeshParallelAddressing.C
+$(partTriMesh)/partTriMeshSimplex.C
+
 $(triSurf)/triSurf.C
 $(triSurf)/triSurfPoints.C
 $(triSurf)/triSurfFacets.C

meshLibrary/cartesian2DMesh/cartesian2DMeshGenerator/cartesian2DMeshGenerator.C

@@ -187,7 +187,6 @@ void cartesian2DMeshGenerator::refBoundaryLayers()
         meshSurfaceEngine mse(mesh_);
         meshSurfaceOptimizer optimizer(mse, *octreePtr_);
 
-        optimizer.optimizeSurface2D();
         optimizer.untangleSurface2D();
     }
 }

meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C

@@ -135,7 +135,7 @@ void cartesianMeshGenerator::surfacePreparation()
     # ifdef DEBUG
     mesh_.write();
     returnReduce(1, sumOp<label>());
-    //::exit(EXIT_SUCCESS);
+    ::exit(EXIT_SUCCESS);
     # endif
 }
 

meshLibrary/utilities/boundaryLayers/boundaryLayerCells.C

@@ -63,9 +63,12 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
     const labelList& boundaryFacePatches = mse.boundaryFacePatches();
     const labelList& faceOwners = mse.faceOwners();
     const labelList& bp = mse.bp();
-    const VRWGraph& pointPatches = mse.pointPatches();
     const VRWGraph& pointFaces = mse.pointFaces();
 
+    const meshSurfacePartitioner& mPart = surfacePartitioner();
+    const VRWGraph& pointPatches = mPart.pointPatches();
+
+    //- mark patches which will be extruded into layer cells
     boolList treatPatches(mesh_.boundaries().size(), false);
     forAll(patchLabels, patchI)
     {
@@ -251,7 +254,7 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
             continue;
 
         //- generate faces of the bnd layer cell
-        DynList<DynList<label, 4>, 6> cellFaces;
+        FixedList<FixedList<label, 4>, 6> cellFaces;
         createNewCellFromEdge(e, pKeyI, pKeyJ, cellFaces);
 
         //- store boundary faces
@@ -336,8 +339,13 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
     }
 
     //- create cells for corner nodes
-    std::map<label, DynList<label, 3> > nodePatches;
+    typedef std::map<std::pair<label, label>, label> mPairToLabelType;
+    typedef std::map<label, mPairToLabelType> mPointsType;
+    typedef std::map<label, DynList<label, 3> > ppType;
+
+    ppType nodePatches;
     labelHashSet parPoint;
+
     if( Pstream::parRun() )
     {
         const labelList& bPoints = mse.boundaryPoints();
@@ -346,17 +354,15 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
         const Map<label>& globalToLocal = mse.globalToLocalBndPointAddressing();
 
         std::map<label, labelLongList> facesToSend;
-        std::map<label, DynList<face, 8> > parPointFaces;
-        std::map<label, DynList<label, 3> > parPointPatches;
-        for
-        (
-            std::map<
-                label, std::map<std::pair<label, label>, label>
-            >::const_iterator iter=otherVrts_.begin();
-            iter!=otherVrts_.end();
-            ++iter
-        )
+
+        typedef std::map<label, DynList<DynList<label, 8>, 8> > ppfType;
+
+        ppfType parPointFaces;
+        ppType parPointPatches;
+
+        forAllConstIter(mPointsType, otherVrts_, iter)
         {
+            //- skip points on feature edges
             if( iter->second.size() == 2 )
                 continue;
 
@@ -373,7 +379,10 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
                     const label prI = pProcs(bpI, i);
 
                     if( facesToSend.find(prI) == facesToSend.end() )
-                        facesToSend.insert(std::make_pair(prI, labelLongList()));
+                        facesToSend.insert
+                        (
+                            std::make_pair(prI, labelLongList())
+                        );
 
                     if( prI < pMin )
                         pMin = prI;
@@ -381,21 +390,25 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
 
                 if( Pstream::myProcNo() == pMin )
                 {
-                    DynList<face, 8> pFaces;
-                    DynList<label, 3> pPatches;
-                    forAllRow(pointFaces, bpI, fI)
+                    DynList<label, 3>& pPatches = parPointPatches[bpI];
+                    pPatches.setSize(pointFaces.sizeOfRow(bpI));
+
+                    DynList<DynList<label, 8>, 8>& pFaces = parPointFaces[bpI];
+                    pFaces.setSize(pPatches.size());
+
+                    forAllRow(pointFaces, bpI, pfI)
                     {
-                        const label bfI = pointFaces(bpI, fI);
-                        pPatches.append(boundaryFacePatches[bfI]);
+                        const label bfI = pointFaces(bpI, pfI);
+                        const face& bf = bFaces[bfI];
 
-                        face bf;
-                        bf.setSize(bFaces[bfI].size());
+                        pPatches[pfI] = boundaryFacePatches[bfI];
+
+                        DynList<label, 8>& bfCopy = pFaces[pfI];
+                        bfCopy.setSize(bf.size());
                         forAll(bf, pI)
-                            bf[pI] = globalPointLabel[bp[bf[pI]]];
-                        pFaces.append(bf);
+                            bfCopy[pI] = globalPointLabel[bp[bf[pI]]];
                     }
-                    parPointFaces.insert(std::make_pair(bpI, pFaces));
-                    parPointPatches.insert(std::make_pair(bpI, pPatches));
+
                     continue;
                 }
 
@@ -409,14 +422,15 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
                 //- 4. global labels of face points
                 stp.append(globalPointLabel[bpI]);
                 stp.append(pointFaces.sizeOfRow(bpI));
-                forAllRow(pointFaces, bpI, fI)
+                forAllRow(pointFaces, bpI, pfI)
                 {
-                    const label bfI = pointFaces(bpI, fI);
+                    const label bfI = pointFaces(bpI, pfI);
                     const face& bf = bFaces[bfI];
+
                     stp.append(bf.size());
                     stp.append(boundaryFacePatches[bfI]);
                     forAll(bf, pI)
-                        stp.append(globalPointLabel[bf[pI]]);
+                        stp.append(globalPointLabel[bp[bf[pI]]]);
                 }
             }
         }
@@ -432,7 +446,7 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
             const label nFaces = receivedData[counter++];
             for(label fI=0;fI<nFaces;++fI)
             {
-                face f(receivedData[counter++]);
+                DynList<label, 8> f(receivedData[counter++]);
                 parPointPatches[bpI].append(receivedData[counter++]);
                 forAll(f, pI)
                     f[pI] = receivedData[counter++];
@@ -441,34 +455,28 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
         }
 
         //- sort faces sharing corners at the parallel boundaries
-        for
-        (
-            std::map<label, DynList<face, 8> >::iterator iter=parPointFaces.begin();
-            iter!=parPointFaces.end();
-            ++iter
-        )
+        forAllIter(ppfType, parPointFaces, iter)
         {
-            DynList<face, 8>& pFaces = iter->second;
+            DynList<DynList<label, 8>, 8>& pFaces = iter->second;
             DynList<label, 3>& fPatches = parPointPatches[iter->first];
             const label gpI = globalPointLabel[iter->first];
 
             for(label i=0;i<pFaces.size();++i)
             {
-                const face& bf = pFaces[i];
-                const edge e = bf.faceEdge(bf.which(gpI));
+                const DynList<label, 8>& bf = pFaces[i];
+                const label pos = bf.containsAtPosition(gpI);
+                const edge e(bf[pos], bf[bf.fcIndex(pos)]);
 
                 for(label j=i+1;j<pFaces.size();++j)
                 {
-                    const face& obf = pFaces[j];
-                    if(
-                        (obf.which(e.start()) >= 0) &&
-                        (obf.which(e.end()) >= 0)
-                    )
+                    const DynList<label, 8>& obf = pFaces[j];
+                    if( obf.contains(e.start()) && obf.contains(e.end()) )
                     {
-                        face add;
-                        add.transfer(pFaces[i+1]);
-                        pFaces[i+1].transfer(pFaces[j]);
-                        pFaces[j].transfer(add);
+                        DynList<label, 8> add;
+                        add = pFaces[i+1];
+                        pFaces[i+1] = pFaces[j];
+                        pFaces[j] = add;
+
                         const label pAdd = fPatches[i+1];
                         fPatches[i+1] = fPatches[j];
                         fPatches[j] = pAdd;
@@ -477,7 +485,7 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
                 }
             }
 
-            DynList<label, 3> patchIDs(fPatches.size());
+            DynList<label, 3> patchIDs;
             forAll(fPatches, fpI)
                 patchIDs.appendIfNotIn(fPatches[fpI]);
 
@@ -485,15 +493,8 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
         }
     }
 
-    //- sort out vertices at one processor
-    for
-    (
-        std::map<
-            label, std::map<std::pair<label, label>, label>
-        >::const_iterator iter=otherVrts_.begin();
-        iter!=otherVrts_.end();
-        ++iter
-    )
+    //- sort out point which are not at inter-processor boundaries
+    forAllConstIter(mPointsType, otherVrts_, iter)
     {
         if( iter->second.size() == 2 )
             continue;
@@ -504,7 +505,7 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
         const label bpI = bp[iter->first];
 
         //- ensure correct orientation
-        labelList pFaces(pointFaces.sizeOfRow(bpI));
+        DynList<label> pFaces(pointFaces.sizeOfRow(bpI));
         forAll(pFaces, fI)
             pFaces[fI] = pointFaces(bpI, fI);
 
@@ -539,12 +540,7 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
     }
 
     //- create layer cells for corner nodes
-    for
-    (
-        std::map<label, DynList<label, 3> >::iterator iter=nodePatches.begin();
-        iter!=nodePatches.end();
-        ++iter
-    )
+    forAllIter(ppType, nodePatches, iter)
     {
         const DynList<label, 3>& patchIDs = iter->second;
         DynList<label, 3> pKeys;
@@ -562,10 +558,10 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
             continue;
 
         # ifdef DEBUGLayer
-        Info << "Creating corner cell at point " << iter->first << endl;
+        Pout << "Creating corner cell at point " << iter->first << endl;
         # endif
 
-        DynList<DynList<label, 4>, 6> cellFaces;
+        FixedList<FixedList<label, 4>, 6> cellFaces;
         createNewCellFromNode(iter->first, pKeys, cellFaces);
 
         //- store boundary faces
@@ -619,6 +615,7 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
     polyMeshGenModifier meshModifier(mesh_);
 
     meshModifier.addCells(cellsToAdd);
+
     cellsToAdd.clear();
     meshModifier.reorderBoundaryFaces();
     meshModifier.replaceBoundary

meshLibrary/utilities/boundaryLayers/boundaryLayers.C

@@ -60,6 +60,14 @@ const meshSurfaceEngine& boundaryLayers::surfaceEngine() const
     return *msePtr_;
 }
 
+const meshSurfacePartitioner& boundaryLayers::surfacePartitioner() const
+{
+    if( !meshPartitionerPtr_ )
+        meshPartitionerPtr_ = new meshSurfacePartitioner(surfaceEngine());
+
+    return *meshPartitionerPtr_;
+}
+
 void boundaryLayers::findPatchesToBeTreatedTogether()
 {
     if( geometryAnalysed_ )
@@ -70,16 +78,17 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
 
     const meshSurfaceEngine& mse = surfaceEngine();
 
-    const VRWGraph& pPatches = mse.pointPatches();
     const pointFieldPMG& points = mesh_.points();
     const faceList::subList& bFaces = mse.boundaryFaces();
     const edgeList& edges = mse.edges();
     const VRWGraph& eFaces = mse.edgeFaces();
     const labelList& boundaryFacePatches = mse.boundaryFacePatches();
 
+    const meshSurfacePartitioner& mPart = surfacePartitioner();
+    const VRWGraph& pPatches = mPart.pointPatches();
+
     //- patches must be treated together if there exist a corner where
     //- more than three patches meet
-    meshSurfacePartitioner mPart(mse);
     const labelHashSet& corners = mPart.corners();
     forAllConstIter(labelHashSet, corners, it)
     {
@@ -487,6 +496,7 @@ boundaryLayers::boundaryLayers
 :
     mesh_(mesh),
     msePtr_(NULL),
+    meshPartitionerPtr_(NULL),
     patchWiseLayers_(true),
     terminateLayersAtConcaveEdges_(false),
     patchNames_(),

meshLibrary/utilities/boundaryLayers/boundaryLayers.H

@@ -38,6 +38,7 @@ SourceFiles
 
 #include "polyMeshGenModifier.H"
 #include "meshSurfaceEngine.H"
+#include "meshSurfacePartitioner.H"
 #include "DynList.H"
 #include "labelLongList.H"
 #include "Map.H"
@@ -52,6 +53,7 @@ namespace Foam
 
 // Forward declarations
 class meshSurfaceEngine;
+class meshSurfacePartitioner;
 
 /*---------------------------------------------------------------------------*\
                            Class boundaryLayers Declaration
@@ -65,6 +67,9 @@ class boundaryLayers
         //- pointer to mesh surface engine
         mutable meshSurfaceEngine* msePtr_;
 
+        //- poiner to meshSurfacePartitioner
+        mutable meshSurfacePartitioner* meshPartitionerPtr_;
+
         //- shall I create patch-wise layers (true)
         //- O-topology layer (false)
         bool patchWiseLayers_;
@@ -99,9 +104,12 @@ class boundaryLayers
         bool geometryAnalysed_;
 
     // Private member functions
-        //- Return reference to meshSurfaceEngine
+        //- Return const reference to meshSurfaceEngine
         const meshSurfaceEngine& surfaceEngine() const;
 
+        //- return const reference to meshSurfacePartitioner
+        const meshSurfacePartitioner& surfacePartitioner() const;
+
         //- find if any other patches need to be treated together
         //- with the given one
         void findPatchesToBeTreatedTogether();
@@ -195,7 +203,7 @@ class boundaryLayers
             const edge& e,
             const label pKeyI,
             const label pKeyJ,
-            DynList<DynList<label, 4>, 6>& cellFaces
+            FixedList<FixedList<label, 4>, 6>& cellFaces
         ) const;
 
         //- creating hex cells near corners
@@ -203,7 +211,7 @@ class boundaryLayers
         (
             const label pointI,
             const DynList<label, 3>& pKeys,
-            DynList<DynList<label, 4>, 6>& cellFaces
+            FixedList<FixedList<label, 4>, 6>& cellFaces
         ) const;
 
         //- create a bnd layer for a given patch
@@ -213,6 +221,7 @@ class boundaryLayers
         inline void clearOut()
         {
             deleteDemandDrivenData(msePtr_);
+            deleteDemandDrivenData(meshPartitionerPtr_);
         }
 
     // Enumerators

meshLibrary/utilities/boundaryLayers/boundaryLayersCreateVertices.C

@@ -36,7 +36,10 @@ Description
 #include "labelledScalar.H"
 
 #include <map>
+
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 //#define DEBUGLayer
 
@@ -54,7 +57,8 @@ void boundaryLayers::findPatchVertices
 ) const
 {
     const meshSurfaceEngine& mse = surfaceEngine();
-    const VRWGraph& pPatches = mse.pointPatches();
+    const meshSurfacePartitioner& mPart = surfacePartitioner();
+    const VRWGraph& pPatches = mPart.pointPatches();
 
     pVertices.setSize(pPatches.size());
     pVertices = NONE;
@@ -97,7 +101,6 @@ void boundaryLayers::findPatchVertices
             if( pVertices[bpI] && (bpAtProcs.sizeOfRow(bpI) != 0) )
                 pVertices[bpI] |= PARALLELBOUNDARY;
     }
-
 }
 
 point boundaryLayers::createNewVertex
@@ -112,10 +115,12 @@ point boundaryLayers::createNewVertex
     const faceList::subList& bFaces = mse.boundaryFaces();
     const vectorField& pNormals = mse.pointNormals();
     const VRWGraph& pFaces = mse.pointFaces();
-    const VRWGraph& pPatches = mse.pointPatches();
     const labelList& boundaryFacePatches = mse.boundaryFacePatches();
     const VRWGraph& pointPoints = mse.pointPoints();
 
+    const meshSurfacePartitioner& mPart = surfacePartitioner();
+    const VRWGraph& pPatches = mPart.pointPatches();
+
     const pointFieldPMG& points = mesh_.points();
 
     # ifdef DEBUGLayer
@@ -164,22 +169,20 @@ point boundaryLayers::createNewVertex
                 }
             }
 
-            const scalar magV = mag(v);
-            if( magV > VSMALL )
-                v /= magV;
+            const scalar magV = mag(v) + VSMALL;
+            v /= magV;
 
             normal -= (normal & v) * v;
 
-            const scalar magN = mag(normal);
-            if( magN > VSMALL )
-                normal /= magN;
+            const scalar magN = mag(normal) + VSMALL;
+            normal /= magN;
 
             forAllRow(pointPoints, bpI, ppI)
             {
                 if( patchVertex[pointPoints(bpI, ppI)] )
                     continue;
 
-                vector vec = points[bPoints[pointPoints(bpI, ppI)]] - p;
+                const vector vec = points[bPoints[pointPoints(bpI, ppI)]] - p;
                 const scalar prod = 0.5 * mag(vec & normal);
 
                 if( prod < dist )
@@ -226,10 +229,9 @@ point boundaryLayers::createNewVertex
 
             //- normal vector is co-linear with that edge
             normal = p - points[bPoints[otherVertex]];
-            dist = 0.5 * mag(normal);
+            dist = 0.5 * mag(normal) + VSMALL;
 
-            if( mag(dist) > VSMALL )
-                normal /= 2.0 * dist;
+            normal /= 2.0 * dist;
         }
         else
         {
@@ -303,7 +305,14 @@ point boundaryLayers::createNewVertex
     Info << "Distance is " << dist << endl;
     # endif
 
-    return p - dist * normal;
+    dist = Foam::max(dist, VSMALL);
+
+    const point newP = p - dist * normal;
+
+    if( help::isnan(newP) || help::isinf(newP) )
+        return p;
+
+    return newP;
 }
 
 void boundaryLayers::createNewVertices(const boolList& treatPatches)
@@ -321,7 +330,6 @@ void boundaryLayers::createNewVertices(const boolList& treatPatches)
     if( Pstream::parRun() )
     {
         mse.pointNormals();
-        mse.pointPatches();
         mse.pointPoints();
     }
 
@@ -400,22 +408,26 @@ void boundaryLayers::createNewVertices(const labelList& patchLabels)
     patchKey_ = -1;
 
     const meshSurfaceEngine& mse = surfaceEngine();
-    const VRWGraph& pPatches = mse.pointPatches();
     const labelList& bPoints = mse.boundaryPoints();
     const label nBndPts = bPoints.size();
 
+    const meshSurfacePartitioner& mPart = surfacePartitioner();
+    const VRWGraph& pPatches = mPart.pointPatches();
+
     //- the following is needed for parallel runs
     //- it is ugly, but must stay for now :(
     mse.boundaryFaces();
     mse.pointNormals();
     mse.pointFaces();
-    mse.boundaryFacePatches();
     mse.pointPoints();
 
     pointFieldPMG& points = mesh_.points();
     boolList treatPatches(mesh_.boundaries().size());
     List<direction> patchVertex(nBndPts);
 
+    //- make sure than the points are never re-allocated during the process
+    points.reserve(points.size() + 2 * nBndPts);
+
     //- generate new layer vertices for each patch
     forAll(patchLabels, patchI)
     {
@@ -449,8 +461,7 @@ void boundaryLayers::createNewVertices(const labelList& patchLabels)
         //- go throught the vertices and create the new ones
         labelLongList procPoints;
         # ifdef USE_OMP
-        # pragma omp parallel for if( nBndPts > 100 ) \
-        schedule(dynamic, Foam::max(10, nBndPts/(2*omp_get_num_threads())))
+        # pragma omp parallel for if( nBndPts > 1000 ) schedule(dynamic, 50)
         # endif
         for(label bpI=0;bpI<nBndPts;++bpI)
         {
@@ -573,7 +584,14 @@ void boundaryLayers::createNewVertices(const labelList& patchLabels)
             newP += newPatchPenetrationVector[0];
             newP += newPatchPenetrationVector[1];
 
-            points.append(newP);
+            if( !help::isnan(newP) && !help::isinf(newP) )
+            {
+                points.append(newP);
+            }
+            else
+            {
+                points.append(p);
+            }
             newLabelForVertex_[pointI] = nPoints_;
             ++nPoints_;
         }
@@ -591,7 +609,15 @@ void boundaryLayers::createNewVertices(const labelList& patchLabels)
                         p + newPatchPenetrationVector[i] +
                         newPatchPenetrationVector[j];
 
-                    points.append(np);
+                    if( !help::isnan(np) && !help::isinf(np) )
+                    {
+                        points.append(np);
+                    }
+                    else
+                    {
+                        points.append(p);
+                    }
+
                     m.insert
                     (
                         std::make_pair
@@ -605,7 +631,15 @@ void boundaryLayers::createNewVertices(const labelList& patchLabels)
             }
 
             //- create new position for the existing point
-            points.append(newP);
+            if( !help::isnan(newP) && !help::isinf(newP) )
+            {