Integration into meshing workflows, checking validility and various bug

fixes

meshLibrary/cartesian2DMesh/cartesian2DMeshGenerator/cartesian2DMeshGenerator.C

@@ -49,6 +49,8 @@ Description
 #include "checkNonMappableCellConnections.H"
 #include "checkBoundaryFacesSharingTwoEdges.H"
 #include "triSurfaceMetaData.H"
+#include "polyMeshGenGeometryModification.H"
+#include "surfaceMeshGeometryModification.H"
 
 //#define DEBUG
 
@@ -167,6 +169,24 @@ void cartesian2DMeshGenerator::generateBoundaryLayers()
 
     bl.addLayerForAllPatches();
 
+    if( modSurfacePtr_ )
+    {
+        polyMeshGenGeometryModification meshMod(mesh_, meshDict_);
+
+        //- revert the mesh into the original space
+        meshMod.revertGeometryModification();
+
+        //- delete modified surface mesh
+        deleteDemandDrivenData(modSurfacePtr_);
+
+        //- delete the octree
+        deleteDemandDrivenData(octreePtr_);
+
+        //- contruct a new octree from the input surface
+        octreePtr_ = new meshOctree(*surfacePtr_, true);
+        meshOctreeCreator(*octreePtr_).createOctreeWithRefinedBoundary(20);
+    }
+
     # ifdef DEBUG
     mesh_.write();
     //::exit(0);
@@ -252,6 +272,7 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
 :
     db_(time),
     surfacePtr_(NULL),
+    modSurfacePtr_(NULL),
     meshDict_
     (
         IOobject
@@ -314,7 +335,18 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
         surfacePtr_ = surfaceWithPatches;
     }
 
-    octreePtr_ = new meshOctree(*surfacePtr_, true);
+    if( meshDict_.found("anisotropicSources") )
+    {
+        surfaceMeshGeometryModification surfMod(*surfacePtr_, meshDict_);
+
+        modSurfacePtr_ = surfMod.modifyGeometry();
+
+        octreePtr_ = new meshOctree(*modSurfacePtr_, true);
+    }
+    else
+    {
+        octreePtr_ = new meshOctree(*surfacePtr_, true);
+    }
 
     meshOctreeCreator(*octreePtr_, meshDict_).createOctreeBoxes();
 
@@ -326,6 +358,7 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
 cartesian2DMeshGenerator::~cartesian2DMeshGenerator()
 {
     deleteDemandDrivenData(surfacePtr_);
+    deleteDemandDrivenData(modSurfacePtr_);
     deleteDemandDrivenData(octreePtr_);
 }
 

meshLibrary/cartesian2DMesh/cartesian2DMeshGenerator/cartesian2DMeshGenerator.H

@@ -61,6 +61,9 @@ class cartesian2DMeshGenerator
         //- pointer to the surface
         const triSurf* surfacePtr_;
 
+        //- pointer to the modified surface mesh
+        const triSurf* modSurfacePtr_;
+
         //- IOdictionary containing information about cell sizes, etc..
         IOdictionary meshDict_;
 

meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C

@@ -47,6 +47,8 @@ Description
 #include "checkNonMappableCellConnections.H"
 #include "checkBoundaryFacesSharingTwoEdges.H"
 #include "triSurfaceMetaData.H"
+#include "polyMeshGenGeometryModification.H"
+#include "surfaceMeshGeometryModification.H"
 
 //#define DEBUG
 
@@ -72,7 +74,7 @@ void cartesianMeshGenerator::createCartesianMesh()
 
     # ifdef DEBUG
     mesh_.write();
-    //::exit(EXIT_SUCCESS);
+    ::exit(EXIT_SUCCESS);
     # endif
 }
 
@@ -220,6 +222,17 @@ void cartesianMeshGenerator::optimiseFinalMesh()
     optimizer.optimizeLowQualityFaces();
     optimizer.untangleMeshFV();
 
+    if( modSurfacePtr_ )
+    {
+        polyMeshGenGeometryModification meshMod(mesh_, meshDict_);
+
+        //- revert the mesh into the original space
+        meshMod.revertGeometryModification();
+
+        //- delete modified surface mesh
+        deleteDemandDrivenData(modSurfacePtr_);
+    }
+
     # ifdef DEBUG
     mesh_.write();
     //::exit(EXIT_SUCCESS);
@@ -286,6 +299,7 @@ cartesianMeshGenerator::cartesianMeshGenerator(const Time& time)
 :
     db_(time),
     surfacePtr_(NULL),
+    modSurfacePtr_(NULL),
     meshDict_
     (
         IOobject
@@ -335,7 +349,18 @@ cartesianMeshGenerator::cartesianMeshGenerator(const Time& time)
         surfacePtr_ = surfaceWithPatches;
     }
 
-    octreePtr_ = new meshOctree(*surfacePtr_);
+    if( meshDict_.found("anisotropicSources") )
+    {
+        surfaceMeshGeometryModification surfMod(*surfacePtr_, meshDict_);
+
+        modSurfacePtr_ = surfMod.modifyGeometry();
+
+        octreePtr_ = new meshOctree(*modSurfacePtr_);
+    }
+    else
+    {
+        octreePtr_ = new meshOctree(*surfacePtr_);
+    }
 
     meshOctreeCreator(*octreePtr_, meshDict_).createOctreeBoxes();
 
@@ -347,6 +372,7 @@ cartesianMeshGenerator::cartesianMeshGenerator(const Time& time)
 cartesianMeshGenerator::~cartesianMeshGenerator()
 {
     deleteDemandDrivenData(surfacePtr_);
+    deleteDemandDrivenData(modSurfacePtr_);
     deleteDemandDrivenData(octreePtr_);
 }
 

meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.H

@@ -62,6 +62,9 @@ class cartesianMeshGenerator
         //- pointer to the surface
         const triSurf* surfacePtr_;
 
+        //- pointer to the modified surface
+        const triSurf* modSurfacePtr_;
+
         //- IOdictionary containing information about cell sizes, etc..
         IOdictionary meshDict_;
 

meshLibrary/tetMesh/tetMeshGenerator/tetMeshGenerator.C

@@ -45,6 +45,8 @@ Description
 #include "triSurfacePatchManipulator.H"
 #include "refineBoundaryLayers.H"
 #include "triSurfaceMetaData.H"
+#include "polyMeshGenGeometryModification.H"
+#include "surfaceMeshGeometryModification.H"
 
 //#define DEBUG
 
@@ -186,6 +188,17 @@ void tetMeshGenerator::optimiseFinalMesh()
     optimizer.optimizeLowQualityFaces();
     optimizer.optimizeMeshFV();
 
+    if( modSurfacePtr_ )
+    {
+        polyMeshGenGeometryModification meshMod(mesh_, meshDict_);
+
+        //- revert the mesh into the original space
+        meshMod.revertGeometryModification();
+
+        //- delete modified surface mesh
+        deleteDemandDrivenData(modSurfacePtr_);
+    }
+
     # ifdef DEBUG
     mesh_.write();
     //::exit(0);
@@ -268,6 +281,7 @@ tetMeshGenerator::tetMeshGenerator(const Time& time)
 :
     runTime_(time),
     surfacePtr_(NULL),
+    modSurfacePtr_(NULL),
     meshDict_
     (
         IOobject
@@ -315,12 +329,20 @@ tetMeshGenerator::tetMeshGenerator(const Time& time)
         surfacePtr_ = surfaceWithPatches;
     }
 
-    octreePtr_ = new meshOctree(*surfacePtr_);
+    if( meshDict_.found("anisotropicSources") )
+    {
+        surfaceMeshGeometryModification surfMod(*surfacePtr_, meshDict_);
+
+        modSurfacePtr_ = surfMod.modifyGeometry();
+
+        octreePtr_ = new meshOctree(*modSurfacePtr_);
+    }
+    else
+    {
+        octreePtr_ = new meshOctree(*surfacePtr_);
+    }
 
-    meshOctreeCreator* octreeCreatorPtr =
-        new meshOctreeCreator(*octreePtr_, meshDict_);
-    octreeCreatorPtr->createOctreeBoxes();
-    deleteDemandDrivenData(octreeCreatorPtr);
+    meshOctreeCreator(*octreePtr_, meshDict_).createOctreeBoxes();
 
     generateMesh();
 }
@@ -331,6 +353,7 @@ tetMeshGenerator::~tetMeshGenerator()
 {
     deleteDemandDrivenData(surfacePtr_);
     deleteDemandDrivenData(octreePtr_);
+    deleteDemandDrivenData(modSurfacePtr_);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

meshLibrary/tetMesh/tetMeshGenerator/tetMeshGenerator.H

@@ -62,6 +62,9 @@ class tetMeshGenerator
         //- pointer to the surface
         const triSurf* surfacePtr_;
 
+        //- pointer to the modified surface mesh
+        const triSurf* modSurfacePtr_;
+
         //- IOdictionary containing information about cell sizes, etc..
         IOdictionary meshDict_;
 

meshLibrary/utilities/anisotropicMeshing/coordinateModification/boxScaling.C

@@ -26,6 +26,7 @@ License
 #include "boxScaling.H"
 #include "addToRunTimeSelectionTable.H"
 #include "boundBox.H"
+#include "plane.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -115,19 +116,11 @@ vector boxScaling::displacement(const point& p) const
     for(direction i=0;i<vector::nComponents;++i)
     {
         const scalar dispVec = lengthVec_[i] * ((1.0/scaleVec_[i]) - 1.0);
+        const scalar t = ((p[i] - pMin_[i]) / lengthVec_[i]);
 
-        if( p[i] < pMin_[i] )
-        {
-            disp[i] = 0.0;
-        }
-        else if( (p[i] >= pMin_[i]) && (p[i] < pMax_[i]) )
-        {
-            disp[i] = ((p[i] - pMin_[i]) / lengthVec_[i]) * dispVec;
-        }
-        else
-        {
-            disp[i] = dispVec;
-        }
+        const scalar tBnd = Foam::max(0.0, Foam::min(t, 1.0));
+
+        disp[i] = tBnd * dispVec;
     }
 
     return disp;
@@ -139,26 +132,48 @@ vector boxScaling::backwardDisplacement(const point& p) const
 
     for(direction i=0;i<vector::nComponents;++i)
     {
-        const scalar dispVec =
-            lengthVec_[i] * (1.0 - scaleVec_[i]);
-
-        if( p[i] < pMin_[i] )
-        {
-            disp[i] = 0.0;
-        }
-        else if( (p[i] >= pMin_[i]) && (p[i] < pMax_[i]) )
-        {
-            disp[i] = -((p[i] - pMin_[i]) / lengthVec_[i]) * dispVec;
-        }
-        else
-        {
-            disp[i] = -dispVec;
-        }
+        const scalar dispVec = lengthVec_[i] * (scaleVec_[i] - 1.0);
+
+        const scalar t = ((p[i] - pMin_[i]) / lengthVec_[i]);
+
+        const scalar tBnd = Foam::max(0.0, Foam::min(t, 1.0));
+
+        disp[i] = tBnd * dispVec;
     }
 
     return disp;
 }
 
+bool boxScaling::combiningPossible() const
+{
+    return true;
+}
+
+void boxScaling::boundingPlanes(PtrList<plane>&pl) const
+{
+    pl.setSize(6);
+    label counter(0);
+    if( Foam::mag(scaleVec_.x() - 1.0) > VSMALL )
+    {
+        pl.set(counter++, new plane(pMin_, vector(1, 0, 0)));
+        pl.set(counter++, new plane(pMax_, vector(1, 0, 1)));
+    }
+
+    if( Foam::mag(scaleVec_.y() - 1.0) > VSMALL )
+    {
+        pl.set(counter++, new plane(pMin_, vector(0, 1, 0)));
+        pl.set(counter++, new plane(pMax_, vector(0, 1, 0)));
+    }
+
+    if( Foam::mag(scaleVec_.z() - 1.0) > VSMALL )
+    {
+        pl.set(counter++, new plane(pMin_, vector(0, 0, 1)));
+        pl.set(counter++, new plane(pMax_, vector(0, 0, 1)));
+    }
+
+    pl.setSize(counter);
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dictionary boxScaling::dict(bool ignoreType) const

meshLibrary/utilities/anisotropicMeshing/coordinateModification/boxScaling.H

@@ -136,6 +136,13 @@ public:
         //- calculate the displacement vector for box scaling
         virtual vector backwardDisplacement(const point&) const;
 
+        //- can this modification object be combined with other ones
+        virtual bool combiningPossible() const;
+
+        //- return that "bounding planes" of the scaling region for
+        //- the given object
+        virtual void boundingPlanes(PtrList<plane>&) const;
+
     //- Return as dictionary of entries
         dictionary dict(bool ignoreType = false) const;
 

meshLibrary/utilities/anisotropicMeshing/coordinateModification/coneScaling.C

@@ -26,6 +26,7 @@ License
 #include "coneScaling.H"
 #include "addToRunTimeSelectionTable.H"
 #include "boundBox.H"
+#include "plane.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -139,32 +140,21 @@ vector coneScaling::displacement(const point& p) const
     const scalar rCone = (1.0 - tBnd) * r0_ + tBnd * r1_;
 
     //- calculate displacement in the axial direction
-    if( t > 1.0 )
-    {
-        //- translate in the axial direction
-        disp += axialVec_ * ((1.0/axialScaling_) - 1.0);
-    }
-    else if( t > 0.0 )
-    {
-        //- scale the distance in the axial direction
-        disp += t * axialVec_ * ((1.0/axialScaling_) - 1.0);
-    }
+    disp += tBnd * axialVec_ * ((1.0/axialScaling_) - 1.0);
+
 
     //- calculate displacement in the radial direction
+
     if( r > VSMALL )
     {
+        const scalar tRadial = r / rCone;
+        const scalar tRadialBnd = Foam::min(1.0, tRadial);
+
         const scalar rScale = (1.0/radialScaling_) - 1.0;
         const vector dispRadial = rCone * (rVec / r) * rScale;
-        if( r > rCone )
-        {
-            //- translate in the radial direction
-            disp += dispRadial;
-        }
-        else
-        {
-            //- scale the distance in the radial direction
-            disp += (r / rCone) * dispRadial;
-        }
+
+        //- scale the distance in the radial direction
+        disp += tRadialBnd * dispRadial;
     }
 
     return disp;
@@ -186,38 +176,50 @@ vector coneScaling::backwardDisplacement(const point& p) const
     const scalar rCone = (1.0 - tBnd) * r0_ + tBnd * r1_;
 
     //- calculate displacement in the axial direction
-    const scalar aScale = 1.0 - axialScaling_;
-    if( t > 1.0 )
-    {
-        //- translate in the axial direction
-        disp -= axialVec_ * aScale;
-    }
-    else if( t > 0.0 )
-    {
-        //- scale the distance in the axial direction
-        disp -= t * axialVec_ * aScale;
-    }
+    const scalar aScale = axialScaling_ - 1.0;
+    disp += tBnd * axialVec_ * aScale;
 
     //- calculate displacement in the radial direction
     if( r > VSMALL )
     {
-        const scalar rScale = 1.0 - radialScaling_;
+        const scalar tRadial = r / rCone;
+        const scalar tRadialBnd = Foam::min(1.0, tRadial);
+
+        const scalar rScale = radialScaling_ - 1.0;
         const vector dispRadial = rCone * (rVec / r) * rScale;
-        if( r > rCone )
-        {
-            //- translate in the radial direction
-            disp -= dispRadial;
-        }
-        else
-        {
-            //- scale the distance in the radial direction
-            disp -= (r / rCone) * dispRadial;
-        }
+
+        //- scale the distance in the radial direction
+        disp += tRadialBnd * dispRadial;
     }
 
     return disp;
 }
 
+bool coneScaling::combiningPossible() const
+{
+    if( Foam::mag(radialScaling_ - 1.0) > VSMALL )
+        return false;
+
+    return true;
+}
+
+void coneScaling::boundingPlanes(PtrList<plane>& pl) const
+{
+    vector n = axialVec_ / (mag(axialVec_) + VSMALL);
+
+    if( Foam::mag(axialScaling_ - 1.0) > VSMALL )
+    {
+        pl.setSize(2);
+
+        pl.set(0, new plane(p0_, n));
+        pl.set(1, new plane(p1_, n));
+    }
+    else
+    {
+        pl.clear();
+    }
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dictionary coneScaling::dict(bool ignoreType) const
@@ -350,6 +352,11 @@ void coneScaling::operator=(const dictionary& d)
     }
     else
     {
+        WarningIn
+        (
+            "void coneScaling::operator=(const dictionary& d)"
+        ) << "Entry radialScaling is not specified!" << endl;
+
         radialScaling_ = 1.0;
     }
 
@@ -360,6 +367,11 @@ void coneScaling::operator=(const dictionary& d)
     }
     else
     {
+        WarningIn
+        (
+            "void coneScaling::operator=(const dictionary& d)"
+        ) << "Entry axialScaling is not specified!" << endl;
+
         axialScaling_ = 1.0;
     }
 

meshLibrary/utilities/anisotropicMeshing/coordinateModification/coneScaling.H

@@ -140,6 +140,13 @@ public:
         //- calculate the displacement vector for box scaling
         virtual vector backwardDisplacement(const point&) const;
 
+        //- can this modification object be combined with other ones
+        virtual bool combiningPossible() const;
+
+        //- return that "bounding planes" of the scaling region for
+        //- the given object
+        virtual void boundingPlanes(PtrList<plane>&) const;
+
     //- Return as dictionary of entries
         dictionary dict(bool ignoreType = false) const;
 

meshLibrary/utilities/anisotropicMeshing/coordinateModification/coordinateModification.H

@@ -38,6 +38,7 @@ SourceFiles
 #include "word.H"
 #include "point.H"
 #include "dictionary.H"
+#include "PtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,6 +47,7 @@ namespace Foam
 
 // Forward declarations
 class boundBox;
+class plane;
 
 /*---------------------------------------------------------------------------*\
                    Class coordinateModification Declaration
@@ -115,6 +117,13 @@ public:
         //- calculate the displacement vector for the backward modification
         virtual vector backwardDisplacement(const point&) const = 0;
 
+        //- can this modification object be combined with other ones
+        virtual bool combiningPossible() const = 0;
+
+        //- return that "bounding planes" of the scaling region for
+        //- the given object
+        virtual void boundingPlanes(PtrList<plane>&) const = 0;
+
 
       // Access
 

meshLibrary/utilities/anisotropicMeshing/coordinateModification/coordinateModifier.C

@@ -24,10 +24,123 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "coordinateModifier.H"
+#include "plane.H"
 
 namespace Foam
 {
 
+void coordinateModifier::checkForValidInverse() const
+{
+    if( modifiers_.size() >= 1 )
+    {
+        //- the if the modifiers allow combinations
+        forAll(modifiers_, modI)
+            if( !modifiers_[modI].combiningPossible() )
+            {
+                FatalErrorIn
+                (
+                    "void coordinateModifier::checkForValidInverse() const"
+                ) << modifiers_[modI].name() << " cannot be combined with"
+                  << " other anisotropic sources. The operation"
+                  << " cannot be reverted!" << exit(FatalError);
+            }
+
+        //- check if the modifications overlap
+        forAll(modifiers_, modI)
+        {
+            PtrList<plane> bndPlanes;
+            modifiers_[modI].boundingPlanes(bndPlanes);
+
+            # ifdef DEBUGCoordinateModifier
+            Info << "Checking planes for object " << modifiers_[modI].name()
+                 << " which are " << bndPlanes << endl;
+            # endif
+
+            for(label modJ=modI+1;modJ<modifiers_.size();++modJ)
+            {
+                PtrList<plane> otherBndPlanes;
+                modifiers_[modJ].boundingPlanes(otherBndPlanes);
+
+                # ifdef DEBUGCoordinateModifier
+                Info << "Bnd planes planes for " << modifiers_[modJ].name()
+                     << " are " << otherBndPlanes << endl;
+                # endif
+
+                for(label i=0;i<bndPlanes.size();i+=2)
+                {
+                    const plane& pl = bndPlanes[i];
+
+                    for(label j=0;j<otherBndPlanes.size();j+=2)
+                    {
+                        const plane& opl = otherBndPlanes[j];
+
+                        const scalar dn = mag(pl.normal() & opl.normal());
+                        if( dn > SMALL )
+                        {
+                            if( dn < (1.0 - SMALL) )
+                            {
+                                FatalErrorIn
+                                (
+                                    "void coordinateModifier::"
+                                    "checkForValidInverse() const"
+                                ) << "Bounding planes of the objects "
+                                  << modifiers_[modI].name()
+                                  << " and " << modifiers_[modJ].name()