diff --git a/applications/utilities/preProcessing/PDR/Allwclean b/applications/utilities/preProcessing/PDR/Allwclean
index 9f1e7393322b44da5d6fd2253a86c4817c61fbff..4478b88ce6123d565bd97050c1a5a6bd764a67c4 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 2131eed2f8dc738a95ed2155ade3a4bf27955b29..10618728d884d7d57f5cc1a00644a55979ccf96c 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 0000000000000000000000000000000000000000..ccf95e298a1f7e8f1fbf02b86d5f782464c9b7f9
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/Make/files
@@ -0,0 +1,7 @@
+PDRfitMeshParams.C
+PDRfitMeshPlane.C
+PDRfitMeshPlanes.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 0000000000000000000000000000000000000000..64d6ba56258273f7da3b85ec0dfa36f3afa5f51d
--- /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 0000000000000000000000000000000000000000..2ef5dc2d92b77050c9924dde21119439793945cd
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMesh.C
@@ -0,0 +1,346 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 <http://www.gnu.org/licenses/>.
+
+Applications
+ PDRfitMesh
+
+Description
+ Scans extends of obstacles to establish reasonable estimates
+ for generating a PDRblockMeshDict.
+
+SourceFiles
+ PDRfitMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "IOdictionary.H"
+
+#include "PDRfitMeshParams.H"
+#include "PDRfitMeshPlanes.H"
+#include "PDRlegacy.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
+ (
+ "legacy",
+ "Force use of legacy obstacles table"
+ );
+
+ argList::addBoolOption
+ (
+ "dry-run",
+ "Read obstacles, write VTK and display output"
+ );
+
+ argList::addBoolOption
+ (
+ "print-dict",
+ "Print PDRblockMeshDict equivalent and exit"
+ );
+
+ #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);
+
+ if (args.found("legacy"))
+ {
+ pars.legacyObsSpec = true;
+ }
+
+
+ // Params for fitMesh
+ // - like setSystemRunTimeDictionaryIO
+
+ IOobject paramIO = IOobject::selectIO
+ (
+ IOobject
+ (
+ "PDRfitMeshDict",
+ runTime.system(),
+ runTime,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ ),
+ args.getOrDefault<fileName>("params", "")
+ );
+
+ if (paramIO.typeHeaderOk<IOdictionary>(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<PDRobstacle> obstacles, cylinders;
+
+ // Read in obstacles
+ if (pars.legacyObsSpec)
+ {
+ PDRobstacle::legacyReadFiles
+ (
+ pars.obsfile_dir, pars.obsfile_names,
+ boundBox::invertedBox,
+ obstacles,
+ cylinders
+ );
+ }
+ else
+ {
+ PDRobstacle::readFiles
+ (
+ pars.obsfile_dir, pars.obsfile_names,
+ boundBox::invertedBox,
+ 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 (dryrun)
+ {
+ PDRobstacle::generateVtk(casepath/"VTK", obstacles, cylinders);
+ }
+
+
+ //
+ // The fitting routines
+ //
+
+ PDRfitMeshPlane::verbose_ = args.found("verbose");
+
+ Vector<PDRblock::gridControl> griding =
+ PDRfitMeshPlanes().calcGriding
+ (
+ {obstacles, cylinders},
+ fitParams,
+ cellWidth
+ );
+
+
+ //
+ // Output
+ //
+
+ if (writable == writeType::DRY_RUN)
+ {
+ InfoErr
+ << nl
+ << "dry-run: 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<Ostream> 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",
+ PDRblock::expansionNames_[PDRblock::EXPAND_RELATIVE]
+ );
+ 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 0000000000000000000000000000000000000000..b29f67f2cb034b7cefa667a223abd37c328c576c
--- /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
+//
+// <cellWidth> and <cellWidthFactor> 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.
+edgeLayers 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 0000000000000000000000000000000000000000..a83577764d233a02d373f08bffbe27b503d9ec4d
--- /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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRfitMeshParams.H"
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::PDRfitMeshParams::readBounds(const dictionary& dict)
+{
+ scalar value;
+
+ if (dict.readIfPresent("xmin", value))
+ {
+ minBounds.x() = optionalData<scalar>(value);
+ }
+ if (dict.readIfPresent("ymin", value))
+ {
+ minBounds.y() = optionalData<scalar>(value);
+ }
+ if (dict.readIfPresent("zmin", value))
+ {
+ minBounds.z() = optionalData<scalar>(value);
+ }
+
+ if (dict.readIfPresent("xmax", value))
+ {
+ maxBounds.x() = optionalData<scalar>(value);
+ }
+ if (dict.readIfPresent("ymax", value))
+ {
+ maxBounds.y() = optionalData<scalar>(value);
+ }
+ if (dict.readIfPresent("zmax", value))
+ {
+ maxBounds.z() = optionalData<scalar>(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("edgeLayers", edgeLayers);
+
+ 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 0000000000000000000000000000000000000000..b27846e27d6092772e544efd7b4e8db2510add5b
--- /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 <http://www.gnu.org/licenses/>.
+
+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<optionalData<scalar>> 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 edgeLayers = 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/PDRfitMeshPlane.C b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlane.C
new file mode 100644
index 0000000000000000000000000000000000000000..740a50b062f5e34be0bdf51c5104b76a26618051
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlane.C
@@ -0,0 +1,831 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRfitMeshPlane.H"
+#include "PDRfitMeshParams.H"
+#include "PDRobstacle.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+bool Foam::PDRfitMeshPlane::verbose_ = false;
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::PDRfitMeshPlane::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 ranges half a 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::PDRfitMeshPlane::clear()
+// {
+// weightedPos_.clear(); totalArea_.clear();
+// weightedPos2_.clear(); totalArea2_.clear();
+// }
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::PDRfitMeshPlane::reset()
+{
+ limits.min() = GREAT;
+ limits.max() = -GREAT;
+
+ hard_min = -GREAT;
+
+ nsteps_ = 0;
+ stepWidth_ = 0;
+
+ weightedPos_.clear();
+ totalArea_.clear();
+
+ weightedPos2_.clear();
+ totalArea2_.clear();
+}
+
+
+void Foam::PDRfitMeshPlane::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 * sqr(stepWidth_);
+ }
+ 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::PDRfitMeshPlane::adjustLimits
+(
+ const scalar point0,
+ const scalar point1
+)
+{
+ limits.add(point0);
+ limits.add(point1);
+}
+
+
+void Foam::PDRfitMeshPlane::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::PDRfitMeshPlane::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(scalar r, label n)
+{
+ if (mag(r-1.0) < SMALL)
+ {
+ return n;
+ }
+
+ return (pow(r, n) - 1.0)/(r - 1.0);
+}
+
+
+inline static scalar end_val(scalar width, scalar ratio, label n)
+{
+ return width / rn1r1(ratio, n);
+}
+
+
+static scalar fit_slope(scalar left, scalar width, label n, scalar& r)
+{
+ const scalar wol = width / left;
+
+ // Initial value
+ if (left < width / n)
+ {
+ r = 1.01;
+ }
+ else
+ {
+ r = 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) / (rm1 * 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::PDRfitMeshPlane::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;
+ }
+ }
+ }
+
+
+ // Initialy only the two outer boundaries and absolute limits
+
+ DynamicList<scalar> plane(nsteps_+1);
+
+ plane.resize(4);
+ plane[0] = -GREAT;
+ plane[1] = Foam::max(limits.min(), hard_min);
+ plane[2] = limits.max();
+ plane[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.
+
+ while (plane.size() < totalArea_.size())
+ {
+ scalar largest = 0;
+ label n_big = 0;
+
+ 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;
+ }
+
+ // 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.
+
+ 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
+
+ plane.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
+
+ plane.append(pos);
+ }
+ }
+ }
+
+
+ // Now sort the positions
+ Foam::sort(plane);
+
+ if ((plane[2] - plane[1]) < 0.4 * cellWidth)
+ {
+ // Remove [1] if too close to [2]
+ plane.remove(1);
+ }
+
+ if ((plane[plane.size()-2] - plane[plane.size()-3] ) < 0.4 * cellWidth)
+ {
+ // Remove [N-3] if too close to [N-2]
+ plane.remove(plane.size()-3);
+ }
+
+
+ label n_planes = plane.size();
+
+
+ // 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 that originally
+ // anticipated.
+
+ // So over-dimension the arrays by a larger factor (eg, 100)
+ // which is generally enough.
+
+ // FIXME: Needs revisiting (2020-12-16)
+
+ scalarList width(20*n_planes);
+
+ labelList i_steps(width.size());
+
+ //
+ // Set up initially uniform zones
+ //
+ for (label ii=1; ii < (n_planes-2); ++ii)
+ {
+ width[ii] = plane[ii+1] - plane[ii];
+ i_steps[ii] = width[ii] / stepWidth_ * pars.areaWidthFactor + 0.5;
+
+ while (i_steps[ii] < 1)
+ {
+ // Planes too close - merge them
+ --n_planes;
+
+ plane[ii] = 0.5 * (plane[ii+1] + plane[ii]);
+ width[ii-1] = plane[ii] - plane[ii-1];
+
+ i_steps[ii-1] =
+ width[ii-1] / stepWidth_ * pars.areaWidthFactor + 0.5;
+
+ if (ii < (n_planes - 2))
+ {
+ for (label iii=ii+1; iii < n_planes-1; ++iii)
+ {
+ plane[iii] = plane[iii+1];
+ }
+ width[ii] = plane[ii+1] - plane[ii];
+
+ i_steps[ii] =
+ width[ii] / stepWidth_ * pars.areaWidthFactor + 0.5;
+ }
+ else
+ {
+ i_steps[ii] = 1; // Dummy value to terminate while
+ }
+ }
+ }
+
+
+ // Initially all steps have equal width (ratio == 1).
+
+ scalarList first(width.size());
+
+ for (label ii=1; ii < (n_planes-2); ++ii)
+ {
+ first[ii] = width[ii] / i_steps[ii];
+ }
+ scalarList last(first);
+
+ //
+ // Now adjust the inter-cell ratios to reduce cell-size changes
+ // between zones
+ //
+
+ scalarList ratio(width.size(), scalar(1));
+
+ // For diagnostics
+ charList marker(width.size(), char('x'));
+ marker[0] = '0';
+
+ if (n_planes > 5)
+ {
+ // The outer zones can be adjusted to fit the adjacent steps
+
+ last[1] = first[2];
+ first[n_planes-3] = last[n_planes-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 < (n_planes-2); ++ii)
+ {
+ if (i_steps[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],
+ i_steps[ii],
+ r_fit
+ )
+ / first[ii+1]
+ )
+ )
+ < mag(log(r_fit))
+ )
+ || (ii == n_planes-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],
+ i_steps[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/(i_steps[ii]-1));
+
+ marker[ii] = 'M';
+ }
+ }
+
+ first[ii] = end_val(width[ii], ratio[ii], i_steps[ii] );
+ last[ii] = end_val(width[ii], 1.0/ratio[ii], i_steps[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 (edgeLayers) cell widths from min.
+ }
+
+
+ if (plane[1] > hard_min + first[1] / pars.outerRatio)
+ {
+ i_steps[0] = pars.edgeLayers;
+ ratio[0] = 1.0 / pars.outerRatio;
+ plane[0] = plane[1] - first[1] * rn1r1( pars.outerRatio, i_steps[0]);
+
+ // Reduce the number of steps if we have extended below hard_min
+ while (plane[0] < hard_min && i_steps[0] > 1)
+ {
+ --i_steps[0];
+ plane[0] = plane[1] - first[1] * rn1r1(pars.outerRatio, i_steps[0]);
+ }
+
+ if (i_steps[0] < pars.edgeLayers)
+ {
+ // We had to adjust for hard_min, so now fit exaxtly
+ plane[0] = hard_min;
+ }
+ }
+ else
+ {
+ // No room for outer zone. Remove it.
+ --n_planes;
+
+ for (label ii = 0; ii < n_planes; ++ii)
+ {
+ plane[ii] = plane[ii+1];
+ i_steps[ii] = i_steps[ii+1];
+ ratio[ii] = ratio[ii+1];
+ first[ii] = first[ii+1];
+ last[ii] = last[ii+1];
+ }
+ plane[0] = hard_min;
+ }
+
+ width[0] = plane[1] - plane[0];
+ first[0] = end_val(width[0],ratio[0],i_steps[0]);
+ last[0] = end_val(width[0],1.0/ratio[0],i_steps[0]);
+
+
+ // Now do the upper edge zone
+ {
+ const label np2 = n_planes - 2;
+
+ ratio[np2] = pars.outerRatio;
+ i_steps[np2] = pars.edgeLayers;
+
+ plane[np2+1] =
+ plane[np2] + last[np2-1] * rn1r1(pars.outerRatio, i_steps[np2] );
+
+ width[np2] = plane[np2+1] - plane[np2];
+
+ first[np2] = last[np2-1];
+ last[np2] = first[np2] * pow(pars.outerRatio, i_steps[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 < (n_planes-1); ++ii)
+ {
+ if (width[ii] > max_zone && i_steps[ii] > 1)
+ {
+ // No. of extra zones
+ const label n_extra = width[ii] / max_zone - 1;
+
+ // Make space for extra zones
+ for (label ij = n_planes-1; ij > ii; --ij)
+ {
+ width[ij+n_extra] = width[ij];
+ ratio[ij+n_extra] = ratio[ij];
+ first[ij+n_extra] = first[ij];
+ last[ij+n_extra] = last[ij];
+ i_steps[ij+n_extra] = i_steps[ij];
+ plane[ij+n_extra] = plane[ij];
+ }
+ n_planes += n_extra;
+
+ const scalar sub_w = width[ii] / (n_extra + 1);
+ scalar bdy = plane[ii];
+
+ last[ii+n_extra] = last[ii];
+ i_steps[ii+n_extra] = i_steps[ii]; // will be decremented
+ width[ii+n_extra] = width[ii]; // will be decremented
+
+ // Current position
+ scalar coord = plane[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+n_extra; ++ij)
+ {
+ label ist = 0;
+ bdy += sub_w;
+ do
+ {
+ coord += step;
+ step *= ratio[ii];
+ ++ist;
+ } while ((coord + 0.5 * step) < bdy);
+
+ // Found a cell boundary at coord that is close to bdy
+
+ // Create this sub-zone
+ last[ij] = step / ratio[ii];
+ first[ij+1] = step;
+ plane[ij+1] = coord;
+ width[ij] = plane[ij+1] - plane[ij];
+ i_steps[ij] = ist;
+ ratio[ij+1] = ratio[ij];
+
+ // Decrement what remaons for the last zub-zone
+ i_steps[ii+n_extra] -= ist;
+ width[ii+n_extra] -= width[ij];
+ }
+
+ ii += n_extra;
+ }
+ }
+ }
+
+
+ if (verbose_)
+ {
+ printf("Zone Cell widths Ratios\n");
+ printf("start first last left mid right\n");
+ for (label ii = 0; ii < (n_planes-1); ++ii)
+ {
+ printf
+ (
+ "%6.2f %6.2f %6.2f %c %6.2f %6.2f %6.2f\n",
+ plane[ii], first[ii], last[ii],
+ marker[ii],
+ first[ii]/last[ii-1], ratio[ii], first[ii+1]/last[ii]
+ );
+ }
+ printf("%6.2f\n", plane[n_planes-1]);
+ }
+
+
+ // Copy back into a gridControl form
+
+ PDRblock::gridControl grid;
+
+ grid.resize(n_planes);
+ for (label i = 0; i < n_planes; ++i)
+ {
+ grid[i] = plane[i];
+ }
+ for (label i = 0; i < n_planes-1; ++i)
+ {
+ grid.divisions_[i] = i_steps[i];
+ }
+ for (label i = 0; i < n_planes-1; ++i)
+ {
+ grid.expansion_[i] = ratio[i];
+ }
+
+ return grid;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlane.H b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlane.H
new file mode 100644
index 0000000000000000000000000000000000000000..35cdfba0bc408852e5ac3bef2fba3216d5d578f7
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlane.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::PDRfitMeshPlane
+
+Description
+ Handling for PDRfitMesh
+
+SourceFiles
+ PDRfitMeshPlane.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRfitMeshPlane_H
+#define PDRfitMeshPlane_H
+
+#include "scalarList.H"
+#include "DynamicList.H"
+#include "MinMax.H"
+#include "PDRblock.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class PDRfitMeshParams;
+class PDRobstacle;
+
+/*---------------------------------------------------------------------------*\
+ Class PDRfitMeshPlane Declaration
+\*---------------------------------------------------------------------------*/
+
+//- A single mesh plane
+class PDRfitMeshPlane
+{
+ // Private Data
+
+ DynamicList<scalar> weightedPos_, totalArea_;
+ DynamicList<scalar> weightedPos2_, totalArea2_;
+
+
+ // Private Member Functions
+
+ //- Area-averaged position at index
+ scalar averagePosition(const label i) const
+ {
+ return
+ (
+ totalArea_[i] < VSMALL
+ ? 0
+ : (weightedPos_[i] / totalArea_[i])
+ );
+ }
+
+ //- Area-averaged half-position at index
+ 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_;
+
+ //- Reuse gridControl
+ typedef PDRblock::gridControl gridControl;
+
+
+ // Data Members
+
+ //- Lower and upper limits
+ scalarMinMax limits;
+
+ scalar hard_min;
+
+ scalar minFaceArea = 0;
+
+ label nsteps_;
+ scalar stepWidth_;
+
+
+ // Constructors
+
+ //- Default construct
+ PDRfitMeshPlane()
+ :
+ limits(GREAT, -GREAT),
+ hard_min(-GREAT),
+ nsteps_(0),
+ stepWidth_(0)
+ {}
+
+
+ // Member Functions
+
+ //- 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/PDRfitMeshPlanes.C b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlanes.C
new file mode 100644
index 0000000000000000000000000000000000000000..e6a7cec08bdb4f0d9d83004ae0a6bb7dba7ea9aa
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlanes.C
@@ -0,0 +1,466 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRfitMeshPlanes.H"
+#include "PDRfitMeshParams.H"
+#include "PDRobstacle.H"
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::PDRfitMeshPlanes::prepare
+(
+ std::initializer_list<UList<PDRobstacle>> obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+)
+{
+ this->reset();
+
+ // Scan for obstacle limits
+ for (const UList<PDRobstacle>& list : obstacles)
+ {
+ this->scanLimits(list);
+ }
+
+ // 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::PDRfitMeshPlanes::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::PDRfitMeshPlanes::cellVolume() const
+{
+ scalar vol = 1;
+
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ vol *= (*this)[cmpt].stepWidth_;
+ }
+
+ return vol;
+}
+
+
+void Foam::PDRfitMeshPlanes::reset()
+{
+ subgridLen_ = 0;
+
+ for (PDRfitMeshPlane& pln : *this)
+ {
+ pln.reset();
+ }
+}
+
+
+void Foam::PDRfitMeshPlanes::resize
+(
+ const scalar cellWidth,
+ const PDRfitMeshParams& pars
+)
+{
+ for (PDRfitMeshPlane& pln : *this)
+ {
+ pln.resize(cellWidth, pars);
+ }
+}
+
+
+
+
+void Foam::PDRfitMeshPlanes::print(Ostream& os) const
+{
+ for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+ {
+ os << vector::componentNames[cmpt] << ' ';
+ (*this)[cmpt].print(os);
+ }
+}
+
+
+void Foam::PDRfitMeshPlanes::scanLimits
+(
+ const UList<PDRobstacle>& 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::PDRfitMeshPlanes::scanAreas
+(
+ const UList<PDRobstacle>& 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<scalar> 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::PDRfitMeshPlanes::scanAreas
+(
+ const UList<PDRobstacle>& obstacles
+)
+{
+ (void)scanAreas(obstacles, -GREAT);
+}
+
+
+void Foam::PDRfitMeshPlanes::scanAreas
+(
+ std::initializer_list<UList<PDRobstacle>> obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+)
+{
+ resize(mag(cellWidth), fitParams);
+
+ const scalar minSubgridLen =
+ (cellWidth < 0 ? (mag(cellWidth) * fitParams.widthFactor) : 0);
+
+ subgridLen_ = 0;
+
+ for (const UList<PDRobstacle>& list : obstacles)
+ {
+ subgridLen_ += scanAreas(list, minSubgridLen);
+ }
+}
+
+
+Foam::Vector<Foam::PDRblock::gridControl>
+Foam::PDRfitMeshPlanes::calcGriding
+(
+ std::initializer_list<UList<PDRobstacle>> obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+)
+{
+ prepare(obstacles, fitParams, cellWidth);
+
+ scalar innerVol = this->volume();
+
+ // Hard limit at the ground
+ z().hard_min = fitParams.ground;
+
+ if (z().limits.min() < fitParams.ground)
+ {
+ z().limits.min() = fitParams.ground;
+ }
+
+ if (verbose())
+ {
+ this->print(Info);
+ Info<< "cellWidth = " << cellWidth << nl;
+ Info<< "inner volume: " << innerVol << nl;
+ Info<< "subgrid length: " << subgridLen_ << nl;
+ }
+
+
+ scalar cw = 0;
+ scalar prev_cw = 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 OBS_PER_CELL 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;
+ }
+ }
+
+
+ // 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.
+
+ if (cellWidth < 0)
+ {
+ cellWidth = 0.5 * (cw + prev_cw);
+ }
+ Info<< nl << "Final cell width: " << cellWidth << nl << nl;
+
+
+ // Now we fit the mesh to the planes and write out the result
+
+ this->resize(cellWidth, fitParams);
+ cellWidth = cbrt(this->cellVolume()) / fitParams.areaWidthFactor;
+
+ scanAreas(obstacles, fitParams, cellWidth);
+
+ const scalar max_zone =
+ (
+ (this->z().limits.mag() + fitParams.edgeLayers * cellWidth)
+ * fitParams.maxZoneToHeight
+ );
+
+
+ Vector<PDRblock::gridControl> 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/PDRfitMeshPlanes.H b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlanes.H
new file mode 100644
index 0000000000000000000000000000000000000000..26b7e8ff8e6a8906aacf5b2a76187ff882df3c59
--- /dev/null
+++ b/applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshPlanes.H
@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::PDRfitMeshPlanes
+
+Description
+ Handling for PDRfitMesh
+
+SourceFiles
+ PDRfitMeshPlanes.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRfitMeshPlanes_H
+#define PDRfitMeshPlanes_H
+
+#include "PDRfitMeshPlane.H"
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class PDRfitMeshPlanes Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRfitMeshPlanes
+:
+ public Vector<PDRfitMeshPlane>
+{
+ // Private Data
+
+ //- The evaluated sub-grid lengths
+ scalar subgridLen_ = 0;
+
+
+ // Private Member Functions
+
+ //- Use verbosity
+ static inline bool verbose()
+ {
+ return PDRfitMeshPlane::verbose_;
+ }
+
+ //- Prepare limits, scan areas
+ void prepare
+ (
+ std::initializer_list<UList<PDRobstacle>> obstacles,
+ const PDRfitMeshParams& fitParams,
+ scalar cellWidth
+ );
+
+public:
+
+ // Default Constructors
+
+
+ // Member Functions
+
+ //- Calculate the grid controls for the given obstacles
+ //- and parameters
+ Vector<PDRblock::gridControl>
+ calcGriding
+ (
+ std::initializer_list<UList<PDRobstacle>> 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<PDRobstacle>& obstacles);
+
+ //- Populate with areas/positions of the obstacles,
+ //- and check for sub-grid obstacles
+ // Sets the internal
+ void scanAreas
+ (
+ std::initializer_list<UList<PDRobstacle>> 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<PDRobstacle>& obstacles,
+ const scalar minSubgridLen
+ );
+
+ //- Populate with areas/positions of the obstacles,
+ //- without checks for sub-grid obstacles
+ void scanAreas(const UList<PDRobstacle>& 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 4ba931c358b451fb1cfd1f65aca348f73eed2c90..cd1922f2b35736211f4380615bd2b9c1f7393483 100644
--- a/applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C
+++ b/applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C
@@ -83,6 +83,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"
@@ -90,8 +92,6 @@ int main(int argc, char* argv[])
IOdictionary setFieldsDict(dictIO);
- const bool dryrun = args.found("dry-run");
-
const fileName& casepath = runTime.globalPath();
pars.timeName = "0";
diff --git a/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRfitMeshDict b/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRfitMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..fb634d2b92c030b9b683a0513571444a9de8a182
--- /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 3609db612c76d61df431682f2b360ec8079ef3d6..239dd273c4384edbac0a716a85162655a80a420c 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