Improvements of 2D mesh generator

meshLibrary/Make/files

@@ -81,6 +81,7 @@ checkNonMappableCellConnections = $(topology)/checkNonMappableCellConnections
 
 triSurfaceTools = utilities/triSurfaceTools
 triSurfaceCleanupDuplicates = $(triSurfaceTools)/triSurfaceCleanupDuplicates
+triSurfaceCleanupDuplicateTriangles = $(triSurfaceTools)/triSurfaceCleanupDuplicateTriangles
 triSurfaceCopyParts = $(triSurfaceTools)/triSurfaceCopyParts
 triSurfaceCurvatureEstimator = $(triSurfaceTools)/triSurfaceCurvatureEstimator
 triSurfacePartitioner = $(triSurfaceTools)/triSurfacePartitioner
@@ -439,6 +440,9 @@ $(surfaceIntersectionsOctree)/surfaceIntersectionsOctreeIntersections.C
 $(triSurfaceCleanupDuplicates)/triSurfaceCleanupDuplicates.C
 $(triSurfaceCleanupDuplicates)/triSurfaceCleanupDuplicatesFunctions.C
 
+$(triSurfaceCleanupDuplicateTriangles)/triSurfaceCleanupDuplicateTriangles.C
+$(triSurfaceCleanupDuplicateTriangles)/triSurfaceCleanupDuplicateTrianglesFunctions.C
+
 $(triSurfaceCopyParts)/triSurfaceCopyParts.C
 
 $(triSurfacePartitioner)/triSurfacePartitioner.C

meshLibrary/cartesian2DMesh/cartesian2DMeshGenerator/cartesian2DMeshGenerator.C

@@ -29,6 +29,7 @@ Description
 #include "polyMeshGen2DEngine.H"
 #include "triSurf.H"
 #include "triSurfacePatchManipulator.H"
+#include "triSurfaceCleanupDuplicateTriangles.H"
 #include "demandDrivenData.H"
 #include "Time.H"
 #include "meshOctreeCreator.H"
@@ -277,6 +278,9 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
 
     if( surfacePtr_->featureEdges().size() != 0 )
     {
+        //- get rid of duplicate triangles as they cause strange problems
+        triSurfaceCleanupDuplicateTriangles(const_cast<triSurf&>(*surfacePtr_));
+
         //- create surface patches based on the feature edges
         //- and update the meshDict based on the given data
         triSurfacePatchManipulator manipulator(*surfacePtr_);
@@ -287,6 +291,8 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
         //- delete the old surface and assign the new one
         deleteDemandDrivenData(surfacePtr_);
         surfacePtr_ = surfaceWithPatches;
+
+        surfacePtr_->writeSurface("surfWithPatches.fms");
     }
 
     octreePtr_ = new meshOctree(*surfacePtr_, true);

meshLibrary/utilities/boundaryLayers/boundaryLayerCells.C

@@ -243,6 +243,15 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
         const label pKeyI = patchKey_[patchI];
         const label pKeyJ = patchKey_[patchJ];
 
+        if( pKeyI < 0 || pKeyJ < 0 )
+        {
+            continue;
+            FatalErrorIn
+            (
+                "void boundaryLayers::createLayerCells(const labelList&)"
+            ) << "Patch key is negative at concave edge" << abort(FatalError);
+        }
+
         if( pKeyI == pKeyJ )
             continue;
 
@@ -333,6 +342,10 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
             faceAtPatches.append(labelPair(patchI, patchJ));
         }
 
+        # ifdef DEBUGLayer
+        Info << "Adding new cell at edge " << cellFaces << endl;
+        # endif
+
         //- append cell to the queue
         cellsToAdd.appendGraph(cellFaces);
     }
@@ -628,10 +641,6 @@ void boundaryLayers::createLayerCells(const labelList& patchLabels)
     //- delete meshSurfaceEngine
     this->clearOut();
 
-    # ifdef DEBUGLayer
-    mesh_.addressingData().checkMesh(true);
-    # endif
-
     Info << "Finished creating layer cells" << endl;
 }
 

meshLibrary/utilities/boundaryLayers/boundaryLayers.C

@@ -93,7 +93,7 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
     {
         const label bpI = it.key();
 
-        if( mPart.numberOfFeatureEdgesAtPoint(bpI) != 3 )
+        if( mPart.numberOfFeatureEdgesAtPoint(bpI) > 3 )
         {
             labelHashSet commonPatches;
             DynList<label> allPatches;
@@ -144,7 +144,11 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
 
         //- check if the any of the face vertices is tangled
         const edge& e = edges[eI];
-        if( invertedVertices.found(e[0]) || invertedVertices.found(e[1]) )
+        if
+        (
+            !is2DMesh_ &&
+            (invertedVertices.found(e[0]) || invertedVertices.found(e[1]))
+        )
             continue;
 
         const label patch0 = boundaryFacePatches[eFaces(eI, 0)];
@@ -164,7 +168,11 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
 
             const face& f1 = bFaces[eFaces(eI, 0)];
             const face& f2 = bFaces[eFaces(eI, 1)];
-            if( !help::isSharedEdgeConvex(points, f1, f2) )
+            if
+            (
+                !help::isSharedEdgeConvex(points, f1, f2) ||
+                (help::angleBetweenFaces(points, f1, f2) > 1.5 * M_PI)
+            )
             {
                 ++edgeClassification[pp].second();
             }
@@ -210,7 +218,11 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
                 continue;
 
             const edge& e = edges[beI];
-            if( invertedVertices.found(e[0]) || invertedVertices.found(e[1]) )
+            if
+            (
+                !is2DMesh_ &&
+                (invertedVertices.found(e[0]) || invertedVertices.found(e[1]))
+            )
                 continue;
 
             //- do not send data if the face on other processor
@@ -287,6 +299,7 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
                 }
             }
 
+            const face& bf = bFaces[eFaces(beI, 0)];
             const label patch0 = boundaryFacePatches[eFaces(beI, 0)];
             const label patch1 = otherProcPatches[beI];
 
@@ -303,7 +316,10 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
 
             if(
                 (otherFaceProc[beI] > Pstream::myProcNo()) &&
-                !help::isSharedEdgeConvex(points, bFaces[eFaces(beI, 0)], f)
+                (
+                    !help::isSharedEdgeConvex(points, bf, f) ||
+                    (help::angleBetweenFaces(points, bf, f) > 1.5 * M_PI)
+                )
             )
             {
                 ++edgeClassification[pp].second();
@@ -443,6 +459,34 @@ void boundaryLayers::findPatchesToBeTreatedTogether()
 
     Info << "Patch names " << patchNames_ << endl;
     Info << "Treat patches with patch " << treatPatchesWithPatch_ << endl;
+
+    label layerI(0), subsetId;
+    boolList usedPatch(treatPatchesWithPatch_.size(), false);
+    const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
+
+    forAll(treatPatchesWithPatch_, patchI)
+    {
+        if( usedPatch[patchI] || (boundaries[patchI].patchSize() == 0) )
+            continue;
+
+        Info << "Adding layer subset " << layerI << " for patch " << patchI << endl;
+        usedPatch[patchI] = true;
+        subsetId = mesh_.addFaceSubset("layer_"+help::scalarToText(layerI));
+        ++layerI;
+
+        forAll(treatPatchesWithPatch_[patchI], i)
+        {
+            const label cPatch = treatPatchesWithPatch_[patchI][i];
+            usedPatch[cPatch] = true;
+
+            label start = boundaries[cPatch].patchStart();
+            const label size = boundaries[cPatch].patchSize();
+            for(label i=0;i<size;++i)
+                mesh_.addFaceToSubset(subsetId, start++);
+        }
+    }
+
+    mesh_.write();
     # endif
 
     geometryAnalysed_ = true;
@@ -498,6 +542,7 @@ boundaryLayers::boundaryLayers
     meshPartitionerPtr_(NULL),
     patchWiseLayers_(true),
     terminateLayersAtConcaveEdges_(false),
+    is2DMesh_(false),
     patchNames_(),
     treatedPatch_(),
     treatPatchesWithPatch_(),
@@ -597,6 +642,8 @@ void boundaryLayers::activate2DMode()
             if( treatedPatch_[patches[i]] )
                 patches.removeElement(i);
     }
+
+    is2DMesh_ = true;
 }
 
 void boundaryLayers::addLayerForAllPatches()

meshLibrary/utilities/boundaryLayers/boundaryLayers.H

@@ -76,6 +76,9 @@ class boundaryLayers
         //- shall the layers be terminated at concave edges (true)
         bool terminateLayersAtConcaveEdges_;
 
+        //- is it a 2D mesh
+        bool is2DMesh_;
+
         //- patch names
         wordList patchNames_;
 

meshLibrary/utilities/boundaryLayers/boundaryLayersI.H

@@ -40,12 +40,14 @@ label boundaryLayers::findNewNodeLabel
     const label pKey
 ) const
 {
-    if( otherVrts_.find(pointI) != otherVrts_.end() )
+    const std::map
+    <
+        label, std::map<std::pair<label, label>, label>
+    >::const_iterator it = otherVrts_.find(pointI);
+
+    if( it != otherVrts_.end() )
     {
-        const std::map
-        <
-            label, std::map<std::pair<label, label>, label>
-        >::const_iterator it = otherVrts_.find(pointI);
+
         const std::map<std::pair<label, label>, label>& m = it->second;
 
         std::map<std::pair<label, label>, label>::const_iterator mit;
@@ -67,7 +69,7 @@ label boundaryLayers::findNewNodeLabel
                     continue;
                 if( mit->first.first == mit->first.second )
                     continue;
-                
+
                 return mit->second;
             }
         }
@@ -88,22 +90,22 @@ inline void boundaryLayers::createNewCellFromEdge
         otherVrts_.find(e.start())->second;
     const std::map<std::pair<label, label>, label>& me =
         otherVrts_.find(e.end())->second;
-    
+
     # ifdef DEBUGLayer
-    Info << "Creating cell for edge " << edgeI << " with nodes " << e << endl;
+    Info << "Creating cell for edge with nodes " << e << endl;
     Info << "pKeyI " << pKeyI << endl;
     Info << "pKeyJ " << pKeyJ << endl;
     std::map<std::pair<label, label>, label>::const_iterator iter;
     for(iter=ms.begin();iter!=ms.end();++iter)
-        Info << "1. Pair (" << iter->first.first << ", " 
+        Info << "1. Pair (" << iter->first.first << ", "
             << iter->first.second << ") has value " << iter->second << endl;
     for(iter=me.begin();iter!=me.end();++iter)
-        Info << "2. Pair (" << iter->first.first << ", " 
+        Info << "2. Pair (" << iter->first.first << ", "
             << iter->first.second << ") has value " << iter->second << endl;
     # endif
-    
+
     label p0s(-1), p1s(-1), ns(-1), p0e(-1), p1e(-1), ne(-1);
-    
+
     if( ms.size() == 2 )
     {
         p0s = ms.find(std::pair<label, label>(pKeyI, pKeyI))->second;
@@ -172,7 +174,7 @@ inline void boundaryLayers::createNewCellFromEdge
             }
         }
     }
-    
+
     //- F0
     cellFaces[0][0] = e.end();
     cellFaces[0][1] = e.start();
@@ -208,6 +210,22 @@ inline void boundaryLayers::createNewCellFromEdge
     cellFaces[5][1] = p0s;
     cellFaces[5][2] = ns;
     cellFaces[5][3] = p1s;
+
+    # ifdef DEBUGLayer
+    forAll(cellFaces, fI)
+    {
+        forAll(cellFaces[fI], pI)
+        {
+            if
+            (
+                cellFaces[fI][pI] < 0 ||
+                cellFaces[fI][pI] >= mesh_.points().size()
+            )
+                FatalError << "Invalid point indices found!"
+                           << abort(FatalError);
+        }
+    }
+    # endif
 }
 
 inline void boundaryLayers::createNewCellFromNode
@@ -219,7 +237,7 @@ inline void boundaryLayers::createNewCellFromNode
 {
     const std::map<std::pair<label, label>, label>& m =
         otherVrts_.find(pointI)->second;
-    
+
     //- create labels before creating cells
     const label n = newLabelForVertex_[pointI];
     const label p00 =
@@ -238,7 +256,7 @@ inline void boundaryLayers::createNewCellFromNode
         pr.second = pKeys[0];
     }
     const label p01 = m.find(pr)->second;
-    
+
     pr.first = pKeys[0];
     pr.second = pKeys[2];
     if( m.find(pr) == m.end() )
@@ -247,7 +265,7 @@ inline void boundaryLayers::createNewCellFromNode
         pr.second = pKeys[0];
     }
     const label p02 = m.find(pr)->second;
-    
+
     pr.first = pKeys[1];
     pr.second = pKeys[2];
     if( m.find(pr) == m.end() )
@@ -256,9 +274,9 @@ inline void boundaryLayers::createNewCellFromNode
         pr.second = pKeys[1];
     }
     const label p12 = m.find(pr)->second;
-    
+
     //- create the cell and append it
-    
+
     //- F0
     cellFaces[0][0] = pointI;
     cellFaces[0][1] = p02;

meshLibrary/utilities/helperFunctions/helperFunctionsGeometryQueries.H

@@ -66,7 +66,15 @@ namespace help
     bool isinf(const ListType&);
 
     //- check if the faces share a convex edge
-    bool isSharedEdgeConvex
+    inline bool isSharedEdgeConvex
+    (
+        const pointField& points,
+        const face& f1,
+        const face& f2
+    );
+
+    //- angle between the two faces in radians
+    inline scalar angleBetweenFaces
     (
         const pointField& points,
         const face& f1,

meshLibrary/utilities/helperFunctions/helperFunctionsGeometryQueriesI.H

@@ -77,8 +77,8 @@ inline bool isSharedEdgeConvex
     const face& f2
 )
 {
-    face triOwn(3);
-    face triNei(3);
+    DynList<label, 3> triOwn(3);
+    DynList<label, 3> triNei(3);
 
     forAll(f1, pI)
     {
@@ -98,7 +98,7 @@ inline bool isSharedEdgeConvex
 
         forAll(triOwn, pJ)
         {
-            if( triNei.which(triOwn[pJ]) == -1 )
+            if( triNei.contains(triOwn[pJ]) )
             {
                 tetrahedron<point, point> tet
                 (
@@ -120,6 +120,98 @@ inline bool isSharedEdgeConvex
     return true;
 }
 
+inline scalar angleBetweenFaces
+(
+    const pointField& points,
+    const face& f1,
+    const face& f2
+)
+{
+    DynList<label, 3> triOwn(3);
+    DynList<label, 3> triNei(3);
+
+    scalar angle(0.0);
+    label counter(0);
+
+    forAll(f1, pI)
+    {
+        label pos = f2.which(f1[pI]);
+        if( pos == -1 )
+            continue;
+
+        triNei[0] = f2[pos];
+        triNei[1] = f2.nextLabel(pos);
+        triNei[2] = f2.prevLabel(pos);
+
+        triOwn[0] = f1[pI];
+        triOwn[1] = f1.nextLabel(pI);
+        triOwn[2] = f1.prevLabel(pI);
+
+        scalar vol(0.0);
+
+        forAll(triOwn, pJ)
+        {
+            if( triNei.contains(triOwn[pJ]) )
+            {
+                tetrahedron<point, point> tet
+                (
+                    points[triNei[0]],
+                    points[triNei[1]],
+                    points[triNei[2]],
+                    points[triOwn[pJ]]
+                );
+
+                vol = tet.mag();
+                break;
+            }
+        }
+
+        vector nOwn
+        (
+            (points[triOwn[1]] - points[triOwn[0]]) ^
+            (points[triOwn[2]] - points[triOwn[0]])
+        );
+        nOwn /= (mag(nOwn) + VSMALL);
+
+        vector nNei
+        (
+            (points[triNei[1]] - points[triNei[0]]) ^
+            (points[triNei[2]] - points[triNei[0]])
+        );
+        nNei /= (mag(nNei) + VSMALL);
+
+        const scalar dot = Foam::max(-1.0, Foam::min(nOwn & nNei, 1.0));
+
+        if( vol > -VSMALL )
+        {
+            //- the angle is in the interval [Pi, 2Pi>
+            const scalar ang = Foam::acos(dot);
+
+            angle += ang + M_PI;
+            ++counter;
+        }
+        else
+        {
+            //- the angle is in the interval [0, Pi>
+            const scalar ang = Foam::acos(-dot);
+
+            angle += ang;
+            ++counter;
+        }
+    }
+
+    if( counter == 0 )
+    {
+        FatalErrorIn
+        (
+            "scalar angleBetweenFaces"
+            "(const pointField&, const face&, const face&)"
+        ) << "Faces do no share an edge" << abort(FatalError);
+    }
+
+    return angle / counter;
+}
+
 inline faceList mergePatchFaces
 (
     const List< DynList<label> >& pfcs,

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

@@ -700,7 +700,8 @@ void meshSurfaceOptimizer::optimizeSurface2D(const label nIterations)
     {
         Info << "." << flush;
 
-        smoothLaplacianFC(movedPoints, procBndPoints, false);
+        //smoothLaplacianFC(movedPoints, procBndPoints, false);
+        smoothSurfaceOptimizer(movedPoints);
 
         //- move the points which are not at minimum z coordinate
         mesh2DEngine.correctPoints();

meshLibrary/utilities/triSurfaceTools/triSurfaceCleanupDuplicateTriangles/triSurfaceCleanupDuplicateTriangles.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | cfMesh: A library for mesh generation
+   \\    /   O peration     |
+    \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
+     \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of cfMesh.
+
+    cfMesh is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    cfMesh is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "triSurfaceCleanupDuplicateTriangles.H"
+#include "triSurf.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+triSurfaceCleanupDuplicateTriangles::triSurfaceCleanupDuplicateTriangles
+(
+    triSurf& surf
+)
+:
+    surf_(surf),
+    newTriangleLabel_()
+{
+    checkDuplicateTriangles();
+}
+
+triSurfaceCleanupDuplicateTriangles::~triSurfaceCleanupDuplicateTriangles()
+{}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

meshLibrary/utilities/triSurfaceTools/triSurfaceCleanupDuplicateTriangles/triSurfaceCleanupDuplicateTriangles.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | cfMesh: A library for mesh generation
+   \\    /   O peration     |
+    \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
+     \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of cfMesh.
+
+    cfMesh is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    cfMesh is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    triSurfaceCleanupDuplicates
+
+Description
+    Removes duplicate triangles from the surface.
+
+SourceFiles
+    triSurfaceCleanupDuplicateTriangles.C
+    triSurfaceCleanupDuplicateTrianglesFunctions.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triSurfaceCleanupDuplicateTriangles_H
+#define triSurfaceCleanupDuplicateTriangles_H
+
+#include "VRWGraph.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class triSurf;
+
+/*---------------------------------------------------------------------------*\