diff --git a/applications/utilities/preProcessing/PDR/Allwclean b/applications/utilities/preProcessing/PDR/Allwclean
index 9f1e739332..4478b88ce6 100755
--- a/applications/utilities/preProcessing/PDR/Allwclean
+++ b/applications/utilities/preProcessing/PDR/Allwclean
@@ -1,7 +1,11 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
+#------------------------------------------------------------------------------
wclean libso pdrFields
+
+wclean PDRfitMesh
+
wclean PDRsetFields
#------------------------------------------------------------------------------
diff --git a/applications/utilities/preProcessing/PDR/Allwmake b/applications/utilities/preProcessing/PDR/Allwmake
index 2131eed2f8..10618728d8 100755
--- a/applications/utilities/preProcessing/PDR/Allwmake
+++ b/applications/utilities/preProcessing/PDR/Allwmake
@@ -5,6 +5,9 @@ cd "${0%/*}" || exit # Run from this directory
#------------------------------------------------------------------------------
wmake $targetType pdrFields
+
+wmake $targetType PDRfitMesh
+
wmake $targetType PDRsetFields
#------------------------------------------------------------------------------
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/Make/files b/applications/utilities/preProcessing/PDR/PDRfitMesh/Make/files
new file mode 100644
index 0000000000..9a61526f53
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/Make/files
@@ -0,0 +1,7 @@
+PDRfitMeshParams.C
+PDRfitMeshScan.C
+PDRfitMeshScans.C
+
+PDRfitMesh.C
+
+EXE = $(FOAM_APPBIN)/PDRfitMesh
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/Make/options b/applications/utilities/preProcessing/PDR/PDRfitMesh/Make/options
new file mode 100644
index 0000000000..64d6ba5625
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/Make/options
@@ -0,0 +1,15 @@
+EXE_INC = \
+ -I../pdrFields/lnInclude \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/surfMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/mesh/blockMesh/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lfileFormats \
+ -lsurfMesh \
+ -lmeshTools \
+ -lblockMesh \
+ -lpdrFields
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMesh.C b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMesh.C
new file mode 100644
index 0000000000..b5a47a5ca9
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMesh.C
@@ -0,0 +1,339 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of 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 .
+
+Applications
+ PDRfitMesh
+
+Description
+ Scans extents of obstacles to establish reasonable estimates
+ for generating a PDRblockMeshDict.
+
+SourceFiles
+ PDRfitMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "Fstream.H"
+#include "IOdictionary.H"
+
+#include "PDRfitMeshScans.H"
+#include "PDRsetFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char* argv[])
+{
+ argList::addNote
+ (
+ "Processes a set of geometrical obstructions to determine"
+ " reasonable estimates for generating a PDRblockMeshDict"
+ );
+ argList::noParallel();
+ argList::noFunctionObjects();
+
+ argList::addOption("dict", "file", "Alternative PDRsetFieldsDict");
+ argList::addOption("params", "file", "Alternative PDRfitMeshDict");
+ argList::addBoolOption
+ (
+ "overwrite",
+ "Overwrite existing system/PDRblockMeshDict"
+ );
+
+ argList::addBoolOption("verbose", "Increase verbosity");
+
+ argList::addBoolOption
+ (
+ "dry-run",
+ "Equivalent to -print-dict"
+ );
+
+ argList::addBoolOption
+ (
+ "print-dict",
+ "Print PDRblockMeshDict equivalent and exit"
+ );
+
+ argList::addBoolOption
+ (
+ "write-vtk",
+ "Write obstacles as VTK files"
+ );
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+
+ const bool dryrun = args.found("dry-run");
+ const bool printDict = args.found("print-dict");
+
+ const word dictName("PDRsetFieldsDict");
+ #include "setSystemRunTimeDictionaryIO.H"
+
+ Info<< "Reading " << dictIO.name() << nl << endl;
+
+ IOdictionary setFieldsDict(dictIO);
+
+ const fileName& casepath = runTime.globalPath();
+
+ // Program parameters (globals)
+ pars.read(setFieldsDict);
+
+ // Params for fitMesh
+ // - like setSystemRunTimeDictionaryIO
+
+ IOobject paramIO = IOobject::selectIO
+ (
+ IOobject
+ (
+ "PDRfitMeshDict",
+ runTime.system(),
+ runTime,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ ),
+ args.getOrDefault("params", "")
+ );
+
+ if (paramIO.typeHeaderOk(true))
+ {
+ Info<< "Using PDRfitMesh parameters from "
+ << runTime.relativePath(paramIO.objectPath()) << nl
+ << endl;
+ }
+ else
+ {
+ paramIO.readOpt() = IOobject::NO_READ;
+ Info<< "No PDRfitMeshDict found, using defaults" << nl
+ << endl;
+ }
+
+
+ IOdictionary fitParamsDict(paramIO, dictionary());
+
+ PDRfitMeshParams fitParams(fitParamsDict);
+
+
+ // Essential parameters
+
+ scalar cellWidth = 0;
+
+ if
+ (
+ !setFieldsDict.readIfPresent("cellWidth", cellWidth)
+ && !fitParamsDict.readIfPresent("cellWidth", cellWidth)
+ )
+ {
+ FatalErrorInFunction
+ << "No cellWidth specified in any dictionary" << nl
+ << exit(FatalError);
+ }
+
+
+ // Storage for obstacles and cylinder-like obstacles
+ DynamicList obstacles, cylinders;
+
+ // Read in obstacles
+ if (pars.legacyObsSpec)
+ {
+ PDRobstacle::legacyReadFiles
+ (
+ pars.obsfile_dir, pars.obsfile_names,
+ boundBox::invertedBox, // ie, no bounds
+ obstacles,
+ cylinders
+ );
+ }
+ else
+ {
+ PDRobstacle::readFiles
+ (
+ pars.obsfile_dir, pars.obsfile_names,
+ boundBox::invertedBox, // ie, no bounds
+ obstacles,
+ cylinders
+ );
+ }
+
+
+ //
+ // Output names
+ //
+
+ IOobject outputIO
+ (
+ "PDRblockMeshDict",
+ runTime.system(),
+ runTime,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false // unregistered
+ );
+
+
+ enum class writeType { DRY_RUN, OKAY, OVERWRITE, NO_CLOBBER };
+
+ writeType writable = writeType::OKAY;
+
+ if (dryrun || printDict)
+ {
+ writable = writeType::DRY_RUN;
+ }
+ else if (isFile(outputIO.objectPath()))
+ {
+ if (args.found("overwrite"))
+ {
+ writable = writeType::OVERWRITE;
+ }
+ else
+ {
+ writable = writeType::NO_CLOBBER;
+ InfoErr
+ << nl
+ << "File exists: "
+ << runTime.relativePath(outputIO.objectPath()) << nl
+ << "Move out of the way or specify -overwrite" << nl << nl
+ << "Exiting" << nl << nl;
+
+ return 1;
+ }
+ }
+ else
+ {
+ writable = writeType::OKAY;
+ }
+
+
+ if (args.found("write-vtk"))
+ {
+ PDRobstacle::generateVtk(casepath/"VTK", obstacles, cylinders);
+ }
+
+
+ //
+ // The fitting routines
+ //
+
+ // Collapse into a single list of obstacles
+ obstacles.append(std::move(cylinders));
+
+ PDRfitMeshScan::verbose(args.found("verbose"));
+
+ Vector griding =
+ PDRfitMeshScans().calcGriding(obstacles, fitParams, cellWidth);
+
+
+ //
+ // Output
+ //
+
+ if (writable == writeType::DRY_RUN)
+ {
+ InfoErr << nl;
+ if (!printDict)
+ {
+ InfoErr
+ << "dry-run: ";
+ }
+ InfoErr
+ << "Displaying equivalent PDRblockMeshDict" << nl
+ << nl;
+ }
+ else if (writable == writeType::OKAY)
+ {
+ InfoErr
+ << nl
+ << "Write "
+ << runTime.relativePath(outputIO.objectPath()) << nl;
+ }
+ else if (writable == writeType::OVERWRITE)
+ {
+ InfoErr
+ << nl
+ << "Overwrite existing "
+ << runTime.relativePath(outputIO.objectPath()) << nl;
+ }
+ // NO_CLOBBER already handled
+
+
+ {
+ autoPtr outputFilePtr;
+
+ if
+ (
+ writable == writeType::OKAY
+ || writable == writeType::OVERWRITE
+ )
+ {
+ outputFilePtr.reset(new OFstream(outputIO.objectPath()));
+ }
+
+ Ostream& os = bool(outputFilePtr) ? *outputFilePtr : Info();
+
+ outputIO.writeHeader(os, IOdictionary::typeName_());
+
+ os.writeEntry
+ (
+ "expansion",
+ PDRfitMeshScan::expansionName()
+ );
+ os << nl;
+
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ griding[cmpt].writeDict(os, cmpt);
+ }
+
+ const dictionary* outerDict;
+
+ if
+ (
+ (outerDict = setFieldsDict.findDict("outer")) != nullptr
+ || (outerDict = fitParamsDict.findDict("outer")) != nullptr
+ )
+ {
+ outerDict->writeEntry(os);
+ }
+ else
+ {
+ // Define our own "outer"
+ os.beginBlock("outer");
+ os.writeEntry("type", "none");
+ os.endBlock();
+ }
+
+ IOobject::writeEndDivider(os);
+ }
+
+ Info<< nl << "\nEnd\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshDict b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshDict
new file mode 100644
index 0000000000..682832db2e
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshDict
@@ -0,0 +1,126 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v2012 |
+| \\ / A nd | Website: www.openfoam.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object PDRfitMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Tuning parameters (defaults) for PDRfitMesh
+// - these normally do not need much changing
+
+
+// Note the following are normally read from PDRsetFieldsDict
+//
+// - cellWidth
+// - cellWidthFactor
+// - outer (contains optional zmin)
+//
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ----------------------------------------------------------------------
+// Finding planes
+// ----------------------------------------------------------------------
+
+// User defined bounds
+bounds
+{
+ // xmin -1e16;
+ // xmax 1e16;
+
+ // ymin -1e16;
+ // ymax 1e16;
+
+ // zmin -1e16;
+ // zmax 1e16;
+}
+
+
+outer
+{
+ //- Ground datum (hard zmin)
+ // zmin -1e16;
+}
+
+
+// ----------------------------------------------------------------------
+// Finding planes
+// ----------------------------------------------------------------------
+
+
+// Ignore faces smaller than this multiple of cell face area
+minFaceArea 5.0;
+
+// Only fit a plane if the face area at the coordinate is at least
+// this times the cell face area
+minAreaRatio 20.0;
+
+// Size of search zones for face areas will be this * the cell width.
+// So faces closer than this zone size might be grouped together
+areaWidthFactor 0.7;
+
+
+// Very long zones produce bad outer boundary shape from
+// makePDRmeshBlocks, so we subdivide a zone if its length is greater
+// than this * the height of the area of cuboid vells
+maxZoneToHeight 2.0;
+
+
+
+// ----------------------------------------------------------------------
+// Choosing cell width
+//
+// and read from PDRsetFieldsDict
+// ----------------------------------------------------------------------
+
+// Minimum number of cells in any direction
+nCellsMin 10;
+
+// Width of obstacles must be less than this * cell width to added
+// into subgrid length
+widthFactor 1.0;
+
+// Optimal average number of obstacles per cell
+obsPerCell 2.0;
+
+//- Maximum cellWidth when auto-sizing
+maxCellWidth 1e16;
+
+// Do not use a cellWidth more than this times the initial estimate
+maxWidthEstimate 5.0;
+
+// Assume converged if optimised width changes by less than this
+maxWidthRatio 1.2;
+
+// Maximum iterations in optimising cellWidth
+maxIterations 5;
+
+
+// ----------------------------------------------------------------------
+// Defining outer region
+//
+// ----------------------------------------------------------------------
+
+// Fraction of rectangular cell layers each side of the central region.
+nEdgeLayers 5.0;
+
+
+// FUTURE?
+//
+// // Ratio of the outer extension distance to the average radius of the
+// // congested region
+// ratioOuterToCongRad 20.0;
+//
+// // Cell size ratio - should be same as that hard-wired in makePDRmeshBlocks
+// outerRatio 1.2;
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.C b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.C
new file mode 100644
index 0000000000..0680c3ef50
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.C
@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRfitMeshParams.H"
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::PDRfitMeshParams::readBounds(const dictionary& dict)
+{
+ scalar value;
+
+ if (dict.readIfPresent("xmin", value))
+ {
+ minBounds.x() = optionalData(value);
+ }
+ if (dict.readIfPresent("ymin", value))
+ {
+ minBounds.y() = optionalData(value);
+ }
+ if (dict.readIfPresent("zmin", value))
+ {
+ minBounds.z() = optionalData(value);
+ }
+
+ if (dict.readIfPresent("xmax", value))
+ {
+ maxBounds.x() = optionalData(value);
+ }
+ if (dict.readIfPresent("ymax", value))
+ {
+ maxBounds.y() = optionalData(value);
+ }
+ if (dict.readIfPresent("zmax", value))
+ {
+ maxBounds.z() = optionalData(value);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::PDRfitMeshParams::PDRfitMeshParams(const dictionary& dict)
+{
+ read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::PDRfitMeshParams::read(const dictionary& dict)
+{
+ const dictionary* dictptr;
+
+ // Geometric limits
+
+ if ((dictptr = dict.findDict("bounds")) != nullptr)
+ {
+ readBounds(*dictptr);
+ }
+
+ if ((dictptr = dict.findDict("outer")) != nullptr)
+ {
+ const dictionary& d = *dictptr;
+
+ d.readIfPresent("zmin", ground);
+ }
+
+
+ // Finding planes
+
+ dict.readIfPresent("minFaceArea", minFaceArea);
+ dict.readIfPresent("minAreaRatio", minAreaRatio);
+ dict.readIfPresent("areaWidthFactor", areaWidthFactor);
+ dict.readIfPresent("maxZoneToHeight", maxZoneToHeight);
+
+
+ // Choosing cell width
+
+ dict.readIfPresent("nCellsMin", nCellsMin);
+ dict.readIfPresent("widthFactor", widthFactor);
+ dict.readIfPresent("obsPerCell", obsPerCell);
+ dict.readIfPresent("maxCellWidth", maxCellWidth);
+ dict.readIfPresent("maxWidthEstimate", maxWidthEstimate);
+ dict.readIfPresent("maxWidthRatio", maxWidthRatio);
+ dict.readIfPresent("maxIterations", maxIterations);
+
+
+ // Outer region
+
+ dict.readIfPresent("nEdgeLayers", nEdgeLayers);
+
+ dict.readIfPresent("outerRatio", outerRatio);
+
+ // dict.readIfPresent("outerRadius", outerRadius);
+};
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.H b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.H
new file mode 100644
index 0000000000..f1a73b69e4
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.H
@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of 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 .
+
+Class
+ Foam::PDRfitMeshParams
+
+Description
+ Parameters for PDRfitMesh
+
+SourceFiles
+ PDRfitMeshParams.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRfitMeshParams_H
+#define PDRfitMeshParams_H
+
+#include "labelList.H"
+#include "scalarList.H"
+#include "dictionary.H"
+#include "optionalData.H"
+#include "MinMax.H"
+#include "Vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class PDRfitMeshParams Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRfitMeshParams
+{
+ // Private Member Functions
+
+ //- Read optional bounds
+ void readBounds(const dictionary& dict);
+
+
+public:
+
+ // Data Members
+
+ //- Optional user bounds (xmin, xmax, ymin, ymax, zmin, zmax)
+ Vector> minBounds, maxBounds;
+
+ //- Ground datum (hard zmin)
+ scalar ground = -1e16;
+
+
+ // Finding planes
+
+ //- Ignore faces smaller than this multiple of cell face area
+ scalar minFaceArea = 5.0;
+
+ //- Only fit a plane if the face area at the coordinate is at
+ //- least this times the cell face area
+ scalar minAreaRatio = 20.0;
+
+ //- Size of search zones for face areas will be this * the
+ //- cell width.
+ // Faces closer than this zone size may be grouped together
+ scalar areaWidthFactor = 0.7;
+
+ // Very long zones produce bad outer boundary shape from
+ // makePDRmeshBlocks, so we subdivide a zone if its length is
+ // greater than this * the height of the area of cuboid vells
+ scalar maxZoneToHeight = 2.0;
+
+
+ // Choosing cell width
+
+ // Minimum number of cells in any direction
+ label nCellsMin = 10;
+
+ // Width of obstacles must be less than this * cell width to
+ // added into subgrid length
+ scalar widthFactor = 1.0;
+
+ // Optimal average number of obstacles per cell
+ scalar obsPerCell = 2.0;
+
+ //- Maximum cellWidth when auto-sizing
+ scalar maxCellWidth = 1e16;
+
+ //- Do not use a cellWidth more than this times the initial estimate
+ scalar maxWidthEstimate = 5.0;
+
+ //- Converged if optimised width changes by less than this amount
+ scalar maxWidthRatio = 1.2;
+
+ //- Maximum iterations in optimising cellWidth
+ label maxIterations = 5;
+
+
+ // Outer region
+
+ //- Fraction of rectangular cell layers on each side of the
+ //- central region
+ scalar nEdgeLayers = 5;
+
+ //- Cell size ratio - should be same as used for PDRblockMesh
+ scalar outerRatio = 1.2;
+
+//FUTURE? // Ratio of the outer extension distance to the average
+//FUTURE? // radius of the congested region
+//FUTURE? scalar outerRadius = 20.0;
+
+
+ // Constructors
+
+ //- Default construct
+ PDRfitMeshParams() = default;
+
+ //- Construct and read dictionary
+ explicit PDRfitMeshParams(const dictionary& dict);
+
+
+ // Member Functions
+
+ //- Read program parameters from dictionary
+ void read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.C b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.C
new file mode 100644
index 0000000000..22d3ccbef7
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.C
@@ -0,0 +1,859 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRfitMeshScan.H"
+#include "PDRfitMeshParams.H"
+#include "PDRblock.H"
+#include "PDRobstacle.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+bool Foam::PDRfitMeshScan::verbose_ = false;
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::word Foam::PDRfitMeshScan::expansionName()
+{
+ return PDRblock::expansionNames_[PDRblock::EXPAND_RELATIVE];
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::PDRfitMeshScan::addCoord
+(
+ const scalar coord,
+ const scalar weight
+)
+{
+ if (!limits_.contains(coord))
+ {
+ return;
+ }
+
+ const scalar coordOffset = (coord - limits_.min());
+
+ if (coordOffset < 0)
+ {
+ FatalErrorInFunction
+ << "Negative coordinate! " << coordOffset << nl
+ << exit(FatalError);
+ }
+
+ // We divide the range into a series of steps, and decide
+ // which step this coord is in
+
+
+ // Note: would be simpler to record the totals only in the
+ // half steps in this routine then get the full-step values
+ // at the end by adding half steps
+
+
+ // Add the area to the total already found in this step
+ // and update the weighted average position of these obstacle faces
+
+ const label stepi1 = 2*floor(coordOffset / stepWidth_);
+
+
+ // Also do half-step offset in case several faces
+ // are around the step boundary.
+ //
+ // Stored in odd-numbered elements
+
+ const label stepi2 =
+ Foam::max
+ (
+ 0,
+ 1 + 2*floor((coordOffset - 0.5*stepWidth_) / stepWidth_)
+ );
+
+
+ // Last, keep track of totals in half-step ranges
+
+ const label stepih = floor(coordOffset / stepWidth_);
+
+
+ const label maxSize = Foam::max(stepi1, stepi2);
+
+ if (weightedPos_.size() <= maxSize)
+ {
+ weightedPos_.resize(maxSize+1, Zero);
+ totalArea_.resize(maxSize+1, Zero);
+
+ weightedPos2_.resize(maxSize+1, Zero);
+ totalArea2_.resize(maxSize+1, Zero);
+ }
+
+ weightedPos_[stepi1] += (coord * weight);
+ totalArea_[stepi1] += weight;
+
+ weightedPos_[stepi2] += (coord * weight);
+ totalArea_[stepi2] += weight;
+
+ weightedPos2_[stepih] += (coord * weight);
+ totalArea2_[stepih] += weight;
+}
+
+
+// void Foam::PDRfitMeshScan::clear()
+// {
+// weightedPos_.clear(); totalArea_.clear();
+// weightedPos2_.clear(); totalArea2_.clear();
+// }
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::PDRfitMeshScan::reset()
+{
+ limits_.min() = GREAT;
+ limits_.max() = -GREAT;
+
+ hard_min_ = -GREAT;
+
+ nsteps_ = 0;
+ stepWidth_ = 0;
+ minFaceArea_ = 0;
+
+ weightedPos_.clear();
+ totalArea_.clear();
+
+ weightedPos2_.clear();
+ totalArea2_.clear();
+}
+
+
+void Foam::PDRfitMeshScan::resize
+(
+ const scalar cellWidth,
+ const PDRfitMeshParams& pars
+)
+{
+ if (limits_.valid())
+ {
+ nsteps_ =
+ (
+ limits_.mag() / (mag(cellWidth) * pars.areaWidthFactor)
+ + 0.5
+ );
+
+ nsteps_ = max(nsteps_, pars.nCellsMin);
+
+ stepWidth_ = limits_.mag() / nsteps_;
+
+ minFaceArea_ = pars.minFaceArea;
+ }
+ else
+ {
+ nsteps_ = 0;
+ stepWidth_ = 0;
+ minFaceArea_ = 0;
+
+ WarningInFunction
+ << "No valid limits defined" << endl;
+ }
+
+ weightedPos_.clear();
+ totalArea_.clear();
+
+ weightedPos2_.clear();
+ totalArea2_.clear();
+
+ // resize
+ weightedPos_.resize(2*nsteps_+1, Zero);
+ totalArea_.resize(2*nsteps_+1, Zero);
+
+ weightedPos2_.resize(2*nsteps_+1, Zero);
+ totalArea2_.resize(2*nsteps_+1, Zero);
+}
+
+
+void Foam::PDRfitMeshScan::adjustLimits
+(
+ const scalar point0,
+ const scalar point1
+)
+{
+ limits_.add(point0);
+ limits_.add(point1);
+}
+
+
+void Foam::PDRfitMeshScan::addArea
+(
+ const scalar area,
+ const scalar point0,
+ const scalar point1
+)
+{
+ if (!nsteps_)
+ {
+ FatalErrorInFunction
+ << "No step-size defined" << nl
+ << exit(FatalError);
+ }
+
+ if (area < minFaceArea_ * sqr(stepWidth_))
+ {
+ return;
+ }
+
+ addCoord(point0, area);
+ addCoord(point1, area);
+}
+
+
+void Foam::PDRfitMeshScan::print(Ostream& os) const
+{
+ if (nsteps_)
+ {
+ os << "steps:" << nsteps_;
+ }
+
+ if (limits_.valid())
+ {
+ os << " limits:" << limits_;
+ }
+
+ scalar totalArea = 0;
+ for (const scalar areaValue : totalArea_)
+ {
+ totalArea += areaValue;
+ }
+
+ os << " area:" << totalArea;
+ os << nl;
+}
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Evaluates (r**n - 1) / (r - 1) including in the limit r=1
+inline static scalar rn1r1(const scalar r, const label n)
+{
+ if (mag(r - 1.0) < SMALL)
+ {
+ return n;
+ }
+
+ return (pow(r, n) - 1.0)/(r - 1.0);
+}
+
+
+inline static scalar end_val
+(
+ const scalar width,
+ const scalar ratio,
+ const label n
+)
+{
+ return width / rn1r1(ratio, n);
+}
+
+
+static scalar fit_slope
+(
+ const scalar left,
+ const scalar width,
+ const label n,
+ scalar& r
+)
+{
+ const scalar wol = width / left;
+
+ // Initial value
+ r = (left < width / n) ? 1.01 : 0.99;
+
+ for (int nIter = 0; nIter < 25; ++nIter)
+ {
+ scalar rm1 = r - 1.0;
+ scalar rn = pow(r, n);
+ scalar f = (rn - 1.0) * r / rm1 - wol;
+ scalar fprime = ((n * rm1 - 1) * rn + 1) / sqr(rm1);
+
+ scalar new_r = r - f / fprime;
+ scalar delta = mag(new_r - r);
+
+ r = 0.5 * (r + new_r);
+
+ // InfoErr
+ // << "New r: " << new_r << ' '
+ // << f << ' ' << fprime << ' ' << delta << ' ' << nIter << nl;
+
+ if (delta <= 1e-3)
+ {
+ break;
+ }
+ }
+
+ return end_val(width, 1.0/r, n);
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::PDRblock::gridControl
+Foam::PDRfitMeshScan::calcGridControl
+(
+ const scalar cellWidth,
+ const scalar max_zone,
+ const PDRfitMeshParams& pars
+)
+{
+ // Prune out small areas
+ {
+ const scalar minArea = sqr(cellWidth) * pars.minAreaRatio;
+ if (verbose())
+ {
+ Info<< "Checking planes. Min-area: " << minArea << nl;
+ }
+
+ for (scalar& areaValue : totalArea_)
+ {
+ if (areaValue < minArea)
+ {
+ areaValue = 0;
+ }
+ }
+
+ for (scalar& areaValue : totalArea2_)
+ {
+ if (areaValue < minArea)
+ {
+ areaValue = 0;
+ }
+ }
+ }
+
+
+ // Initially only the two outer boundaries and absolute limits
+
+ DynamicList initialCoord(nsteps_+1);
+
+ initialCoord.resize(4);
+ initialCoord[0] = -GREAT;
+ initialCoord[1] = Foam::max(limits_.min(), hard_min_);
+ initialCoord[2] = limits_.max();
+ initialCoord[3] = GREAT;
+
+
+ // Select largest areas for fitting.
+ // - done manually (instead of via sort) to handle full/half-steps
+ // - find the largest area, add that to the list and remove its
+ // area (and immediate surroundings) from the input.
+
+ // We have the "full steps overlapping so that we do not miss
+ // a group of adjacent faces that are split over a boundary.
+ // When we have identified a plane, we need to zero that
+ // step, and also the overlapping ones to avoid
+ // double-counting. when we have done that twice, we might
+ // have completely removed the half-step that was between
+ // them. That is why we also have the areas stored in
+ // half-steps, so that is not lost.
+
+ while (initialCoord.size()-4 < totalArea_.size())
+ {
+ label n_big = 0;
+ scalar largest = 0;
+
+ // Find largest area, using "real" scanned areas
+ // thus checking the stored even locations.
+
+ for (label ii=0; ii < 2*nsteps_; ++ii)
+ {
+ if (largest < totalArea_[ii])
+ {
+ largest = totalArea_[ii];
+ n_big = ii;
+ }
+ if (largest < totalArea2_[ii])
+ {
+ largest = totalArea2_[ii];
+ n_big = -ii-1;
+ }
+ }
+
+ if (mag(largest) < SMALL)
+ {
+ break;
+ }
+
+ if (n_big >= 0)
+ {
+ // Full step
+
+ const scalar pos = averagePosition(n_big);
+
+ // Remove from future checks
+ totalArea_[n_big] = 0;
+
+ if (pos > hard_min_ + 0.4 * cellWidth)
+ {
+ // Only accept if not close to the ground
+
+ initialCoord.append(pos);
+
+ // Remove this and adjacent areas, which overlap this one
+ if (n_big > 1)
+ {
+ totalArea_[n_big-1] = 0;
+ }
+ if (n_big < 2*nsteps_ - 1)
+ {
+ totalArea_[n_big+1] = 0;
+ }
+
+ // Remove half-steps contained by this step
+ totalArea2_[n_big] = 0;
+ totalArea2_[n_big+1] = 0;
+ }
+ }
+ else
+ {
+ // Half step
+ n_big = -n_big-1;
+
+ const scalar pos = averageHalfPosition(n_big);
+
+ // Remove from future checks
+ totalArea2_[n_big] = 0;
+
+ // On the first passes, this area will be included in full steps,
+ // but it may eventually only remain in the half-step,
+ // which could be the next largest
+
+ if (pos > hard_min_ + 0.4 * cellWidth)
+ {
+ // Only accept if not close to the ground
+
+ initialCoord.append(pos);
+ }
+ }
+ }
+
+ // Now sort the positions
+ Foam::sort(initialCoord);
+
+ // Avoid small zones near the begin/end positions
+
+ if
+ (
+ (initialCoord[2] - initialCoord[1])
+ < (0.4 * cellWidth)
+ )
+ {
+ // Remove [1] if too close to [2]
+ initialCoord.remove(1);
+ }
+
+ if
+ (
+ (
+ initialCoord[initialCoord.size()-2]
+ - initialCoord[initialCoord.size()-3]
+ )
+ < (0.4 * cellWidth)
+ )
+ {
+ // Remove [N-3] if too close to [N-2]
+ initialCoord.remove(initialCoord.size()-3);
+ }
+
+
+ // This is somewhat questionable (horrible mix of C and C++ addressing).
+ //
+ // If we have specified a max_zone (max number of cells per segment)
+ // we will actually generate more divisions than originally
+ // anticipated.
+
+ // So over-dimension arrays by a larger factor (eg, 20-100)
+ // which is presumably good enough.
+
+ // FIXME: Needs revisiting (2020-12-16)
+
+ label nCoords = initialCoord.size();
+
+ scalarList coord(50*initialCoord.size(), Zero);
+ SubList(coord, initialCoord.size()) = initialCoord;
+
+ scalarList width(coord.size(), Zero);
+ labelList nDivs(width.size(), Zero);
+
+ //
+ // Set up initially uniform zones
+ //
+ for (label ii=1; ii < (nCoords-2); ++ii)
+ {
+ width[ii] = coord[ii+1] - coord[ii];
+ nDivs[ii] = width[ii] / stepWidth_ * pars.areaWidthFactor + 0.5;
+
+ while (nDivs[ii] < 1)
+ {
+ // Positions too close - merge them
+ --nCoords;
+
+ coord[ii] = 0.5 * (coord[ii+1] + coord[ii]);
+ width[ii-1] = coord[ii] - coord[ii-1];
+ nDivs[ii-1] =
+ width[ii-1] / stepWidth_ * pars.areaWidthFactor + 0.5;
+
+ if (ii < (nCoords - 2))
+ {
+ for (label iii=ii+1; iii < nCoords-1; ++iii)
+ {
+ coord[iii] = coord[iii+1];
+ }
+ width[ii] = coord[ii+1] - coord[ii];
+
+ nDivs[ii] =
+ width[ii] / stepWidth_ * pars.areaWidthFactor + 0.5;
+ }
+ else
+ {
+ nDivs[ii] = 1; // Dummy value to terminate while
+ }
+ }
+ }
+
+
+ // Initially all steps have equal width (ratio == 1)
+
+ scalarList first(width.size(), Zero);
+
+ for (label ii=1; ii < (nCoords-2); ++ii)
+ {
+ first[ii] = width[ii] / nDivs[ii];
+ }
+
+ scalarList last(first);
+ scalarList ratio(width.size(), scalar(1));
+
+ // For diagnostics
+ charList marker(width.size(), char('x'));
+ marker[0] = '0';
+
+
+ // Now adjust the inter-cell ratios to reduce cell-size changes
+ // between zones
+
+ if (nCoords > 5)
+ {
+ // The outer zones can be adjusted to fit the adjacent steps
+
+ last[1] = first[2];
+ first[nCoords-3] = last[nCoords-4];
+
+ // We adjust the cell size ratio in each zone so that the
+ // cell sizes at the ends of the zones fit as well as
+ // possible with their neighbours. Work forward through the
+ // zones and then repeat a few times so that the adjustment
+ // settles down.
+
+ for (int nIter = 0; nIter < 2; ++nIter)
+ {
+ for (label ii=1; ii < (nCoords-2); ++ii)
+ {
+ if (nDivs[ii] == 1)
+ {
+ // If only one step, no grading is possible
+
+ marker[ii] = '1'; // For Diagnostics
+ continue;
+ }
+ marker[ii] = 'L'; // For Diagnostics
+
+ // Try making the in-zone steps equal to the step at
+ // the left end to the last of the previous zone.
+
+ // If the step at the right is less than this, then
+ // good... (fF this is the penultimate zone, we do
+ // not need to worry about the step on the right
+ // because the extent of the last zone can be
+ // adjusted later to fit.)
+
+ scalar r_fit = 1;
+
+ if
+ (
+ (
+ (ii > 1)
+ && mag
+ (
+ log
+ (
+ fit_slope
+ (
+ last[ii-1],
+ width[ii],
+ nDivs[ii],
+ r_fit
+ )
+ / first[ii+1]
+ )
+ )
+ < mag(log(r_fit))
+ )
+ || (ii == nCoords-3)
+ )
+ {
+ ratio[ii] = r_fit;
+ }
+ else
+ {
+ marker[ii] = 'R';
+
+ // otherwise try making the in-zone steps equal to
+ // the step at the right end to the last of the
+ // previous zone. If the step at the right is less
+ // thhan this, then good... (If this is the second
+ // zone, we do not need to worry about the step on
+ // the left because the extent of the first zone
+ // can be adjusted later to fit.)
+
+ if
+ (
+ mag
+ (
+ log
+ (
+ fit_slope
+ (
+ first[ii+1],
+ width[ii],
+ nDivs[ii],
+ r_fit
+ )
+ / last[ii-1]
+ )
+ )
+ < mag(log(r_fit))
+
+ || (ii == 1)
+ )
+ {
+ ratio[ii] = 1.0 / r_fit;
+ }
+ else
+ {
+ ratio[ii] =
+ pow(first[ii+1] / last[ii-1], 1.0/(nDivs[ii]-1));
+
+ marker[ii] = 'M';
+ }
+ }
+
+ first[ii] = end_val(width[ii], ratio[ii], nDivs[ii]);
+ last[ii] = end_val(width[ii], 1.0/ratio[ii], nDivs[ii]);
+ }
+ };
+
+
+ // Adjust the ratio and the width in the first zone to
+ // grade from cellWidth at the outside to the first cell of
+ // the next zone. Increase the number of steps to keep edge
+ // at least two (nEdgeLayers) cell widths from min.
+ }
+
+
+ if (coord[1] > hard_min_ + first[1] / pars.outerRatio)
+ {
+ nDivs[0] = pars.nEdgeLayers;
+ ratio[0] = 1.0 / pars.outerRatio;
+ coord[0] = coord[1] - first[1] * rn1r1(pars.outerRatio, nDivs[0]);
+
+ // Reduce the number of steps if we have extended below hard_min_
+ while (coord[0] < hard_min_ && nDivs[0] > 1)
+ {
+ --nDivs[0];
+ coord[0] = coord[1] - first[1] * rn1r1(pars.outerRatio, nDivs[0]);
+ }
+
+ if (nDivs[0] < label(pars.nEdgeLayers))
+ {
+ // We had to adjust for hard_min_, so now fit exaxtly
+ coord[0] = hard_min_;
+ }
+ }
+ else
+ {
+ // No room for outer zone. Remove it.
+ --nCoords;
+
+ for (label ii = 0; ii < nCoords; ++ii)
+ {
+ coord[ii] = coord[ii+1];
+ nDivs[ii] = nDivs[ii+1];
+ ratio[ii] = ratio[ii+1];
+ first[ii] = first[ii+1];
+ last[ii] = last[ii+1];
+ }
+ coord[0] = hard_min_;
+ }
+
+ width[0] = coord[1] - coord[0];
+ first[0] = end_val(width[0],ratio[0],nDivs[0]);
+ last[0] = end_val(width[0],1.0/ratio[0],nDivs[0]);
+
+
+ // Now do the upper edge zone
+ {
+ const label np2 = nCoords - 2;
+
+ ratio[np2] = pars.outerRatio;
+ nDivs[np2] = pars.nEdgeLayers;
+
+ coord[np2+1] =
+ coord[np2] + last[np2-1] * rn1r1(pars.outerRatio, nDivs[np2]);
+
+ width[np2] = coord[np2+1] - coord[np2];
+
+ first[np2] = last[np2-1];
+ last[np2] = first[np2] * pow(pars.outerRatio, nDivs[np2] - 1);
+ }
+
+
+ // If we have a zone along the length of the plant that is much
+ // longer than the width or height, then makePDRMeshBlocks does
+ // not produce a very good outer boundary shape.
+ // So here we divide the zone into several zones of roughly equal width.
+ // A bit commplicated because of the increasing or decreasing step
+ // sizes across the zone.
+
+ if (max_zone > 0)
+ {
+ for (label ii = 0; ii < (nCoords-1); ++ii)
+ {
+ if (width[ii] > max_zone && nDivs[ii] > 1)
+ {
+ // No. of extra zones
+ const label nExtra = width[ii] / max_zone - 1;
+
+ // Make space for extra zones
+ for (label ij = nCoords-1; ij > ii; --ij)
+ {
+ width[ij+nExtra] = width[ij];
+ ratio[ij+nExtra] = ratio[ij];
+ first[ij+nExtra] = first[ij];
+ last[ij+nExtra] = last[ij];
+ nDivs[ij+nExtra] = nDivs[ij];
+ coord[ij+nExtra] = coord[ij];
+ }
+ nCoords += nExtra;
+
+ const scalar subWidth = width[ii] / (nExtra + 1);
+ scalar bdy = coord[ii];
+
+ last[ii+nExtra] = last[ii];
+ nDivs[ii+nExtra] = nDivs[ii]; // will be decremented
+ width[ii+nExtra] = width[ii]; // will be decremented
+
+ // Current position
+ scalar here = coord[ii];
+ scalar step = first[ii];
+
+ // Look for cell boundaries close to where we want
+ // to put the zone boundaries
+ for (label ij = ii; ij < ii+nExtra; ++ij)
+ {
+ label ist = 0;
+ bdy += subWidth;
+ do
+ {
+ here += step;
+ step *= ratio[ii];
+ ++ist;
+ } while ((here + 0.5 * step) < bdy);
+
+ // Found a cell boundary at 'here' that is close to bdy
+
+ // Create this sub-zone
+ last[ij] = step / ratio[ii];
+ first[ij+1] = step;
+ coord[ij+1] = here;
+ width[ij] = coord[ij+1] - coord[ij];
+ nDivs[ij] = ist;
+ ratio[ij+1] = ratio[ij];
+
+ // Decrement what remains for the last sub-zone
+ nDivs[ii+nExtra] -= ist;
+ width[ii+nExtra] -= width[ij];
+ }
+
+ ii += nExtra;
+ }
+ }
+ }
+
+
+ if (verbose())
+ {
+ printf("Zone Cell widths Ratios\n");
+ printf("start first last left mid right\n");
+ for (label ii = 0; ii < (nCoords-1); ++ii)
+ {
+ printf
+ (
+ "%6.2f %6.2f %6.2f %c %6.2f %6.2f %6.2f\n",
+ coord[ii], first[ii], last[ii],
+ marker[ii],
+ first[ii]/last[ii-1], ratio[ii], first[ii+1]/last[ii]
+ );
+ }
+ printf("%6.2f\n", coord[nCoords-1]);
+ }
+
+
+ // Copy back into a gridControl form
+
+ PDRblock::gridControl grid;
+
+ grid.resize(nCoords);
+ for (label i = 0; i < nCoords; ++i)
+ {
+ grid[i] = coord[i];
+ }
+ for (label i = 0; i < nCoords-1; ++i)
+ {
+ grid.divisions_[i] = nDivs[i];
+ }
+ for (label i = 0; i < nCoords-1; ++i)
+ {
+ grid.expansion_[i] = ratio[i];
+ }
+
+ return grid;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.H b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.H
new file mode 100644
index 0000000000..3965242c05
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.H
@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of 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 .
+
+Class
+ Foam::PDRfitMeshScan
+
+Description
+ Scanning of obstacles in a single direction for PDRfitMesh
+
+SourceFiles
+ PDRfitMeshScan.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRfitMeshScan_H
+#define PDRfitMeshScan_H
+
+#include "scalarList.H"
+#include "DynamicList.H"
+#include "MinMax.H"
+#include "PDRblock.H"
+#include "PDRfitMeshParams.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class PDRfitMeshScan Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRfitMeshScan
+{
+ // Private Data
+
+ //- Lower and upper limits
+ scalarMinMax limits_;
+
+ //- Absolute lower limit (eg, for the ground plane)
+ scalar hard_min_ = -GREAT;
+
+ //- Number of steps
+ label nsteps_ = 0;
+
+ //- Step width = limits / nsteps
+ scalar stepWidth_ = 0;
+
+ //- Minimum per cell face area. Defined via PDRfitMeshParams
+ scalar minFaceArea_ = 0;
+
+ //- Area-weighted position
+ DynamicList weightedPos_;
+
+ //- Area-weights for position
+ DynamicList totalArea_;
+
+ //- Area-weighted half-cell position
+ DynamicList weightedPos2_;
+
+ //- Area-weights for half-cell position
+ DynamicList totalArea2_;
+
+
+ // Private Member Functions
+
+ //- Area-averaged position at index
+ inline scalar averagePosition(const label i) const
+ {
+ return
+ (
+ totalArea_[i] < VSMALL
+ ? 0
+ : (weightedPos_[i] / totalArea_[i])
+ );
+ }
+
+ //- Area-averaged half-position at index
+ inline scalar averageHalfPosition(const label i) const
+ {
+ return
+ (
+ totalArea2_[i] < VSMALL
+ ? 0
+ : (weightedPos2_[i] / totalArea2_[i])
+ );
+ }
+
+
+ //- Add coordinate and weight. Already sanity checked
+ void addCoord(const scalar coord, const scalar weight);
+
+
+public:
+
+ // Public Data
+
+ //- Output verbosity
+ static bool verbose_;
+
+ //- Use gridControl as per PDRblock, with EXPAND_RELATIVE
+ typedef PDRblock::gridControl gridControl;
+
+
+ // Constructors
+
+ //- Default construct
+ PDRfitMeshScan()
+ :
+ limits_(GREAT, -GREAT),
+ hard_min_(-GREAT),
+ nsteps_(0),
+ stepWidth_(0)
+ {}
+
+
+ // Static Member Functions
+
+ //- Get verbosity
+ static inline bool verbose() noexcept
+ {
+ return verbose_;
+ }
+
+ //- Set verbosity
+ static inline bool verbose(const bool val) noexcept
+ {
+ bool old(verbose_);
+ verbose_ = val;
+ return old;
+ }
+
+ //- The expansion name
+ static word expansionName();
+
+
+ // Member Functions
+
+ //- The min/max limits
+ const scalarMinMax& limits() const noexcept
+ {
+ return limits_;
+ }
+
+ //- The min/max limits
+ scalarMinMax& limits() noexcept
+ {
+ return limits_;
+ }
+
+ //- The step-width
+ scalar stepWidth() const noexcept
+ {
+ return stepWidth_;
+ }
+
+ //- Set the hard-min
+ void hard_min(const scalar val) noexcept
+ {
+ hard_min_ = val;
+ }
+
+
+ //- Reset to initial state
+ void reset();
+
+ //- Adjust limits to accommodate the specified point coordinates
+ void adjustLimits(const scalar point0, const scalar point1);
+
+ //- Define number of steps (and step-size) according to the
+ //- current limits and the specified parameters
+ void resize(const scalar cellWidth, const PDRfitMeshParams& pars);
+
+ void addArea
+ (
+ const scalar weight,
+ const scalar point0,
+ const scalar point1
+ );
+
+ //- Calculate equivalent grid control from the scanned planes
+ PDRblock::gridControl
+ calcGridControl
+ (
+ const scalar cellWidth,
+ const scalar max_zone,
+ const PDRfitMeshParams& pars
+ );
+
+ void print(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.C b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.C
new file mode 100644
index 0000000000..ece39777f3
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.C
@@ -0,0 +1,464 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of 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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRfitMeshScans.H"
+#include "PDRobstacle.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::PDRfitMeshScans::prepare
+(
+ const UList& obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+)
+{
+ reset();
+
+ // Scan for obstacle limits
+ scanLimits(obstacles);
+
+ // Clip with optional user bounds
+
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ scalarMinMax& limits = (*this)[cmpt].limits();
+
+ if (fitParams.minBounds[cmpt].has_value())
+ {
+ limits.min() = fitParams.minBounds[cmpt].value();
+ }
+ if (fitParams.maxBounds[cmpt].has_value())
+ {
+ limits.max() = fitParams.maxBounds[cmpt].value();
+ }
+ }
+
+ // Get obstacle locations and subgrid length (if any)
+ scanAreas(obstacles, fitParams, cellWidth);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::scalar Foam::PDRfitMeshScans::volume() const
+{
+ scalar vol = 1;
+
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ const scalarMinMax& limits = (*this)[cmpt].limits();
+
+ if (limits.valid())
+ {
+ vol *= limits.mag();
+ }
+ else
+ {
+ return 0;
+ }
+ }
+
+ return vol;
+}
+
+
+Foam::scalar Foam::PDRfitMeshScans::cellVolume() const
+{
+ scalar vol = 1;
+
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ vol *= (*this)[cmpt].stepWidth();
+ }
+
+ return vol;
+}
+
+
+void Foam::PDRfitMeshScans::reset()
+{
+ subgridLen_ = 0;
+
+ for (PDRfitMeshScan& pln : *this)
+ {
+ pln.reset();
+ }
+}
+
+
+void Foam::PDRfitMeshScans::resize
+(
+ const scalar cellWidth,
+ const PDRfitMeshParams& pars
+)
+{
+ for (PDRfitMeshScan& pln : *this)
+ {
+ pln.resize(cellWidth, pars);
+ }
+}
+
+
+
+
+void Foam::PDRfitMeshScans::print(Ostream& os) const
+{
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ os << vector::componentNames[cmpt] << ' ';
+ (*this)[cmpt].print(os);
+ }
+}
+
+
+void Foam::PDRfitMeshScans::scanLimits
+(
+ const UList& obstacles
+)
+{
+ for (const PDRobstacle& obs : obstacles)
+ {
+ switch (obs.typeId)
+ {
+ case PDRobstacle::CYLINDER:
+ case PDRobstacle::DIAG_BEAM:
+ {
+ (*this)[obs.orient].adjustLimits
+ (
+ obs.pt[obs.orient],
+ obs.pt[obs.orient] + obs.len()
+ );
+ break;
+ }
+
+ case PDRobstacle::CUBOID_1:
+ case PDRobstacle::LOUVRE_BLOWOFF:
+ case PDRobstacle::CUBOID:
+ case PDRobstacle::WALL_BEAM:
+ case PDRobstacle::GRATING:
+ case PDRobstacle::RECT_PATCH:
+ {
+ (*this).x().adjustLimits(obs.x(), obs.x() + obs.span.x());
+ (*this).y().adjustLimits(obs.y(), obs.y() + obs.span.y());
+ (*this).z().adjustLimits(obs.z(), obs.z() + obs.span.z());
+ break;
+ }
+ }
+ }
+}
+
+
+Foam::scalar
+Foam::PDRfitMeshScans::scanAreas
+(
+ const UList& obstacles,
+ const scalar minSubgridLen
+)
+{
+ scalar subgridLen = 0;
+
+ const bool checkSubgrid = (minSubgridLen > 0);
+
+ if (verbose())
+ {
+ Info<< "Scan " << obstacles.size()
+ << " obstacles. Check for subgrid < " << minSubgridLen
+ << nl;
+ }
+
+ // When sorting flat elements, convert vector -> List
+ List gridSpan(3);
+
+ for (const PDRobstacle& obs : obstacles)
+ {
+ switch (obs.typeId)
+ {
+ case PDRobstacle::CYLINDER:
+ {
+ // #ifdef FULLDEBUG
+ // Info<< "Add " << obs.info() << nl;
+ // #endif
+
+ (*this)[obs.orient].addArea
+ (
+ sqr(0.5*obs.dia()),
+ obs.pt[obs.orient],
+ obs.pt[obs.orient] + obs.len()
+ );
+
+ // Store total length of subgrid obstcales
+ if (checkSubgrid && obs.dia() < minSubgridLen)
+ {
+ subgridLen += obs.len();
+ }
+ break;
+ }
+
+ case PDRobstacle::DIAG_BEAM:
+ {
+ // #ifdef FULLDEBUG
+ // Info<< "Add " << obs.info() << nl;
+ // #endif
+
+ (*this)[obs.orient].addArea
+ (
+ (obs.wa * obs.wb),
+ obs.pt[obs.orient],
+ obs.pt[obs.orient] + obs.len()
+ );
+
+ // Store total length of subgrid obstcales
+ if (checkSubgrid && Foam::max(obs.wa, obs.wb) < minSubgridLen)
+ {
+ subgridLen += obs.len();
+ }
+ break;
+ }
+
+ case PDRobstacle::CUBOID_1:
+ case PDRobstacle::LOUVRE_BLOWOFF:
+ case PDRobstacle::CUBOID:
+ case PDRobstacle::WALL_BEAM:
+ case PDRobstacle::GRATING:
+ case PDRobstacle::RECT_PATCH:
+ {
+ // #ifdef FULLDEBUG
+ // Info<< "Add " << obs.info() << nl;
+ // #endif
+
+ // Can be lazy here, addArea() filters out zero areas
+
+ (*this)[vector::X].addArea
+ (
+ (obs.span.y() * obs.span.z()),
+ obs.x(),
+ obs.x() + obs.span.x()
+ );
+
+ (*this)[vector::Y].addArea
+ (
+ (obs.span.z() * obs.span.x()),
+ obs.y(),
+ obs.y() + obs.span.y()
+ );
+
+ (*this)[vector::Z].addArea
+ (
+ (obs.span.x() * obs.span.y()),
+ obs.z(),
+ obs.z() + obs.span.z()
+ );
+
+ if (checkSubgrid)
+ {
+ gridSpan[0] = obs.span.x();
+ gridSpan[1] = obs.span.y();
+ gridSpan[2] = obs.span.z();
+ Foam::sort(gridSpan);
+
+ // Ignore zero dimension when considering subgrid
+
+ // Store total length of subgrid obstcales
+ if (gridSpan[1] < minSubgridLen)
+ {
+ subgridLen += gridSpan.last();
+ }
+ }
+ break;
+ }
+ }
+ }
+
+ return subgridLen;
+}
+
+
+void Foam::PDRfitMeshScans::scanAreas
+(
+ const UList& obstacles
+)
+{
+ (void)scanAreas(obstacles, -GREAT);
+}
+
+
+void Foam::PDRfitMeshScans::scanAreas
+(
+ const UList& obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+)
+{
+ resize(mag(cellWidth), fitParams);
+
+ const scalar minSubgridLen =
+ (cellWidth < 0 ? (mag(cellWidth) * fitParams.widthFactor) : 0);
+
+ subgridLen_ = scanAreas(obstacles, minSubgridLen);
+}
+
+
+Foam::Vector
+Foam::PDRfitMeshScans::calcGriding
+(
+ const UList& obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+)
+{
+ prepare(obstacles, fitParams, cellWidth);
+
+ const scalar innerVol = volume();
+
+ // Set hard lower limit at the ground
+
+ z().hard_min(fitParams.ground);
+
+ if (z().limits().min() < fitParams.ground)
+ {
+ z().limits().min() = fitParams.ground;
+ }
+
+ if (verbose())
+ {
+ print(Info);
+ Info<< "cellWidth = " << cellWidth << nl;
+ Info<< "inner volume: " << innerVol << nl;
+ Info<< "subgrid length: " << subgridLen_ << nl;
+ }
+
+
+ // Read in obstacles and record face planes Also record length of
+ // subgrid obstacles so that we can optimise cell size to get
+ // desired average no, of obstacles oper cell. Determinuing which
+ // obstacles are subgrid is iself dependent on the cell size.
+ // Therefore we iterate (if cell_width is -ve) If cell_width is
+ // +ve just use the user-supplied value.
+
+ {
+ // The cell-widths during optimisation
+ scalar prev_cw = 0, cw = 0;
+
+ const bool optimiseWidth = (cellWidth < 0);
+
+ for (label nIter = 0; nIter < fitParams.maxIterations; ++nIter)
+ {
+ scanAreas(obstacles, fitParams, cellWidth);
+
+ if (cellWidth < 0)
+ {
+ constexpr scalar relax = 0.7;
+
+ prev_cw = cw;
+
+ if (subgridLen_ < cellWidth)
+ {
+ FatalErrorInFunction
+ << "No sub-grid obstacles found" << endl
+ << exit(FatalError);
+ }
+
+ // Optimise to average subgrid obstacles per cell
+ cw = sqrt(fitParams.obsPerCell * innerVol / subgridLen_);
+ cw = Foam::min(cw, fitParams.maxCellWidth);
+
+ if (cw > cellWidth * fitParams.maxWidthEstimate)
+ {
+ Warning
+ << "Calculated cell width more than"
+ " maxWidthEstimate x estimate." << nl
+ << "Too few sub-grid obstacles?"
+ << nl;
+ }
+
+ Info<< "Current cellwidth: " << cw << nl;
+
+ scalar ratio = cw / cellWidth;
+
+ if (ratio < 1)
+ {
+ ratio = 1/ratio;
+ }
+
+ if (ratio < fitParams.maxWidthRatio)
+ {
+ cellWidth = -cellWidth;
+ }
+ else
+ {
+ cellWidth = relax * (-cw) + (1.0 - relax) * cellWidth;
+ }
+ Info<< "Cell width: " << cellWidth << nl;
+ }
+
+ if (cellWidth > 0)
+ {
+ break;
+ }
+ }
+
+ if (cellWidth < 0)
+ {
+ cellWidth = 0.5 * (cw + prev_cw);
+ }
+
+ if (optimiseWidth)
+ {
+ Info<< nl << "Final cell width: " << cellWidth << nl << nl;
+ }
+ }
+
+
+ // Now we fit the mesh to the planes and write out the result
+
+ resize(cellWidth, fitParams);
+ cellWidth = cbrt(cellVolume()) / fitParams.areaWidthFactor;
+
+ scanAreas(obstacles, fitParams, cellWidth);
+
+ const scalar max_zone =
+ (
+ (z().limits().mag() + fitParams.nEdgeLayers * cellWidth)
+ * fitParams.maxZoneToHeight
+ );
+
+
+ Vector griding;
+
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ griding[cmpt] =
+ (*this)[cmpt].calcGridControl(cellWidth, max_zone, fitParams);
+ }
+
+ return griding;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.H b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.H
new file mode 100644
index 0000000000..ebf974bcd6
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.H
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of 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 .
+
+Class
+ Foam::PDRfitMeshScans
+
+Description
+ Scanning of obstacles in a multiple directions for PDRfitMesh
+
+SourceFiles
+ PDRfitMeshScans.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRfitMeshScans_H
+#define PDRfitMeshScans_H
+
+#include "PDRfitMeshScan.H"
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class PDRobstacle;
+
+/*---------------------------------------------------------------------------*\
+ Class PDRfitMeshScans Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRfitMeshScans
+:
+ public Vector
+{
+ // Private Data
+
+ //- The evaluated sub-grid length
+ scalar subgridLen_ = 0;
+
+
+ // Private Member Functions
+
+ //- Prepare limits, scan areas
+ void prepare
+ (
+ const UList& obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+ );
+
+public:
+
+ // Default Constructors
+
+
+ // Member Functions
+
+ //- Verbosity
+ static inline bool verbose() noexcept
+ {
+ return PDRfitMeshScan::verbose();
+ }
+
+ //- Calculate the grid controls for the given obstacles
+ //- and parameters
+ Vector
+ calcGriding
+ (
+ const UList& obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+ );
+
+
+ // Low-level functions
+
+ //- The volume of the limits
+ scalar volume() const;
+
+ //- The volume of single cell of the respective step-width
+ scalar cellVolume() const;
+
+ //- Reset to initial state
+ void reset();
+
+ //- Define number of steps (and step-size) according to the
+ //- current limits and the specified parameters
+ void resize(const scalar cellWidth, const PDRfitMeshParams& pars);
+
+ //- Adjust directional limits to accommodate obstacles
+ void scanLimits(const UList& obstacles);
+
+ //- Populate with areas/positions of the obstacles,
+ //- and check for sub-grid obstacles
+ // Sets the internal
+ void scanAreas
+ (
+ const UList& obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+ );
+
+ //- Populate with areas/positions of the obstacles,
+ //- and check for sub-grid obstacles
+ //
+ // \return total subgrid lengths
+ scalar scanAreas
+ (
+ const UList& obstacles,
+ const scalar minSubgridLen
+ );
+
+ //- Populate with areas/positions of the obstacles,
+ //- without checks for sub-grid obstacles
+ void scanAreas(const UList& obstacles);
+
+ //- Print information
+ void print(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C b/applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C
index d5d8b6c02a..a753f57da7 100644
--- a/applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C
+++ b/applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C
@@ -82,6 +82,8 @@ int main(int argc, char* argv[])
#include "setRootCase.H"
#include "createTime.H"
+ const bool dryrun = args.found("dry-run");
+
const word dictName("PDRsetFieldsDict");
#include "setSystemRunTimeDictionaryIO.H"
diff --git a/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRfitMeshDict b/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRfitMeshDict
new file mode 100644
index 0000000000..fb634d2b92
--- /dev/null
+++ b/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRfitMeshDict
@@ -0,0 +1,105 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: v2012 |
+| \\ / A nd | Website: www.openfoam.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object PDRfitMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Settings for PDRfitMesh
+
+cellWidth 0.22;
+
+// Tuning parameters (defaults) for PDRfitMesh
+
+bounds
+{
+// xmin -1e16;
+// xmax 1e16;
+//
+// ymin -1e16;
+// ymax 1e16;
+//
+// zmin -1e16;
+// zmax 1e16;
+}
+
+
+findPlanes
+{
+ // A face smaller than this multiple of cell face area will not be
+ // recorded
+ minSigFaceArea 5;
+
+ // Only fit a plane if the face area at the coordinate is at least
+ // this times the cell face area
+ minAreaRatio 20.0;
+
+ // Size of search zones for face areas will be this * the cell width.
+ // So faces closer than this zone size might be grouped together
+ areasetWidthFactor 0.7;
+
+ // Very long zones produce bad outer boundary shape from
+ // makePDRmeshBlocks, so we subdivide a zone if its length is greater
+ // than this * the height of the area of cuboid vells
+ //
+ maxRatioZoneToH 2.0 ;
+}
+
+
+/****************** Choosing cell width ***************/
+
+// Note "cellWidth" and "cellWidthFactor" are read from PDRsetFieldsDict
+// cellWidth
+// {
+// }
+
+// Note "cellWidth" and "cellWidthFactor" are read from PDRsetFieldsDict
+grid
+{
+ auto true;
+ width 0;
+ max great;
+}
+
+
+// Width of obstacles must be less than this * cell width to added
+// into subgrid length
+widthFactor 1.0;
+
+// Optimal average number of obstacles per cell
+obsPerCell 2.0;
+
+// Minimum number of cells in any direction
+minStepsPerDim 10;
+
+// Do not try a cell width more than this times the initial estimate
+maxCwToEst 5.0;
+
+// Assume converged if optimised width changes by less than this
+maxWidthRatio 1.2;
+
+// Maximum iterations in optimizing cell width
+maxPasses 5;
+
+
+/****************** Defining outer expanding region ***************/
+
+// Number of rectangular cell layers each side of the congested region. (double, not int)
+nEdgeLayers 5.0;
+
+// Ratio of the outer extension distance to the average radius of the congested region
+ratioOuterToCongRad 20.0;
+
+// Cell size ratio - should be same as that hard-wired in makePDRmeshBlocks
+outerRatio 1.2;
+
+
+// ************************************************************************* //
diff --git a/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRsetFieldsDict b/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRsetFieldsDict
index 1f101e0986..d6b212b7ad 100644
--- a/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRsetFieldsDict
+++ b/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRsetFieldsDict
@@ -35,6 +35,24 @@ cellWidth 0.22;
// Optional
cellWidthFactor 1.0;
+// Here or in PDRfitMeshDict. Largely pass-through to PDRblockMesh
+
+// Copy for PDRblockMeshDict
+outer
+{
+ type sphere;
+ onGround true;
+ expansion relative;
+
+ ratios 1.1;
+
+ size 3;
+ nCells 10;
+
+ // For PDRfitMesh only:
+ zmin 0;
+}
+
// ------------------
// Advanced