Merge branch 'enhancement-tangledMeshSubsets-t20' into development

meshLibrary/cartesian2DMesh/cartesian2DMeshGenerator/cartesian2DMeshGenerator.C

@@ -214,25 +214,35 @@ void cartesian2DMeshGenerator::renumberMesh()
 
 void cartesian2DMeshGenerator::generateMesh()
 {
-    createCartesianMesh();
+    try
+    {
+        createCartesianMesh();
 
-    surfacePreparation();
+        surfacePreparation();
 
-    mapMeshToSurface();
+        mapMeshToSurface();
 
-    mapEdgesAndCorners();
+        mapEdgesAndCorners();
 
-    optimiseMeshSurface();
+        optimiseMeshSurface();
 
-    generateBoundaryLayers();
+        generateBoundaryLayers();
 
-    optimiseMeshSurface();
+        optimiseMeshSurface();
 
-    refBoundaryLayers();
+        refBoundaryLayers();
 
-    renumberMesh();
+        renumberMesh();
 
-    replaceBoundaries();
+        replaceBoundaries();
+    }
+    catch(...)
+    {
+        WarningIn
+        (
+            "void cartesian2DMeshGenerator::generateMesh()"
+        ) << "Meshing process terminated!" << endl;
+    }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C

@@ -224,12 +224,30 @@ void cartesianMeshGenerator::refBoundaryLayers()
 void cartesianMeshGenerator::optimiseFinalMesh()
 {
     //- untangle the surface if needed
-    meshSurfaceEngine mse(mesh_);
-    meshSurfaceOptimizer(mse, *octreePtr_).optimizeSurface();
+    bool enforceConstraints(false);
+    if( meshDict_.found("enforceGeometryConstraints") )
+    {
+        enforceConstraints =
+            readBool(meshDict_.lookup("enforceGeometryConstraints"));
+    }
+
+    if( true )
+    {
+        meshSurfaceEngine mse(mesh_);
+        meshSurfaceOptimizer surfOpt(mse, *octreePtr_);
+
+        if( enforceConstraints )
+            surfOpt.enforceConstraints();
+
+        surfOpt.optimizeSurface();
+    }
+
     deleteDemandDrivenData(octreePtr_);
 
     //- final optimisation
     meshOptimizer optimizer(mesh_);
+    if( enforceConstraints )
+        optimizer.enforceConstraints();
     optimizer.optimizeMeshFV();
     optimizer.optimizeLowQualityFaces();
     optimizer.untangleMeshFV();
@@ -262,25 +280,35 @@ void cartesianMeshGenerator::renumberMesh()
 
 void cartesianMeshGenerator::generateMesh()
 {
-    createCartesianMesh();
+    try
+    {
+        createCartesianMesh();
 
-    surfacePreparation();
+        surfacePreparation();
 
-    mapMeshToSurface();
+        mapMeshToSurface();
 
-    mapEdgesAndCorners();
+        mapEdgesAndCorners();
 
-    optimiseMeshSurface();
+        optimiseMeshSurface();
 
-    generateBoundaryLayers();
+        generateBoundaryLayers();
 
-    optimiseFinalMesh();
+        optimiseFinalMesh();
 
-    refBoundaryLayers();
+        refBoundaryLayers();
 
-    renumberMesh();
+        renumberMesh();
 
-    replaceBoundaries();
+        replaceBoundaries();
+    }
+    catch(...)
+    {
+        WarningIn
+        (
+            "void cartesianMeshGenerator::generateMesh()"
+        ) << "Meshing process terminated!" << endl;
+    }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

meshLibrary/hexMesh/hexMeshGenerator/hexMeshGenerator.C

@@ -179,8 +179,17 @@ void hexMeshGenerator::refBoundaryLayers()
 
 void hexMeshGenerator::optimiseFinalMesh()
 {
-    //- final optimisation
+    //- check if strict enforcing of constraints is active
+    bool enforceConstraints(false);
+    if( meshDict_.found("enforceGeometryConstraints") )
+    {
+        enforceConstraints =
+            readBool(meshDict_.lookup("enforceGeometryConstraints"));
+    }
+
     meshOptimizer optimizer(mesh_);
+    if( enforceConstraints )
+        optimizer.enforceConstraints();
 
     optimizer.optimizeSurface(*octreePtr_);
 
@@ -216,21 +225,31 @@ void hexMeshGenerator::renumberMesh()
 
 void hexMeshGenerator::generateMesh()
 {
-    generateDualMesh();
+    try
+    {
+        generateDualMesh();
 
-    surfacePreparation();
+        surfacePreparation();
 
-    mapMeshToSurface();
+        mapMeshToSurface();
 
-    optimiseMeshSurface();
+        optimiseMeshSurface();
 
-    generateBoundaryLayers();
+        generateBoundaryLayers();
 
-    optimiseFinalMesh();
+        optimiseFinalMesh();
 
-    renumberMesh();
+        renumberMesh();
 
-    replaceBoundaries();
+        replaceBoundaries();
+    }
+    catch(...)
+    {
+        WarningIn
+        (
+            "void hexMeshGenerator::generateMesh()"
+        ) << "Meshing process terminated!" << endl;
+    }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

meshLibrary/tetMesh/tetMeshGenerator/tetMeshGenerator.C

@@ -169,7 +169,16 @@ void tetMeshGenerator::generateBoudaryLayers()
 void tetMeshGenerator::optimiseFinalMesh()
 {
     //- final optimisation
+    bool enforceConstraints(false);
+    if( meshDict_.found("enforceGeometryConstraints") )
+    {
+        enforceConstraints =
+            readBool(meshDict_.lookup("enforceGeometryConstraints"));
+    }
+
     meshOptimizer optimizer(mesh_);
+    if( enforceConstraints )
+        optimizer.enforceConstraints();
 
     optimizer.optimizeSurface(*octreePtr_);
 
@@ -223,25 +232,35 @@ void tetMeshGenerator::renumberMesh()
 
 void tetMeshGenerator::generateMesh()
 {
-    createTetMesh();
+    try
+    {
+        createTetMesh();
 
-    surfacePreparation();
+        surfacePreparation();
 
-    mapMeshToSurface();
+        mapMeshToSurface();
 
-    mapEdgesAndCorners();
+        mapEdgesAndCorners();
 
-    optimiseMeshSurface();
+        optimiseMeshSurface();
 
-    generateBoudaryLayers();
+        generateBoudaryLayers();
 
-    optimiseFinalMesh();
+        optimiseFinalMesh();
 
-    refBoundaryLayers();
+        refBoundaryLayers();
 
-    renumberMesh();
+        renumberMesh();
 
-    replaceBoundaries();
+        replaceBoundaries();
+    }
+    catch(...)
+    {
+        WarningIn
+        (
+            "void tetMeshGenerator::generateMesh()"
+        ) << "Meshing process terminated!" << endl;
+    }
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

meshLibrary/utilities/checkMeshDict/checkMeshDict.C

@@ -109,6 +109,12 @@ void checkMeshDict::checkBasicSettings() const
             readBool(meshDict_.lookup("checkForGluedMesh"));
         }
     }
+
+    //- check if enforceConstraints is available
+    if( meshDict_.found("enforceGeometryConstraints") )
+    {
+        readBool(meshDict_.lookup("enforceGeometryConstraints"));
+    }
 }
 
 void checkMeshDict::checkPatchCellSize() const

meshLibrary/utilities/smoothers/geometry/meshOptimizer/meshOptimizer.C

@@ -139,7 +139,9 @@ meshOptimizer::meshOptimizer(polyMeshGen& mesh)
 :
     mesh_(mesh),
     vertexLocation_(mesh.points().size(), INSIDE),
-    msePtr_(NULL)
+    msePtr_(NULL),
+    enforceConstraints_(false),
+    badPointsSubsetName_()
 {
     const meshSurfaceEngine& mse = meshSurface();
     const labelList& bPoints = mse.boundaryPoints();
@@ -177,6 +179,15 @@ meshOptimizer::~meshOptimizer()
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+void meshOptimizer::enforceConstraints(const word subsetName)
+{
+    enforceConstraints_ = true;
+
+    badPointsSubsetName_ = subsetName;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 } // End namespace Foam
 
 // ************************************************************************* //

meshLibrary/utilities/smoothers/geometry/meshOptimizer/meshOptimizer.H

@@ -63,77 +63,83 @@ class meshOptimizer
     // Private data
         //- reference to the mesh
         polyMeshGen& mesh_;
-    
+
         //- location of vertex (internal, boundary, edge, corner)
         List<direction> vertexLocation_;
-    
+
         //- mesh surface
         mutable meshSurfaceEngine* msePtr_;
 
+        //- enforce constraints
+        bool enforceConstraints_;
+
+        //- name of the subset contaning tangled points
+        word badPointsSubsetName_;
+
     // Private member functions
         //- return mesh surface
         const meshSurfaceEngine& meshSurface() const;
         void clearSurface();
-    
+
         //- find problematic faces
         label findBadFaces(labelHashSet&, const boolList&) const;
         label findLowQualityFaces(labelHashSet&, const boolList&) const;
-    
+
     // Nested classes
-    
+
         class laplaceSmoother
         {
             // Private data
                 //- reference to the mesh
                 polyMeshGen& mesh_;
-            
+
                 //- location of vertex (internal, boundary, edge, corner)
                 const List<direction>& vertexLocation_;
-            
+
             // Private member functions
                 //- smooth the node using the laplacian smoother
                 //- new position is the average of the neighbouring vertices
                 void laplacian(const labelLongList&, const label);
                 void laplacianSurface(const labelLongList&, const label);
-            
+
                 void laplacianParallel
                 (
                     const labelLongList& procPoints,
                     const bool smoothOnlySurfaceNodes = false
                 );
-                
+
                 //- smooth the node using the laplacian smoother
                 //- new position is the average of the centres of faces attached
                 //- to the vertex
                 void laplacianPC(const labelLongList&, const label);
                 void laplacianPCParallel(const labelLongList& procPoints);
-                
+
                 //- smooth the node using the laplacian smoother
                 //- new position is the average of the centres of faces attached
                 //- to the vertex
                 void laplacianWPC(const labelLongList&, const label);
                 void laplacianWPCParallel(const labelLongList& procPoints);
-                
+
                 //- update geometry after smoothing
                 void updateMeshGeometry(const labelLongList& smoothPoints);
-                
+
                 //- Disallow default bitwise copy construct
                 laplaceSmoother(const laplaceSmoother&);
 
                 //- Disallow default bitwise assignment
                 void operator=(const laplaceSmoother&);
-    
+
         public:
-                
+
             // Constructor
-        
+
                 //- Construct from mesh and vertex locations
                 laplaceSmoother(polyMeshGen&, const List<direction>&);
-        
+
             // Destructor
-        
+
                 ~laplaceSmoother();
-        
+
             // Member Functions
                 //- new position is the average of the neighbouring vertices
                 void optimizeLaplacian(const label nIterations = 1);
@@ -142,7 +148,7 @@ class meshOptimizer
                     const labelHashSet& badFaces,
                     const label nIterations = 1
                 );
-        
+
                 //- new position of surface point is the average of
                 //- the neighbouring surface vertices
                 void optimizeSurfaceLaplacian
@@ -150,7 +156,7 @@ class meshOptimizer
                     const labelHashSet& badFaces,
                     const label nIterations = 1
                 );
-                
+
                 //- new positions are the average of the centres of the cells
                 //- adjacent to the vertex
                 void optimizeLaplacianPC(const label nIterations = 1);
@@ -159,7 +165,7 @@ class meshOptimizer
                     const labelHashSet& badFaces,
                     const label nIterations = 1
                 );
-                
+
                 //- new positions are the average of the centres of the cells
                 //- adjacent to the vertex weighted by cell volumes
                 void optimizeLaplacianWPC(const label nIterations = 1);
@@ -169,13 +175,13 @@ class meshOptimizer
                     const label nIterations = 1
                 );
         };
-        
+
         //- Disallow default bitwise copy construct
         meshOptimizer(const meshOptimizer&);
 
         //- Disallow default bitwise assignment
         void operator=(const meshOptimizer&);
-        
+
     // enumerators
         enum vertexType_
         {
@@ -199,6 +205,12 @@ public:
         ~meshOptimizer();
 
     // Member Functions
+        //- Flag stopping the meshing process if it is not possible
+        //- to untangle the surface without sacrificing geometry constraints
+        //- Points which must be moved away from the required position are
+        //- stored into a point subset
+        void enforceConstraints(const word subsetName="badPoints");
+
         //- smooth surface vertices
         void optimizeSurface(const meshOctree&);
 

meshLibrary/utilities/smoothers/geometry/meshOptimizer/meshOptimizerOptimizeSurface.C

@@ -44,14 +44,17 @@ namespace Foam
 void meshOptimizer::optimizeSurface(const meshOctree& octree)
 {
     Info << "Optimizing positions of surface nodes" << endl;
-    
+
     meshSurfaceEngine& mse = const_cast<meshSurfaceEngine&>(meshSurface());
     meshSurfaceOptimizer surfaceOptimizer(mse, octree);
 
+    if( enforceConstraints_ )
+        surfaceOptimizer.enforceConstraints(badPointsSubsetName_);
+
     surfaceOptimizer.optimizeSurface();
-    
+
     meshSurfaceMapper(mse, octree).mapVerticesOntoSurfacePatches();
-    
+
     clearSurface();
 
     Info << "Finished optimizing positions of surface nodes" << endl;

meshLibrary/utilities/smoothers/geometry/meshOptimizer/optimizeMeshFV.C

@@ -154,7 +154,36 @@ void meshOptimizer::untangleMeshFV()
 
             //- perform optimisation
             if( nBadFaces == 0 )
+            {
                 break;
+            }
+            else if( enforceConstraints_ )
+            {
+                const label subsetId =
+                    mesh_.addPointSubset(badPointsSubsetName_);
+
+                forAllConstIter(labelHashSet, badFaces, it)
+                {
+                    const face& f = faces[it.key()];
+                    forAll(f, pI)
+                        mesh_.addPointToSubset(subsetId, f[pI]);
+                }
+
+                WarningIn
+                (
+                    "void meshOptimizer::untangleMeshFV()"
+                ) << "Writing mesh with " << badPointsSubsetName_
+                  << " subset. These points cannot be untangled"
+                  << " without sacrificing geometry constraints. Exitting.."
+                  << endl;
+
+                returnReduce(1, sumOp<label>());
+                throw std::logic_error
+                (
+                    "void meshOptimizer::untangleMeshFV()"
+                    "Cannot untangle mesh!!"
+                );
+            }
 
             partTetMesh tetMesh(mesh_, badFaces, 0);
             tetMeshOptimisation tmo(tetMesh);

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

@@ -126,7 +126,9 @@ meshSurfaceOptimizer::meshSurfaceOptimizer
     vertexType_(surface.boundaryPoints().size()),
     partitionerPtr_(new meshSurfacePartitioner(surface)),
     deletePartitioner_(true),
-    triMeshPtr_(NULL)
+    triMeshPtr_(NULL),
+    enforceConstraints_(false),
+    badPointsSubsetName_()
 {
     classifySurfaceVertices();
 }
@@ -142,7 +144,9 @@ meshSurfaceOptimizer::meshSurfaceOptimizer
     vertexType_(surfaceEngine_.boundaryPoints().size()),
     partitionerPtr_(&partitioner),
     deletePartitioner_(false),
-    triMeshPtr_(NULL)
+    triMeshPtr_(NULL),
+    enforceConstraints_(false),
+    badPointsSubsetName_()
 {
     classifySurfaceVertices();
 }
@@ -159,6 +163,12 @@ meshSurfaceOptimizer::~meshSurfaceOptimizer()
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+void meshSurfaceOptimizer::enforceConstraints(const word subsetName)
+{
+    enforceConstraints_ = true;
+
+    badPointsSubsetName_ = subsetName;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -79,6 +79,12 @@ class meshSurfaceOptimizer
         //- mesh of surface triangles needed for some smoothers
         mutable partTriMesh* triMeshPtr_;
 
+        //- enforce constraints
+        bool enforceConstraints_;
+
+        //- name of the subset contaning tangled points
+        word badPointsSubsetName_;
+
     // Private member functions
         //- classify surface vertices as PARTITION, EDGE, CORNER
         void classifySurfaceVertices();
@@ -253,6 +259,12 @@ public:
         ~meshSurfaceOptimizer();
 
     // Member Functions
+        //- Flag stopping the meshing process if it is not possible
+        //- to untangle the surface without sacrificing geometry constraints
+        //- Points which must be moved away from the required position are
+        //- stored into a point subset
+        void enforceConstraints(const word subsetName="badPoints");
+
         //- runs a surface smoother on the selected boundary points
         bool untangleSurface
         (

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

@@ -40,6 +40,7 @@ Description
 #include "FIFOStack.H"
 
 #include <map>
+#include <stdexcept>
 
 # ifdef USE_OMP
 #include <omp.h>
@@ -373,7 +374,40 @@ bool meshSurfaceOptimizer::untangleSurface
             nInvertedTria =
                 findInvertedVertices(smoothVertex, nAdditionalLayers);
 
-            if( nInvertedTria == 0 ) break;
+            if( nInvertedTria == 0 )
+            {
+                break;
+            }
+            else if( enforceConstraints_ && !remapVertex )
+            {
+                polyMeshGen& mesh =
+                    const_cast<polyMeshGen&>(surfaceEngine_.mesh());
+
+                const label subsetId =
+                    mesh.addPointSubset(badPointsSubsetName_);
+
+                forAll(smoothVertex, bpI)
+                    if( smoothVertex[bpI] )
+                        mesh.addPointToSubset(subsetId, bPoints[bpI]);
+
+                WarningIn
+                (
+                    "bool meshSurfaceOptimizer::untangleSurface"
+