zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-PDRfitMesh' into 'develop'

Browse Source 
Draft: Integration of PDRfitMesh

See merge request Development/openfoam!415
This commit is contained in:
Mark OLESEN
2025-10-23 22:38:43 +00:00
parent 6b3e65d17d 1ed3fcd4c0
commit 0336bfe699
15 changed files with 2601 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
applications/utilities/preProcessing/PDR/Allwclean
View File

@ -1,7 +1,11 @@
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
#------------------------------------------------------------------------------
wclean libso pdrFields
wclean PDRfitMesh
wclean PDRsetFields
#------------------------------------------------------------------------------

3
applications/utilities/preProcessing/PDR/Allwmake
View File

@ -5,6 +5,9 @@ cd "${0%/*}" || exit # Run from this directory
#------------------------------------------------------------------------------
wmake $targetType pdrFields
wmake $targetType PDRfitMesh
wmake $targetType PDRsetFields
#------------------------------------------------------------------------------

7
applications/utilities/preProcessing/PDR/PDRfitMesh/Make/files Normal file
View File

@ -0,0 +1,7 @@
PDRfitMeshParams.C
PDRfitMeshScan.C
PDRfitMeshScans.C
PDRfitMesh.C
EXE = $(FOAM_APPBIN)/PDRfitMesh

15
applications/utilities/preProcessing/PDR/PDRfitMesh/Make/options Normal file
View File

@ -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

339
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMesh.C Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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<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, // 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<PDRblock::gridControl> 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<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",
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;
}
// ************************************************************************* //

126
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshDict Normal file
View File

@ -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.
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;
// ************************************************************************* //

123
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.C Normal file
View File

@ -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("nEdgeLayers", nEdgeLayers);
dict.readIfPresent("outerRatio", outerRatio);
// dict.readIfPresent("outerRadius", outerRadius);
};
// ************************************************************************* //

158
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshParams.H Normal file
View File

@ -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 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
// ************************************************************************* //

859
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.C Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<scalar> 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<scalar>(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;
}
// ************************************************************************* //

224
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScan.H Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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<scalar> weightedPos_;
//- Area-weights for position
DynamicList<scalar> totalArea_;
//- Area-weighted half-cell position
DynamicList<scalar> weightedPos2_;
//- Area-weights for half-cell position
DynamicList<scalar> 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
// ************************************************************************* //

464
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.C Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "PDRfitMeshScans.H"
#include "PDRobstacle.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::PDRfitMeshScans::prepare
(
const UList<PDRobstacle>& 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<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::PDRfitMeshScans::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::PDRfitMeshScans::scanAreas
(
const UList<PDRobstacle>& obstacles
)
{
(void)scanAreas(obstacles, -GREAT);
}
void Foam::PDRfitMeshScans::scanAreas
(
const UList<PDRobstacle>& 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::PDRblock::gridControl>
Foam::PDRfitMeshScans::calcGriding
(
const UList<PDRobstacle>& 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<PDRblock::gridControl> griding;
for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
{
griding[cmpt] =
(*this)[cmpt].calcGridControl(cellWidth, max_zone, fitParams);
}
return griding;
}
// ************************************************************************* //

154
applications/utilities/preProcessing/PDR/PDRfitMesh/PDRfitMeshScans.H Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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<PDRfitMeshScan>
{
// Private Data
//- The evaluated sub-grid length
scalar subgridLen_ = 0;
// Private Member Functions
//- Prepare limits, scan areas
void prepare
(
const UList<PDRobstacle>& 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<PDRblock::gridControl>
calcGriding
(
const 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
(
const 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
// ************************************************************************* //

2
applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C
View File

@ -82,6 +82,8 @@ int main(int argc, char* argv[])
#include "setRootCase.H"
#include "createTime.H"
const bool dryrun = args.found("dry-run");
const word dictName("PDRsetFieldsDict");
#include "setSystemRunTimeDictionaryIO.H"

105
tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRfitMeshDict Normal file
View File

@ -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;
// ************************************************************************* //

18
tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRsetFieldsDict
View File

@ -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