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

GIT: Merge/resolve conflict

Browse Source
This commit is contained in:
andy
2012-08-03 12:28:20 +01:00
parent d9912c2003 fc1bc5c5e2
commit a3870d65df
316 changed files with 5877 additions and 3138 deletions
Download Patch File Download Diff File Expand all files Collapse all files

788
src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C
View File

@ -27,8 +27,8 @@ License
#include "volFields.H"
#include "regionSplit.H"
#include "fvcVolumeIntegrate.H"
#include "Histogram.H"
#include "mathematicalConstants.H"
#include "stringListOps.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -56,6 +56,281 @@ namespace Foam
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::regionSizeDistribution::writeGraph
(
const coordSet& coords,
const word& valueName,
const scalarField& values
) const
{
const fvMesh& mesh = refCast<const fvMesh>(obr_);
const wordList valNames(1, valueName);
fileName outputPath;
if (Pstream::parRun())
{
outputPath = mesh.time().path()/".."/name_;
}
else
{
outputPath = mesh.time().path()/name_;
}
if (mesh.name() != fvMesh::defaultRegion)
{
outputPath = outputPath/mesh.name();
}
mkDir(outputPath/mesh.time().timeName());
OFstream str
(
outputPath
/ mesh.time().timeName()
/ formatterPtr_().getFileName(coords, valNames)
);
Info<< "Writing distribution of " << valueName << " to " << str.name()
<< endl;
List<const scalarField*> valPtrs(1);
valPtrs[0] = &values;
formatterPtr_().write(coords, valNames, valPtrs, str);
}
void Foam::regionSizeDistribution::writeAlphaFields
(
const regionSplit& regions,
const Map<label>& patchRegions,
const Map<scalar>& regionVolume,
const volScalarField& alpha
) const
{
const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
// Split alpha field
// ~~~~~~~~~~~~~~~~~
// Split into
// - liquidCore : region connected to inlet patches
// - per region a volume : for all other regions
// - backgroundAlpha : remaining alpha
// Construct field
volScalarField liquidCore
(
IOobject
(
alphaName_ + "_liquidCore",
obr_.time().timeName(),
obr_,
IOobject::NO_READ
),
alpha,
fvPatchField<scalar>::calculatedType()
);
volScalarField backgroundAlpha
(
IOobject
(
alphaName_ + "_background",
obr_.time().timeName(),
obr_,
IOobject::NO_READ
),
alpha,
fvPatchField<scalar>::calculatedType()
);
// Knock out any cell not in patchRegions
forAll(liquidCore, cellI)
{
label regionI = regions[cellI];
if (patchRegions.found(regionI))
{
backgroundAlpha[cellI] = 0;
}
else
{
liquidCore[cellI] = 0;
scalar regionVol = regionVolume[regionI];
if (regionVol < maxDropletVol)
{
backgroundAlpha[cellI] = 0;
}
}
}
liquidCore.correctBoundaryConditions();
backgroundAlpha.correctBoundaryConditions();
Info<< "Volume of liquid-core = "
<< fvc::domainIntegrate(liquidCore).value()
<< endl;
Info<< "Volume of background = "
<< fvc::domainIntegrate(backgroundAlpha).value()
<< endl;
Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
liquidCore.write();
Info<< "Writing background field to " << backgroundAlpha.name() << endl;
backgroundAlpha.write();
}
Foam::Map<Foam::label> Foam::regionSizeDistribution::findPatchRegions
(
const polyMesh& mesh,
const regionSplit& regions
) const
{
// Mark all regions starting at patches
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Count number of patch faces (just for initial sizing)
const labelHashSet patchIDs(mesh.boundaryMesh().patchSet(patchNames_));
label nPatchFaces = 0;
forAllConstIter(labelHashSet, patchIDs, iter)
{
nPatchFaces += mesh.boundaryMesh()[iter.key()].size();
}
Map<label> patchRegions(nPatchFaces);
forAllConstIter(labelHashSet, patchIDs, iter)
{
const polyPatch& pp = mesh.boundaryMesh()[iter.key()];
// Collect all regions on the patch
const labelList& faceCells = pp.faceCells();
forAll(faceCells, i)
{
patchRegions.insert
(
regions[faceCells[i]],
Pstream::myProcNo() // dummy value
);
}
}
// Make sure all the processors have the same set of regions
Pstream::mapCombineGather(patchRegions, minEqOp<label>());
Pstream::mapCombineScatter(patchRegions);
return patchRegions;
}
Foam::tmp<Foam::scalarField> Foam::regionSizeDistribution::divide
(
const scalarField& num,
const scalarField& denom
)
{
tmp<scalarField> tresult(new scalarField(num.size()));
scalarField& result = tresult();
forAll(denom, i)
{
if (denom[i] != 0)
{
result[i] = num[i]/denom[i];
}
else
{
result[i] = 0.0;
}
}
return tresult;
}
void Foam::regionSizeDistribution::writeGraphs
(
const word& fieldName, // name of field
const labelList& indices, // index of bin for each region
const scalarField& sortedField, // per region field data
const scalarField& binCount, // per bin number of regions
const coordSet& coords // graph data for bins
) const
{
if (Pstream::master())
{
// Calculate per-bin average
scalarField binSum(nBins_, 0.0);
forAll(sortedField, i)
{
binSum[indices[i]] += sortedField[i];
}
scalarField binAvg(divide(binSum, binCount));
// Per bin deviation
scalarField binSqrSum(nBins_, 0.0);
forAll(sortedField, i)
{
binSqrSum[indices[i]] += Foam::sqr(sortedField[i]);
}
scalarField binDev
(
sqrt(divide(binSqrSum, binCount) - Foam::sqr(binAvg))
);
// Write average
writeGraph(coords, fieldName + "_sum", binSum);
// Write average
writeGraph(coords, fieldName + "_avg", binAvg);
// Write deviation
writeGraph(coords, fieldName + "_dev", binDev);
}
}
void Foam::regionSizeDistribution::writeGraphs
(
const word& fieldName, // name of field
const scalarField& cellField, // per cell field data
const regionSplit& regions, // per cell the region(=droplet)
const labelList& sortedRegions, // valid regions in sorted order
const scalarField& sortedNormalisation,
const labelList& indices, // per region index of bin
const scalarField& binCount, // per bin number of regions
const coordSet& coords // graph data for bins
) const
{
// Sum on a per-region basis. Parallel reduced.
Map<scalar> regionField(regionSum(regions, cellField));
// Extract in region order
scalarField sortedField
(
sortedNormalisation
* extractData
(
sortedRegions,
regionField
)
);
writeGraphs
(
fieldName, // name of field
indices, // index of bin for each region
sortedField, // per region field data
binCount, // per bin number of regions
coords // graph data for bins
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::regionSizeDistribution::regionSizeDistribution
@ -103,8 +378,11 @@ void Foam::regionSizeDistribution::read(const dictionary& dict)
dict.lookup("field") >> alphaName_;
dict.lookup("patches") >> patchNames_;
dict.lookup("threshold") >> threshold_;
dict.lookup("volFraction") >> volFraction_;
dict.lookup("maxDiameter") >> maxDiam_;
minDiam_ = 0.0;
dict.readIfPresent("minDiameter", minDiam_);
dict.lookup("nBins") >> nBins_;
dict.lookup("fields") >> fields_;
word format(dict.lookup("setFormat"));
formatterPtr_ = writer<scalar>::New(format);
@ -163,14 +441,17 @@ void Foam::regionSizeDistribution::write()
: obr_.lookupObject<volScalarField>(alphaName_)
);
Info<< "Volume of alpha = "
Info<< "Volume of alpha = "
<< fvc::domainIntegrate(alpha).value()
<< endl;
const scalar meshVol = gSum(mesh.V());
Info<< "Mesh volume = " << meshVol << endl;
Info<< "Background region volume limit = " << volFraction_*meshVol
<< endl;
const scalar maxDropletVol = 1.0/6.0*pow(maxDiam_, 3);
const scalar delta = (maxDiam_-minDiam_)/nBins_;
Info<< "Mesh volume = " << meshVol << endl;
Info<< "Maximum droplet diameter = " << maxDiam_ << endl;
Info<< "Maximum droplet volume = " << maxDropletVol << endl;
// Determine blocked faces
@ -260,138 +541,115 @@ void Foam::regionSizeDistribution::write()
}
// Determine regions connected to supplied patches
Map<label> patchRegions(findPatchRegions(mesh, regions));
// Sum all regions
Map<Pair<scalar> > regionVolume(regions.nRegions()/Pstream::nProcs());
forAll(alpha, cellI)
{
scalar cellVol = mesh.V()[cellI];
scalar alphaVol = alpha[cellI]*cellVol;
label regionI = regions[cellI];
Map<Pair<scalar> >::iterator fnd = regionVolume.find(regionI);
if (fnd == regionVolume.end())
{
regionVolume.insert
(
regionI,
Pair<scalar>(cellVol, alphaVol)
);
}
else
{
fnd().first() += cellVol;
fnd().second() += alphaVol;
}
}
Pstream::mapCombineGather(regionVolume, ListPlusEqOp<scalar, 2>());
Pstream::mapCombineScatter(regionVolume);
const scalarField alphaVol(alpha.internalField()*mesh.V());
Map<scalar> allRegionVolume(regionSum(regions, mesh.V()));
Map<scalar> allRegionAlphaVolume(regionSum(regions, alphaVol));
Map<label> allRegionNumCells
(
regionSum
(
regions,
labelField(mesh.nCells(), 1.0)
)
);
if (debug)
{
Info<< token::TAB << "Region"
<< token::TAB << "Volume(mesh)"
<< token::TAB << "Volume(" << alpha.name() << "):"
<< token::TAB << "nCells"
<< endl;
scalar meshSumVol = 0.0;
scalar alphaSumVol = 0.0;
label nCells = 0;
forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
Map<scalar>::const_iterator vIter = allRegionVolume.begin();
Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
Map<label>::const_iterator numIter = allRegionNumCells.begin();
for
(
;
vIter != allRegionVolume.end()
&& aIter != allRegionAlphaVolume.end();
++vIter, ++aIter, ++numIter
)
{
Info<< token::TAB << iter.key()
<< token::TAB << iter().first()
<< token::TAB << iter().second() << endl;
Info<< token::TAB << vIter.key()
<< token::TAB << vIter()
<< token::TAB << aIter()
<< token::TAB << numIter()
<< endl;
meshSumVol += iter().first();
alphaSumVol += iter().second();
meshSumVol += vIter();
alphaSumVol += aIter();
nCells += numIter();
}
Info<< token::TAB << "Total:"
<< token::TAB << meshSumVol
<< token::TAB << alphaSumVol << endl;
<< token::TAB << alphaSumVol
<< token::TAB << nCells
<< endl;
Info<< endl;
}
// Mark all regions starting at patches
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Count number of patch faces (just for initial sizing)
label nPatchFaces = 0;
forAll(patchNames_, i)
{
const word& pName = patchNames_[i];
label patchI = mesh.boundaryMesh().findPatchID(pName);
if (patchI == -1)
Info<< "Patch connected regions (liquid core):" << endl;
Info<< token::TAB << "Region"
<< token::TAB << "Volume(mesh)"
<< token::TAB << "Volume(" << alpha.name() << "):"
<< endl;
forAllConstIter(Map<label>, patchRegions, iter)
{
WarningIn("regionSizeDistribution::write()")
<< "Cannot find patch " << pName << ". Valid patches are "
<< mesh.boundaryMesh().names()
<< endl;
}
else
{
nPatchFaces += mesh.boundaryMesh()[patchI].size();
label regionI = iter.key();
Info<< token::TAB << iter.key()
<< token::TAB << allRegionVolume[regionI]
<< token::TAB << allRegionAlphaVolume[regionI] << endl;
}
Info<< endl;
}
Map<label> keepRegions(nPatchFaces);
forAll(patchNames_, i)
{
const word& pName = patchNames_[i];
Info<< "Background regions:" << endl;
Info<< token::TAB << "Region"
<< token::TAB << "Volume(mesh)"
<< token::TAB << "Volume(" << alpha.name() << "):"
<< endl;
Map<scalar>::const_iterator vIter = allRegionVolume.begin();
Map<scalar>::const_iterator aIter = allRegionAlphaVolume.begin();
label patchI = mesh.boundaryMesh().findPatchID(pName);
if (patchI != -1)
{
const polyPatch& pp = mesh.boundaryMesh()[patchI];
// Collect all regions on the patch
const labelList& faceCells = pp.faceCells();
forAll(faceCells, i)
{
keepRegions.insert
(
regions[faceCells[i]],
Pstream::myProcNo()
);
}
}
}
// Make sure all the processors have the same set of regions
Pstream::mapCombineGather(keepRegions, minEqOp<label>());
Pstream::mapCombineScatter(keepRegions);
Info<< "Patch connected regions (liquid core):" << endl;
forAllConstIter(Map<label>, keepRegions, iter)
{
label regionI = iter.key();
Pair<scalar>& vols = regionVolume[regionI];
Info<< token::TAB << iter.key()
<< token::TAB << vols.first()
<< token::TAB << vols.second() << endl;
}
Info<< endl;
Info<< "Background regions:" << endl;
forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
{
if
for
(
!keepRegions.found(iter.key())
&& iter().first() >= volFraction_*meshVol
;
vIter != allRegionVolume.end()
&& aIter != allRegionAlphaVolume.end();
++vIter, ++aIter
)
{
Info<< token::TAB << iter.key()
<< token::TAB << iter().first()
<< token::TAB << iter().second() << endl;
if
(
!patchRegions.found(vIter.key())
&& vIter() >= maxDropletVol
)
{
Info<< token::TAB << vIter.key()
<< token::TAB << vIter()
<< token::TAB << aIter() << endl;
}
}
Info<< endl;
}
Info<< endl;
// Split alpha field
@ -400,185 +658,197 @@ void Foam::regionSizeDistribution::write()
// - liquidCore : region connected to inlet patches
// - per region a volume : for all other regions
// - backgroundAlpha : remaining alpha
writeAlphaFields(regions, patchRegions, allRegionVolume, alpha);
// Construct field
volScalarField liquidCore
(
IOobject
(
alphaName_ + "_liquidCore",
obr_.time().timeName(),
obr_,
IOobject::NO_READ
),
alpha,
fvPatchField<scalar>::calculatedType()
);
volScalarField backgroundAlpha
(
IOobject
(
alphaName_ + "_background",
obr_.time().timeName(),
obr_,
IOobject::NO_READ
),
alpha,
fvPatchField<scalar>::calculatedType()
);
// Knock out any cell not in keepRegions
forAll(liquidCore, cellI)
// Extract droplet-only allRegionVolume, i.e. delete liquid core
// (patchRegions) and background regions from maps.
// Note that we have to use mesh volume (allRegionVolume) and not
// allRegionAlphaVolume since background might not have alpha in it.
forAllIter(Map<scalar>, allRegionVolume, vIter)
{
label regionI = regions[cellI];
if (keepRegions.found(regionI))
label regionI = vIter.key();
if
(
patchRegions.found(regionI)
|| vIter() >= maxDropletVol
)
{
backgroundAlpha[cellI] = 0;
}
else
{
liquidCore[cellI] = 0;
scalar regionVol = regionVolume[regionI].first();
if (regionVol < volFraction_*meshVol)
{
backgroundAlpha[cellI] = 0;
}
allRegionVolume.erase(vIter);
allRegionAlphaVolume.erase(regionI);
allRegionNumCells.erase(regionI);
}
}
liquidCore.correctBoundaryConditions();
backgroundAlpha.correctBoundaryConditions();
Info<< "Volume of liquid-core = "
<< fvc::domainIntegrate(liquidCore).value()
<< endl;
Info<< "Writing liquid-core field to " << liquidCore.name() << endl;
liquidCore.write();
Info<< "Volume of background = "
<< fvc::domainIntegrate(backgroundAlpha).value()
<< endl;
Info<< "Writing background field to " << backgroundAlpha.name() << endl;
backgroundAlpha.write();
// Collect histogram
if (Pstream::master())
if (allRegionVolume.size())
{
DynamicList<scalar> diameters(regionVolume.size());
forAllConstIter(Map<Pair<scalar> >, regionVolume, iter)
// Construct mids of bins for plotting
pointField xBin(nBins_);
scalar x = 0.5*delta;
forAll(xBin, i)
{
if (!keepRegions.found(iter.key()))
{
if (iter().first() < volFraction_*meshVol)
{
scalar v = iter().second();
//scalar diam = Foam::cbrt(v*6/mathematicalConstant::pi);
scalar diam =
Foam::cbrt(v*6/constant::mathematical::pi);
diameters.append(diam);
}
}
xBin[i] = point(x, 0, 0);
x += delta;
}
if (diameters.size())
const coordSet coords("diameter", "x", xBin, mag(xBin));
// Get in region order the alpha*volume and diameter
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const labelList sortedRegions = allRegionAlphaVolume.sortedToc();
scalarField sortedVols
(
extractData
(
sortedRegions,
allRegionAlphaVolume
)
);
// Calculate the diameters
scalarField sortedDiameters(sortedVols.size());
forAll(sortedDiameters, i)
{
scalar maxDiam = max(diameters);
scalar minDiam = 0.0;
Info<< "Maximum diameter:" << maxDiam << endl;
Histogram<List<scalar> > bins
sortedDiameters[i] = Foam::cbrt
(
minDiam,
maxDiam,
nBins_,
diameters
sortedVols[i]
*6/constant::mathematical::pi
);
}
/* 1.7.x
scalarField xBin(nBins_);
// Determine the bin index for all the diameters
labelList indices(sortedDiameters.size());
forAll(sortedDiameters, i)
{
indices[i] = (sortedDiameters[i]-minDiam_)/delta;
}
scalar dx = (maxDiam-minDiam)/nBins_;
scalar x = 0.5*dx;
forAll(bins.counts(), i)
// Calculate the counts per diameter bin
scalarField binCount(nBins_, 0.0);
forAll(sortedDiameters, i)
{
binCount[indices[i]] += 1.0;
}
// Write counts
if (Pstream::master())
{
writeGraph(coords, "count", binCount);
}
// Write to screen
{
Info<< "Bins:" << endl;
Info<< token::TAB << "Bin"
<< token::TAB << "Min diameter"
<< token::TAB << "Count:"
<< endl;
scalar diam = 0.0;
forAll(binCount, binI)
{
xBin[i] = x;
x += dx;
Info<< token::TAB << binI
<< token::TAB << diam
<< token::TAB << binCount[binI] << endl;
diam += delta;
}
Info<< endl;
}
scalarField normalisedCount(bins.counts().size());
forAll(bins.counts(), i)
// Write average and deviation of droplet volume.
writeGraphs
(
"volume", // name of field
indices, // per region the bin index
sortedVols, // per region field data
binCount, // per bin number of regions
coords // graph data for bins
);
// Collect some more field
{
wordList scalarNames(obr_.names(volScalarField::typeName));
labelList selected = findStrings(fields_, scalarNames);
forAll(selected, i)
{
normalisedCount[i] = 1.0*bins.counts()[i];
const word& fldName = scalarNames[selected[i]];
Info<< "Scalar field " << fldName << endl;
const scalarField& fld = obr_.lookupObject
<
volScalarField
>(fldName).internalField();
writeGraphs
(
fldName, // name of field
alphaVol*fld, // per cell field data
regions, // per cell the region(=droplet)
sortedRegions, // valid regions in sorted order
1.0/sortedVols, // per region normalisation
indices, // index of bin for each region
binCount, // per bin number of regions
coords // graph data for bins
);
}
}
{
wordList vectorNames(obr_.names(volVectorField::typeName));
labelList selected = findStrings(fields_, vectorNames);
const coordSet coords
(
"diameter",
"x",
xBin
);
*/
pointField xBin(nBins_);
scalar dx = (maxDiam - minDiam)/nBins_;
scalar x = 0.5*dx;
forAll(bins.counts(), i)
forAll(selected, i)
{
xBin[i] = point(x, 0, 0);
x += dx;
const word& fldName = vectorNames[selected[i]];
Info<< "Vector field " << fldName << endl;
const vectorField& fld = obr_.lookupObject
<
volVectorField
>(fldName).internalField();
// Components
for (direction cmp = 0; cmp < vector::nComponents; cmp++)
{
writeGraphs
(
fldName + vector::componentNames[cmp],
alphaVol*fld.component(cmp),// per cell field data
regions, // per cell the region(=droplet)
sortedRegions, // valid regions in sorted order
1.0/sortedVols, // per region normalisation
indices, // index of bin for each region
binCount, // per bin number of regions
coords // graph data for bins
);
}
// Magnitude
writeGraphs
(
fldName + "mag", // name of field
alphaVol*mag(fld), // per cell field data
regions, // per cell the region(=droplet)
sortedRegions, // valid regions in sorted order
1.0/sortedVols, // per region normalisation
indices, // index of bin for each region
binCount, // per bin number of regions
coords // graph data for bins
);
}
scalarField normalisedCount(bins.counts().size());
forAll(bins.counts(), i)
{
normalisedCount[i] = 1.0*bins.counts()[i];
}
const coordSet coords
(
"diameter",
"x",
xBin,
mag(xBin)
);
const wordList valNames(1, "count");
fileName outputPath;
if (Pstream::parRun())
{
outputPath = mesh.time().path()/".."/name_;
}
else
{
outputPath = mesh.time().path()/name_;
}
if (mesh.name() != fvMesh::defaultRegion)
{
outputPath = outputPath/mesh.name();
}
mkDir(outputPath/mesh.time().timeName());
OFstream str
(
outputPath
/ mesh.time().timeName()
/ formatterPtr_().getFileName(coords, valNames)
);
Info<< "Writing distribution to " << str.name() << endl;
List<const scalarField*> valPtrs(1);
valPtrs[0] = &normalisedCount;
formatterPtr_().write(coords, valNames, valPtrs, str);
}
}
}

145
src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H
View File

@ -28,11 +28,67 @@ Group
grpFieldFunctionObjects
Description
Looks up a field, interpolates it to the faces and determines a connected
region from a patch where the field is above a certain value.
- Writes a field containing all regions starting at given patch
('liquid core')
- All other regions are summed for volume and a histogram is calculated.
Droplet size distribution calculation.
Looks up a void-fraction (alpha) field and splits the mesh into regions
based on where the field is below the threshold value. These
regions ("droplets") can now be analysed.
Regions:
- (debug) write regions as a volScalarField
- (debug) print for all regions the sum of volume and alpha*volume
- print the regions connected to a user-defined set of patches.
(in spray calculation these form the liquid core)
- print the regions with too large volume. These are the 'background'
regions.
Fields:
- write volScalarField alpha_liquidCore : alpha with outside liquid core
set to 0.
alpha_background : alpha with outside background
set to 0.
Histogram:
- determine histogram of diameter (given minDiameter, maxDiameter, nBins)
- write graph of number of droplets per bin
- write graph of sum, average and deviation of droplet volume per bin
- write graph of sum, average and deviation of user-defined fields. For
volVectorFields these are those of the 3 components and the magnitude.
Sample input:
functions
{
regionSizeDistribution
{
type regionSizeDistribution;
outputControl timeStep;
outputInterval 1;
// Field to determine regions from
field alpha;
// Patches that provide the liquid core
patches (inlet);
// Delimit alpha regions
threshold 0.4;
// Fields to sample (no need to include alpha)
fields (p U);
// Number of bins for histogram
nBins 100;
// Max droplet diameter
maxDiameter 0.5e-4;
//// Min droplet diameter (default is 0)
//minDiameter 0;
// Writing format
setFormat gnuplot;
}
}
SourceFiles
regionSizeDistribution.C
@ -44,6 +100,9 @@ SourceFiles
#include "pointFieldFwd.H"
#include "writer.H"
#include "Map.H"
#include "volFieldsFwd.H"
#include "wordReList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -54,6 +113,8 @@ namespace Foam
class objectRegistry;
class dictionary;
class mapPolyMesh;
class regionSplit;
class polyMesh;
/*---------------------------------------------------------------------------*\
Class regionSizeDistribution Declaration
@ -75,23 +136,85 @@ class regionSizeDistribution
word alphaName_;
//- Patches to walk from
wordList patchNames_;
wordReList patchNames_;
//- Clip value
scalar threshold_;
//- Background region volFraction
scalar volFraction_;
//- Maximum droplet diameter
scalar maxDiam_;
//- Minimum droplet diameter
scalar minDiam_;
//- Mumber of bins
label nBins_;
//- Names of fields to sample on regions
wordReList fields_;
//- Output formatter to write
autoPtr<writer<scalar> > formatterPtr_;
// Private Member Functions
template<class Type>
Map<Type> regionSum(const regionSplit&, const Field<Type>&) const;
//- Get data in order
template<class Type>
List<Type> extractData(const UList<label>& keys, const Map<Type>&)
const;
void writeGraph
(
const coordSet& coords,
const word& valueName,
const scalarField& values
) const;
//- Write volfields with the parts of alpha which are not
// droplets (liquidCore, backGround)
void writeAlphaFields
(
const regionSplit& regions,
const Map<label>& keepRegions,
const Map<scalar>& regionVolume,
const volScalarField& alpha
) const;
//- Mark all regions starting at patches
Map<label> findPatchRegions(const polyMesh&, const regionSplit&) const;
//- Helper: divide if denom != 0
static tmp<scalarField> divide(const scalarField&, const scalarField&);
//- Given per-region data calculate per-bin average/deviation and graph
void writeGraphs
(
const word& fieldName, // name of field
const labelList& indices, // index of bin for each region
const scalarField& sortedField, // per region field data
const scalarField& binCount, // per bin number of regions
const coordSet& coords // graph data for bins
) const;
//- Given per-cell data calculate per-bin average/deviation and graph
void writeGraphs
(
const word& fieldName, // name of field
const scalarField& cellField, // per cell field data
const regionSplit& regions, // per cell the region(=droplet)
const labelList& sortedRegions, // valid regions in sorted order
const scalarField& sortedNormalisation,
const labelList& indices, // index of bin for each region
const scalarField& binCount, // per bin number of regions
const coordSet& coords // graph data for bins
) const;
//- Disallow default bitwise copy construct
regionSizeDistribution(const regionSizeDistribution&);
@ -159,6 +282,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "regionSizeDistributionTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

81
src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistributionTemplates.C Normal file
View File

@ -0,0 +1,81 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "regionSizeDistribution.H"
#include "regionSplit.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::Map<Type> Foam::regionSizeDistribution::regionSum
(
const regionSplit& regions,
const Field<Type>& fld
) const
{
// Per region the sum of fld
Map<Type> regionToSum(regions.nRegions()/Pstream::nProcs());
forAll(fld, cellI)
{
label regionI = regions[cellI];
typename Map<Type>::iterator fnd = regionToSum.find(regionI);
if (fnd == regionToSum.end())
{
regionToSum.insert(regionI, fld[cellI]);
}
else
{
fnd() += fld[cellI];
}
}
Pstream::mapCombineGather(regionToSum, plusEqOp<Type>());
Pstream::mapCombineScatter(regionToSum);
return regionToSum;
}
// Get data in sortedToc order
template<class Type>
Foam::List<Type> Foam::regionSizeDistribution::extractData
(
const UList<label>& keys,
const Map<Type>& regionData
) const
{
List<Type> sortedData(keys.size());
forAll(keys, i)
{
sortedData[i] = regionData[keys[i]];
}
return sortedData;
}
// ************************************************************************* //