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

ENH: use MinMax for bounds management in binModels (code simplication)

Browse Source
This commit is contained in:
Mark Olesen
2023-07-06 17:07:08 +02:00
parent dc95242cd2
commit d65e2d89b5
7 changed files with 223 additions and 170 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/functionObjects/field/binField/binModels/binModel/binModel.C
View File

@ -174,7 +174,7 @@ bool Foam::binModel::read(const dictionary& dict)
decomposePatchValues_ = dict.getOrDefault("decomposePatchValues", false);
filePtrs_.setSize(fieldNames_.size());
filePtrs_.resize(fieldNames_.size());
forAll(filePtrs_, i)
{
filePtrs_.set(i, newFileAtStartTime(fieldNames_[i] + "Bin"));

120
src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.C
View File

@ -46,46 +46,69 @@ void Foam::binModels::singleDirectionUniformBin::initialise()
{
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
// Determine extents of patches in a given direction
scalar geomMin = GREAT;
scalar geomMax = -GREAT;
for (const label patchi : patchIDs_)
{
const polyPatch& pp = pbm[patchi];
const scalarField d(pp.faceCentres() & binDir_);
geomMin = min(min(d), geomMin);
geomMax = max(max(d), geomMax);
}
for (const label zonei : cellZoneIDs_)
{
const cellZone& cZone = mesh_.cellZones()[zonei];
const vectorField cz(mesh_.C(), cZone);
const scalarField d(cz & binDir_);
geomMin = min(min(d), geomMin);
geomMax = max(max(d), geomMax);
}
reduce(geomMin, minOp<scalar>());
reduce(geomMax, maxOp<scalar>());
// Slightly boost max so that region of interest is fully within bounds
geomMax = 1.0001*(geomMax - geomMin) + geomMin;
// Use geometry limits if not specified by the user
if (binMin_ == GREAT) binMin_ = geomMin;
if (binMax_ == GREAT) binMax_ = geomMax;
const bool useGeomLimits
(
binLimits_.min() == GREAT
|| binLimits_.max() == GREAT
);
binDx_ = (binMax_ - binMin_)/scalar(nBin_);
if (useGeomLimits)
{
// Determine extents of patches/cells in a given direction
scalarMinMax geomLimits;
if (binDx_ <= 0)
for (const label patchi : patchIDs_)
{
for (const vector& p : pbm[patchi].faceCentres())
{
geomLimits.add(p & binDir_);
}
}
for (const label zonei : cellZoneIDs_)
{
for (const label celli : mesh_.cellZones()[zonei])
{
geomLimits.add(mesh_.C()[celli] & binDir_);
}
}
// Globally consistent
reduce(geomLimits, minMaxOp<scalar>());
if (!geomLimits.good())
{
FatalErrorInFunction
<< "No patches/cellZones provided"
<< exit(FatalError);
}
// Slightly boost max so that region of interest is fully within bounds
// TBD: also adjust min?
const scalar adjust(1e-4*geomLimits.span());
geomLimits.max() += adjust;
// Use geometry limits if not specified by the user
if (binLimits_.min() == GREAT)
{
binLimits_.min() = geomLimits.min();
}
if (binLimits_.max() == GREAT)
{
binLimits_.max() = geomLimits.max();
}
}
binWidth_ = binLimits_.span()/scalar(nBin_);
if (binWidth_ <= 0)
{
FatalErrorInFunction
<< "Max bound must be greater than min bound" << nl
<< " d = " << binDx_ << nl
<< " min = " << binMin_ << nl
<< " max = " << binMax_ << nl
<< " d = " << binWidth_ << nl
<< " min = " << binLimits_.min() << nl
<< " max = " << binLimits_.max() << nl
<< exit(FatalError);
}
}
@ -101,9 +124,8 @@ Foam::binModels::singleDirectionUniformBin::singleDirectionUniformBin
)
:
binModel(dict, mesh, outputPrefix),
binDx_(0),
binMin_(GREAT),
binMax_(GREAT),
binWidth_(0),
binLimits_(GREAT),
binDir_(Zero)
{
read(dict);
@ -125,30 +147,30 @@ bool Foam::binModels::singleDirectionUniformBin::read(const dictionary& dict)
nBin_ = binDict.getCheck<label>("nBin", labelMinMax::ge(1));
Info<< " Employing " << nBin_ << " bins" << endl;
if (binDict.readIfPresent("min", binMin_))
Info<< " Employing " << nBin_ << " bins" << nl;
if (binDict.readIfPresent("min", binLimits_.min()))
{
Info<< " - min : " << binMin_ << endl;
Info<< " - min : " << binLimits_.min() << nl;
}
if (binDict.readIfPresent("max", binMax_))
if (binDict.readIfPresent("max", binLimits_.max()))
{
Info<< " - max : " << binMax_ << endl;
Info<< " - max : " << binLimits_.max() << nl;
}
cumulative_ = binDict.getOrDefault<bool>("cumulative", false);
Info<< " - cumulative : " << cumulative_ << endl;
Info<< " - decomposePatchValues : " << decomposePatchValues_ << endl;
Info<< " - cumulative : " << cumulative_ << nl
<< " - decomposePatchValues : " << decomposePatchValues_ << nl;
binDir_ = binDict.get<vector>("direction");
binDir_.normalise();
if (mag(binDir_) == 0)
if (binDir_.mag() < SMALL)
{
FatalIOErrorInFunction(dict)
<< "Input direction should not be zero valued" << nl
<< " direction = " << binDir_ << nl
<< exit(FatalIOError);
}
binDir_.normalise();
Info<< " - direction : " << binDir_ << nl << endl;
@ -163,11 +185,13 @@ void Foam::binModels::singleDirectionUniformBin::apply()
forAll(fieldNames_, i)
{
const bool ok =
(
processField<scalar>(i)
|| processField<vector>(i)
|| processField<sphericalTensor>(i)
|| processField<symmTensor>(i)
|| processField<tensor>(i);
|| processField<tensor>(i)
);
if (!ok)
{

11
src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBin.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021-2022 OpenCFD Ltd.
Copyright (C) 2021-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -105,13 +105,10 @@ protected:
// Protected Data
//- Distance between bin divisions
scalar binDx_;
scalar binWidth_;
//- Minimum bin bound
scalar binMin_;
//- Maximum bin bound
scalar binMax_;
//- The min/max bounds for the bins
MinMax<scalar> binLimits_;
//- Binning direction
vector binDir_;

67
src/functionObjects/field/binField/binModels/singleDirectionUniformBin/singleDirectionUniformBinTemplates.C
View File

@ -25,7 +25,6 @@ License
\*---------------------------------------------------------------------------*/
template<class Type>
void Foam::binModels::singleDirectionUniformBin::writeFileHeader
(
@ -33,9 +32,9 @@ void Foam::binModels::singleDirectionUniformBin::writeFileHeader
) const
{
writeHeaderValue(os, "bins", nBin_);
writeHeaderValue(os, "start", binMin_);
writeHeaderValue(os, "end", binMax_);
writeHeaderValue(os, "delta", binDx_);
writeHeaderValue(os, "start", binLimits_.min());
writeHeaderValue(os, "end", binLimits_.max());
writeHeaderValue(os, "delta", binWidth_);
writeHeaderValue(os, "direction", binDir_);
// Compute and print bin end points in the binning direction
@ -43,7 +42,7 @@ void Foam::binModels::singleDirectionUniformBin::writeFileHeader
writeCommented(os, "x co-ords :");
forAll(binPoints, pointi)
{
binPoints[pointi] = (binMin_ + (pointi + 1)*binDx_)*binDir_;
binPoints[pointi] = (binLimits_.min() + (pointi + 1)*binWidth_)*binDir_;
os << tab << binPoints[pointi].x();
}
os << nl;
@ -80,7 +79,6 @@ void Foam::binModels::singleDirectionUniformBin::writeFileHeader
{
writeTabbed(os, writeComponents<Type>("patch" + ibin));
}
}
os << endl;
@ -130,27 +128,39 @@ bool Foam::binModels::singleDirectionUniformBin::processField
List<List<Type>> data(nField);
for (auto& binList : data)
{
binList.setSize(nBin_, Zero);
binList.resize(nBin_, Zero);
}
const auto whichBin = [&](const scalar d) -> label
{
if (d >= binLimits_.min() && d <= binLimits_.max())
{
// Find the bin division
label bini = floor
(
(d - binLimits_.min())/binWidth_
);
return min(max(bini, 0), nBin_ - 1);
}
else
{
return -1;
}
};
for (const label zonei : cellZoneIDs_)
{
const cellZone& cZone = mesh_.cellZones()[zonei];
for (const label celli : cZone)
{
const scalar dd = mesh_.C()[celli] & binDir_;
const label bini = whichBin(mesh_.C()[celli] & binDir_);
if (dd < binMin_ || dd > binMax_)
if (bini >= 0)
{
continue;
data[0][bini] += fld[celli];
}
// Find the bin division corresponding to the cell
const label bini =
min(max(floor((dd - binMin_)/binDx_), 0), nBin_ - 1);
data[0][bini] += fld[celli];
}
}
@ -159,25 +169,22 @@ bool Foam::binModels::singleDirectionUniformBin::processField
const polyPatch& pp = mesh_.boundaryMesh()[patchi];
const vectorField np(mesh_.boundary()[patchi].nf());
const auto& pts = pp.faceCentres();
const scalarField dd(pp.faceCentres() & binDir_);
forAll(dd, facei)
forAll(pts, facei)
{
// Avoid faces outside of the bin
if (dd[facei] < binMin_ || dd[facei] > binMax_)
const label bini = whichBin(pts[facei] & binDir_);
if (bini >= 0)
{
continue;
}
const Type& v = fld.boundaryField()[patchi][facei];
// Find the bin division corresponding to the face
const label bini =
min(max(floor((dd[facei] - binMin_)/binDx_), 0), nBin_ - 1);
const Type& v = fld.boundaryField()[patchi][facei];
if (!decomposePatchValues(data, bini, v, np[facei]))
{
data[1][bini] += v;
if (!decomposePatchValues(data, bini, v, np[facei]))
{
data[1][bini] += v;
}
}
}
}

169
src/functionObjects/field/binField/binModels/uniformBin/uniformBin.C
View File

@ -45,70 +45,87 @@ void Foam::binModels::uniformBin::initialise()
{
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
// Determine extents of patches in a given coordinate system
vector geomMin(GREAT, GREAT, GREAT);
vector geomMax(-GREAT, -GREAT, -GREAT);
for (const label patchi : patchIDs_)
// Use geometry limits if not specified by the user
{
const polyPatch& pp = pbm[patchi];
const vectorField ppcs(coordSysPtr_->localPosition(pp.faceCentres()));
// Determine extents of patches/cells
boundBox geomLimits;
for (direction i = 0; i < vector::nComponents; ++i)
for (const label patchi : patchIDs_)
{
geomMin[i] = min(min(ppcs.component(i)), geomMin[i]);
geomMax[i] = max(max(ppcs.component(i)), geomMax[i]);
vectorField pts
(
coordSysPtr_->localPosition(pbm[patchi].faceCentres())
);
MinMax<vector> limits(pts);
geomLimits.add(limits.min());
geomLimits.add(limits.max());
}
}
for (const label zonei : cellZoneIDs_)
{
const cellZone& cZone = mesh_.cellZones()[zonei];
const vectorField d
(
coordSysPtr_->localPosition(vectorField(mesh_.C(), cZone))
);
for (direction i = 0; i < vector::nComponents; ++i)
for (const label zonei : cellZoneIDs_)
{
geomMin[i] = min(min(d.component(i)), geomMin[i]);
geomMax[i] = max(max(d.component(i)), geomMax[i]);
const cellZone& cZone = mesh_.cellZones()[zonei];
const vectorField pts
(
coordSysPtr_->localPosition(vectorField(mesh_.C(), cZone))
);
MinMax<vector> limits(pts);
geomLimits.add(limits.min());
geomLimits.add(limits.max());
}
}
reduce(geomMin, minOp<vector>());
reduce(geomMax, maxOp<vector>());
// Globally consistent
geomLimits.reduce();
for (direction i = 0; i < vector::nComponents; ++i)
{
// Slightly boost max so that region of interest is fully within bounds
geomMax[i] = 1.0001*(geomMax[i] - geomMin[i]) + geomMin[i];
// TBD: could also adjust min?
const vector adjust(1e-4*geomLimits.span());
geomLimits.max() += adjust;
// Use geometry limits if not specified by the user
if (binMinMax_[i][0] == GREAT) binMinMax_[i][0] = geomMin[i];
if (binMinMax_[i][1] == GREAT) binMinMax_[i][1] = geomMax[i];
for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
{
// Use geometry limits if not specified by the user
if (binLimits_.min()[cmpt] == GREAT)
{
binLimits_.min()[cmpt] = geomLimits.min()[cmpt];
}
if (binLimits_.max()[cmpt] == GREAT)
{
binLimits_.max()[cmpt] = geomLimits.max()[cmpt];
}
}
}
if (binMinMax_[i][0] > binMinMax_[i][1])
for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
{
if (binLimits_.min()[cmpt] > binLimits_.max()[cmpt])
{
FatalErrorInFunction
<< "Max bounds must be greater than min bounds" << nl
<< " direction = " << i << nl
<< " min = " << binMinMax_[i][0] << nl
<< " max = " << binMinMax_[i][1] << nl
<< " direction = " << cmpt << nl
<< " min = " << binLimits_.min()[cmpt] << nl
<< " max = " << binLimits_.max()[cmpt] << nl
<< exit(FatalError);
}
//- Compute bin widths in binning directions
binW_[i] = (binMinMax_[i][1] - binMinMax_[i][0])/scalar(nBins_[i]);
binWidth_[cmpt] =
(
(binLimits_.max()[cmpt] - binLimits_.min()[cmpt])
/ scalar(nBins_[cmpt])
);
if (binW_[i] <= 0)
if (binWidth_[cmpt] <= 0)
{
FatalErrorInFunction
<< "Bin widths must be greater than zero" << nl
<< " direction = " << i << nl
<< " min bound = " << binMinMax_[i][0] << nl
<< " max bound = " << binMinMax_[i][1] << nl
<< " bin width = " << binW_[i]
<< " direction = " << cmpt << nl
<< " min bound = " << binLimits_.min()[cmpt] << nl
<< " max bound = " << binLimits_.max()[cmpt] << nl
<< " bin width = " << binWidth_[cmpt] << nl
<< exit(FatalError);
}
}
@ -125,9 +142,13 @@ Foam::labelList Foam::binModels::uniformBin::binAddr(const vectorField& d) const
{
// Avoid elements outside of the bin
bool faceInside = true;
for (direction j = 0; j < vector::nComponents; ++j)
for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
{
if (d[i][j] < binMinMax_[j][0] || d[i][j] > binMinMax_[j][1])
if
(
d[i][cmpt] < binLimits_.min()[cmpt]
|| d[i][cmpt] > binLimits_.max()[cmpt]
)
{
faceInside = false;
break;
@ -138,14 +159,18 @@ Foam::labelList Foam::binModels::uniformBin::binAddr(const vectorField& d) const
{
// Find the bin division corresponding to the element
Vector<label> n(Zero);
for (direction j = 0; j < vector::nComponents; ++j)
for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
{
n[j] = floor((d[i][j] - binMinMax_[j][0])/binW_[j]);
n[j] = min(max(n[j], 0), nBins_[j] - 1);
label bini = floor
(
(d[i][cmpt] - binLimits_.min()[cmpt])/binWidth_[cmpt]
);
n[cmpt] = min(max(bini, 0), nBins_[cmpt] - 1);
}
// Order: (e1, e2, e3), the first varies the fastest
binIndices[i] = n[0] + nBins_[0]*n[1] + nBins_[0]*nBins_[1]*n[2];
binIndices[i] = n.x() + nBins_[0]*n.y() + nBins_[0]*nBins_[1]*n.z();
}
else
{
@ -202,13 +227,8 @@ Foam::binModels::uniformBin::uniformBin
:
binModel(dict, mesh, outputPrefix),
nBins_(Zero),
binW_(Zero),
binMinMax_
(
vector2D(GREAT, GREAT),
vector2D(GREAT, GREAT),
vector2D(GREAT, GREAT)
)
binWidth_(Zero),
binLimits_(vector::uniform(GREAT))
{
read(dict);
}
@ -238,36 +258,40 @@ bool Foam::binModels::uniformBin::read(const dictionary& dict)
{
FatalIOErrorInFunction(binDict)
<< "Number of bins must be greater than zero" << nl
<< " e1 bins = " << nBins_[0] << nl
<< " e2 bins = " << nBins_[1] << nl
<< " e3 bins = " << nBins_[2]
<< " e1 bins = " << nBins_.x() << nl
<< " e2 bins = " << nBins_.y() << nl
<< " e3 bins = " << nBins_.z()
<< exit(FatalIOError);
}
Info<< " - Employing:" << nl
<< " " << nBins_[0] << " e1 bins," << nl
<< " " << nBins_[1] << " e2 bins," << nl
<< " " << nBins_[2] << " e3 bins"
<< " " << nBins_.x() << " e1 bins," << nl
<< " " << nBins_.y() << " e2 bins," << nl
<< " " << nBins_.z() << " e3 bins"
<< endl;
cumulative_ = binDict.getOrDefault<bool>("cumulative", false);
Info<< " - cumulative : " << cumulative_ << endl;
Info<< " - decomposePatchValues : " << decomposePatchValues_ << endl;
if (binDict.found("minMax"))
const dictionary* minMaxDictPtr = binDict.findDict("minMax");
if (minMaxDictPtr)
{
const dictionary& minMaxDict = binDict.subDict("minMax");
const auto& minMaxDict = *minMaxDictPtr;
for (direction i = 0; i < vector::nComponents; ++i)
for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
{
const word ei("e" + Foam::name(i));
const word ei("e" + Foam::name(cmpt));
if (minMaxDict.readIfPresent(ei, binMinMax_[i]))
scalarMinMax range;
if (minMaxDict.readIfPresent(ei, range))
{
Info<< " - " << ei << " min : "
<< binMinMax_[i][0] << nl
<< " - " << ei << " max : "
<< binMinMax_[i][1] << endl;
binLimits_.min()[cmpt] = range.min();
binLimits_.max()[cmpt] = range.max();
Info<< " - " << ei << " min/max : " << range << nl;
}
}
}
@ -284,17 +308,18 @@ void Foam::binModels::uniformBin::apply()
forAll(fieldNames_, i)
{
const bool ok =
(
processField<scalar>(i)
|| processField<vector>(i)
|| processField<sphericalTensor>(i)
|| processField<symmTensor>(i)
|| processField<tensor>(i);
|| processField<tensor>(i)
);
if (!ok)
{
WarningInFunction
<< "Unable to find field " << fieldNames_[i]
<< endl;
<< "Unable to find field " << fieldNames_[i] << endl;
}
}

8
src/functionObjects/field/binField/binModels/uniformBin/uniformBin.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021-2022 OpenCFD Ltd.
Copyright (C) 2021-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -129,10 +129,10 @@ protected:
Vector<label> nBins_;
//- Equidistant bin widths in binning directions
vector binW_;
vector binWidth_;
//- Min-max bounds of bins in binning directions
Vector<vector2D> binMinMax_;
//- The geometric min/max bounds for the bins
MinMax<vector> binLimits_;
//- Face index to bin index addressing
labelList faceToBin_;

16
src/functionObjects/field/binField/binModels/uniformBin/uniformBinTemplates.C
View File

@ -39,9 +39,9 @@ void Foam::binModels::uniformBin::writeFileHeader
for (direction i = 0; i < vector::nComponents; ++i)
{
writeHeaderValue(os, "e" + Foam::name(i) + " bins", nBins_[i]);
writeHeaderValue(os, " start", binMinMax_[i][0]);
writeHeaderValue(os, " end", binMinMax_[i][1]);
writeHeaderValue(os, " delta", binW_[i]);
writeHeaderValue(os, " start", binLimits_.min()[i]);
writeHeaderValue(os, " end", binLimits_.max()[i]);
writeHeaderValue(os, " delta", binWidth_[i]);
writeHeaderValue(os, " direction", R.col(i));
}
writeCommented(os, "bin end co-ordinates:");
@ -50,12 +50,12 @@ void Foam::binModels::uniformBin::writeFileHeader
// Compute and print bin end points in binning directions
for (direction i = 0; i < vector::nComponents; ++i)
{
scalar binEnd = binMinMax_[i][0];
scalar binEnd = binLimits_.min()[i];
writeCommented(os, "e"+Foam::name(i)+" co-ords :");
for (label j = 0; j < nBins_[i]; ++j)
{
binEnd += binW_[i];
binEnd += binWidth_[i];
os << tab << binEnd;
}
os << nl;
@ -124,7 +124,7 @@ bool Foam::binModels::uniformBin::processField(const label fieldi)
List<List<Type>> data(nField);
for (auto& binList : data)
{
binList.setSize(nBin_, Zero);
binList.resize(nBin_, Zero);
}
for (const label zonei : cellZoneIDs_)
@ -135,7 +135,7 @@ bool Foam::binModels::uniformBin::processField(const label fieldi)
{
const label bini = cellToBin_[celli];
if (bini != -1)
if (bini >= 0)
{
data[0][bini] += fld[celli];
}
@ -153,7 +153,7 @@ bool Foam::binModels::uniformBin::processField(const label fieldi)
pp.start() - mesh_.nInternalFaces() + facei;
const label bini = faceToBin_[localFacei];
if (bini != -1)
if (bini >= 0)
{
const Type& v = fld.boundaryField()[patchi][facei];