From 1ed3fcd4c0b79ba829d279d7582409c8c5fb0839 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Mon, 21 Dec 2020 09:02:29 +0100 Subject: [PATCH] ENH: integration of PDRfitMesh (#1961) - serves as a pre-processor for PDRblockMesh. Scans obstacles to determine an initial cell distribution for PDRblockMesh --- .../utilities/preProcessing/PDR/Allwclean | 4 + .../utilities/preProcessing/PDR/Allwmake | 3 + .../preProcessing/PDR/PDRfitMesh/Make/files | 7 + .../preProcessing/PDR/PDRfitMesh/Make/options | 15 + .../preProcessing/PDR/PDRfitMesh/PDRfitMesh.C | 339 +++++++ .../PDR/PDRfitMesh/PDRfitMeshDict | 126 +++ .../PDR/PDRfitMesh/PDRfitMeshParams.C | 123 +++ .../PDR/PDRfitMesh/PDRfitMeshParams.H | 158 ++++ .../PDR/PDRfitMesh/PDRfitMeshScan.C | 859 ++++++++++++++++++ .../PDR/PDRfitMesh/PDRfitMeshScan.H | 224 +++++ .../PDR/PDRfitMesh/PDRfitMeshScans.C | 464 ++++++++++ .../PDR/PDRfitMesh/PDRfitMeshScans.H | 154 ++++ .../PDR/PDRsetFields/PDRsetFields.C | 4 +- .../simplePipeCage/system/PDRfitMeshDict | 105 +++ .../simplePipeCage/system/PDRsetFieldsDict | 18 + 15 files changed, 2601 insertions(+), 2 deletions(-) create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/Make/files create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/Make/options create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMesh.C create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshDict create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.C create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.H create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.C create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.H create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.C create mode 100644 applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.H create mode 100644 tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRfitMeshDict 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 4ba931c358..cd1922f2b3 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 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 3609db612c..239dd273c4 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