refineMesh: Added customizable direction fields to multiDirRefinement

Contribution provided by Bruno Santos
Resolves feature request http://www.openfoam.org/mantisbt/view.php?id=2004

applications/utilities/mesh/manipulation/refineMesh/refineMeshDict

@@ -22,8 +22,12 @@ set c0;
 //   x,y,z axis. Specify in globalCoeffs section below.
 // - patchLocal : coordinate system different for every cell. Specify in
 //   patchLocalCoeffs section below.
+// - fieldBased : uses the list of field names from the directions list for
+//   selecting the directions to cut. Meant to be used with geometricCut, but
+//   can also be used with useHexTopology.
 coordinateSystem global;
 //coordinateSystem patchLocal;
+//coordinateSystem fieldBased;
 
 // .. and its coefficients. x,y in this case. (normal direction is calculated
 // as tan1^tan2)
@@ -39,7 +43,7 @@ patchLocalCoeffs
     tan1 (1 0 0);
 }
 
-// List of directions to refine
+// List of directions to refine, if global or patchLocal
 directions
 (
     tan1
@@ -47,6 +51,15 @@ directions
     normal
 );
 
+// List of directions to refine, if "fieldBased". Keep in mind that these
+// fields must be of type "vectorField", not "volVectorField".
+//directions
+//(
+//    radialDirectionFieldName
+//    angularDirectionFieldName
+//    heightDirectionFieldName
+//);
+
 // Whether to use hex topology. This will
 // - if patchLocal: all cells on selected patch should be hex
 // - split all hexes in 2x2x2 through the middle of edges.

src/dynamicMesh/meshCut/cellCuts/cellCuts.C

@@ -38,7 +38,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(cellCuts, 0);
+    defineTypeNameAndDebug(cellCuts, 0);
 }
 
 
@@ -120,7 +120,7 @@ void Foam::cellCuts::writeUncutOBJ
     const label cellI
 ) const
 {
-    //- Cell edges
+    // Cell edges
     OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");
 
     Pout<< "Writing cell for time " <<  mesh().time().timeName()
@@ -135,7 +135,7 @@ void Foam::cellCuts::writeUncutOBJ
         labelList(1, cellI)
     );
 
-    //- Loop cutting cell in two
+    // Loop cutting cell in two
     OFstream cutStream(dir / "cellCuts_" + name(cellI) + ".obj");
 
     Pout<< "Writing raw cuts on cell for time " <<  mesh().time().timeName()
@@ -180,7 +180,7 @@ void Foam::cellCuts::writeOBJ
     const labelList& anchors
 ) const
 {
-    //- Cell edges
+    // Cell edges
     OFstream cutsStream(dir / "cell_" + name(cellI) + ".obj");
 
     Pout<< "Writing cell for time " <<  mesh().time().timeName()
@@ -196,7 +196,7 @@ void Foam::cellCuts::writeOBJ
     );
 
 
-    //- Loop cutting cell in two
+    // Loop cutting cell in two
     OFstream loopStream(dir / "cellLoop_" + name(cellI) + ".obj");
 
     Pout<< "Writing loop for time " <<  mesh().time().timeName()
@@ -207,7 +207,7 @@ void Foam::cellCuts::writeOBJ
     writeOBJ(loopStream, loopPoints, vertI);
 
 
-    //- Anchors for cell
+    // Anchors for cell
     OFstream anchorStream(dir / "anchors_" + name(cellI) + ".obj");
 
     Pout<< "Writing anchors for time " <<  mesh().time().timeName()
@@ -2240,7 +2240,7 @@ void Foam::cellCuts::setFromCellCutter
 
                 WarningInFunction
                     << "Found loop on cell " << cellI
-                    << " that resulted in an unexpected bad cut."
+                    << " that resulted in an unexpected bad cut." << nl
                     << "    Suggestions:" << nl
                     << "      - Turn on the debug switch for 'cellCuts' to get"
                     << " geometry files that identify this cell." << nl

src/dynamicMesh/meshCut/directions/directions.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,13 +57,12 @@ const Foam::NamedEnum<Foam::directions::directionType, 3>
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// For debugging
 void Foam::directions::writeOBJ(Ostream& os, const point& pt)
 {
     os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
 }
 
-// For debugging
+
 void Foam::directions::writeOBJ
 (
     Ostream& os,
@@ -81,7 +80,6 @@ void Foam::directions::writeOBJ
 }
 
 
-// Dump to file.
 void Foam::directions::writeOBJ
 (
     const fileName& fName,
@@ -137,7 +135,6 @@ void Foam::directions::check2D
 }
 
 
-// Get direction on all cells
 Foam::vectorField Foam::directions::propagateDirection
 (
     const polyMesh& mesh,
@@ -283,34 +280,35 @@ Foam::directions::directions
     List<vectorField>(wordList(dict.lookup("directions")).size())
 {
     const wordList wantedDirs(dict.lookup("directions"));
+    const word coordSystem(dict.lookup("coordinateSystem"));
 
     bool wantNormal = false;
     bool wantTan1 = false;
     bool wantTan2 = false;
+    label nDirs = 0;
 
-    forAll(wantedDirs, i)
+    if (coordSystem != "fieldBased")
     {
-        directionType wantedDir = directionTypeNames_[wantedDirs[i]];
-
-        if (wantedDir == NORMAL)
+        forAll(wantedDirs, i)
         {
-            wantNormal = true;
-        }
-        else if (wantedDir == TAN1)
-        {
-            wantTan1 = true;
-        }
-        else if (wantedDir == TAN2)
-        {
-            wantTan2 = true;
+            directionType wantedDir = directionTypeNames_[wantedDirs[i]];
+
+            if (wantedDir == NORMAL)
+            {
+                wantNormal = true;
+            }
+            else if (wantedDir == TAN1)
+            {
+                wantTan1 = true;
+            }
+            else if (wantedDir == TAN2)
+            {
+                wantTan2 = true;
+            }
         }
     }
 
 
-    label nDirs = 0;
-
-    const word coordSystem(dict.lookup("coordinateSystem"));
-
     if (coordSystem == "global")
     {
         const dictionary& globalDict = dict.subDict("globalCoeffs");
@@ -424,12 +422,29 @@ Foam::directions::directions
             this->operator[](nDirs++) = tan2Dirs;
         }
     }
+    else if (coordSystem == "fieldBased")
+    {
+        forAll(wantedDirs, i)
+        {
+            operator[](nDirs++) =
+                vectorIOField
+                (
+                    IOobject
+                    (
+                        mesh.instance()/wantedDirs[i],
+                        mesh,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    )
+                );
+        }
+    }
     else
     {
         FatalErrorInFunction
             << "Unknown coordinate system "
             << coordSystem << endl
-            << "Known types are global and patchLocal"
+            << "Known types are global, patchLocal and fieldBased"
             << exit(FatalError);
     }
 }

src/dynamicMesh/meshCut/directions/directions.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,7 +30,8 @@ Description
 
     Used in splitting cells.
     Either all cells have similar refinement direction ('global') or
-    direction is dependent on local cell geometry. Controlled by dictionary.
+    direction is dependent on local cell geometry, or loads selected fields
+    by name ('fieldBased'). Controlled by dictionary.
 
 SourceFiles
     directions.C
@@ -65,7 +66,9 @@ class directions
 :
     public List<vectorField>
 {
+
 public:
+
     // Data types
 
         //- Enumeration listing the possible coordinate directions.
@@ -76,6 +79,7 @@ public:
             NORMAL
         };
 
+
 private:
 
         static const NamedEnum<directionType, 3> directionTypeNames_;
@@ -83,7 +87,6 @@ private:
 
     // Private Member Functions
 
-
         //- For debugging. Write point coordinate.
         static void writeOBJ(Ostream& os, const point& pt);
 
@@ -141,7 +144,6 @@ public:
             const dictionary& dict,
             const twoDPointCorrector* correct2DPtr = NULL
         );
-
 };
 
 

tutorials/mesh/refineMesh/refineFieldDirs/Allclean

+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+rm -rf 0 > /dev/null 2>&1
+
+cleanCase
+
+#------------------------------------------------------------------------------

tutorials/mesh/refineMesh/refineFieldDirs/Allrun

+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+wmake calcRadiusField
+wclean calcRadiusField
+
+runApplication blockMesh
+
+##### Procedure for special refinement over Z
+
+# We need the 0 folder to exist for these steps
+mkdir 0
+
+# Refine over Z, in 6 passes
+for index in 1 2 3 4 5 6; do
+
+  runApplication calcRadiusField
+  mv log.calcRadiusField log.calcRadiusField.tier$index
+
+  runApplication topoSet -dict system/topoSetDict.tier$index
+  mv log.topoSet log.topoSet.tier$index
+
+  ## foamToVTK -cellSet tier$index
+
+  runApplication refineMesh -dict system/refineMeshDict.tier$index -overwrite
+  mv log.refineMesh log.refineMesh.tier$index
+
+  rm -r 0/*
+
+done
+
+# Refine over cylindrical coordinates, in 3 passes
+for index in 1 2 3; do
+
+  runApplication calcRadiusField -calcDirections
+  mv log.calcRadiusField log.calcRadiusField.range$index
+
+  runApplication topoSet -dict system/topoSetDict.range$index
+  mv log.topoSet log.topoSet.range$index
+
+  ## foamToVTK -cellSet tier$index
+
+  runApplication refineMesh -dict system/refineMeshDict.range$index \
+      -overwrite
+  mv log.refineMesh log.refineMesh.range$index
+
+  rm -r 0/*
+
+done
+
+#------------------------------------------------------------------------------

tutorials/mesh/refineMesh/refineFieldDirs/calcRadiusField/Make/files

+calcRadiusField.C
+
+EXE = $(FOAM_USER_APPBIN)/calcRadiusField

tutorials/mesh/refineMesh/refineFieldDirs/calcRadiusField/Make/options

+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume

tutorials/mesh/refineMesh/refineFieldDirs/calcRadiusField/calcRadiusField.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is derived from OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Write the volScalarField "radiusFieldXY" that has the distance to the
+    origin over X,Y.
+
+    And also write the direction fields based on the option "-calcDirections".
+    The resulting fields are:
+      - radialDirection
+      - angularDirection
+      - heightDirection
+
+    Derived from:
+      $FOAM_UTILITIES/postProcessing/miscellaneous/writeCellCentres
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "timeSelector.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "vectorIOField.H"
+#include "volFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    timeSelector::addOptions();
+    #include "addRegionOption.H"
+
+    argList::addBoolOption
+    (
+        "calcDirections",
+        "calculate the direction fields as well"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    instantList timeDirs = timeSelector::select0(runTime, args);
+    const bool calcDirections = args.optionFound("calcDirections");
+
+    #include "createNamedMesh.H"
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+
+        Info<< "Time = " << runTime.timeName() << endl;
+
+        // Check for new mesh
+        mesh.readUpdate();
+
+        Info<< "Writing radius field over X,Y in "
+            <<  runTime.timeName() << endl;
+
+        volScalarField radiusFieldXY
+        (
+            IOobject
+            (
+                "radiusFieldXY",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            sqrt
+            (
+                mesh.C().component(0)*mesh.C().component(0)
+              + mesh.C().component(1)*mesh.C().component(1)
+            )
+        );
+        radiusFieldXY.write();
+
+
+        if(calcDirections)
+        {
+
+            vectorIOField radialDirection
+            (
+                IOobject
+                (
+                    "radialDirection",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh.C()/magSqr(mesh.C())
+            );
+            radialDirection.replace(vector::Z, scalar(0.0));
+            radialDirection /= sqrt(magSqr(radialDirection));
+            radialDirection.write();
+
+
+            const tensor transform2Tangencial
+            (
+                0, -1, 0,
+                1,  0, 0,
+                0,  0, 1
+            );
+            vectorIOField angularDirection
+            (
+                IOobject
+                (
+                    "angularDirection",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                transform2Tangencial & mesh.C()
+            );
+            angularDirection.replace(vector::Z, scalar(0.0));
+            angularDirection /= sqrt(magSqr(angularDirection));
+            angularDirection.write();
+
+            vectorIOField heightDirection
+            (
+                IOobject
+                (
+                    "heightDirection",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                radialDirection ^ angularDirection
+            );
+            heightDirection.write();
+        }
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

tutorials/mesh/refineMesh/refineFieldDirs/system/blockMeshDict

+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant/polyMesh";
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+/*
+
+ A ------------------------------- B
+  /  H                          /
+ C\                            /
+     \                        /
+        \                    /
+           \                /     G
+              \            /
+                 \        /
+                    \    /
+                       \/
+                         D
+*/
+
+Ax  0.438912;    Ay  0.000000;
+Bx  18.28800;    By  0.000000;
+Cx  0.310358;    Cy  -0.310358;
+Dx  12.931569;   Dy  -12.931569;
+
+Hx  0.405502;    Hy  -0.167964;
+Gx  16.895909;   Gy  -6.998515;
+
+Z_DB_low  -0.88706;
+Z_AC_low  -1.07;
+Z_high  4.39208;
+
+vertices
+(
+   ($Cx  $Cy  $Z_AC_low)  //0
+   ($Dx  $Dy  $Z_DB_low)  //1
+   ($Bx  $By  $Z_DB_low)  //2
+   ($Ax  $Ay  $Z_AC_low)  //3
+   ($Cx  $Cy  $Z_high)    //4
+   ($Dx  $Dy  $Z_high)    //5
+   ($Bx  $By  $Z_high)    //6
+   ($Ax  $Ay  $Z_high)    //7
+);
+
+blocks
+(
+   hex (0 1 2 3 4 5 6 7)    (47 10 4) simpleGrading (41.6669 1 1)
+);
+
+edges
+(
+  arc 0 3 ($Hx $Hy $Z_AC_low)
+  arc 4 7 ($Hx $Hy $Z_high)
+
+  arc 1 2 ($Gx $Gy $Z_DB_low)
+  arc 5 6 ($Gx $Gy $Z_high)
+);
+
+patches
+(
+   patch maxX
+   (
+       (1 2 6 5)
+   )
+   patch minX
+   (
+       (0 4 7 3)
+   )
+   patch maxY
+   (
+       (3 7 6 2)
+   )
+   patch minY
+   (
+       (1 5 4 0)
+   )
+   patch maxZ
+   (
+       (4 5 6 7)
+   )
+   patch minZ
+   (
+       (0 3 2 1)
+   )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //