Merge branch 'master' into develop

This commit is contained in:
Andrew Heather
2016-09-08 12:00:46 +01:00
32 changed files with 741 additions and 212 deletions

View File

@ -41,7 +41,6 @@
p_rghDDtEqn = p_rghDDtEqn =
( (
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
== ==
fvOptions(psi, p_rgh, rho.name()) fvOptions(psi, p_rgh, rho.name())
); );
@ -52,6 +51,7 @@
fvScalarMatrix p_rghEqn fvScalarMatrix p_rghEqn
( (
p_rghDDtEqn() p_rghDDtEqn()
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p_rgh) - fvm::laplacian(rhorAUf, p_rgh)
); );

View File

@ -118,6 +118,8 @@ tr.memlist
.OFPlainTable tr td { .OFPlainTable tr td {
height: 20px; height: 20px;
padding-left: 5px; padding-left: 5px;
}
div.line, div.line,
span.comment, span.comment,
span.keyword, span.keyword,

View File

@ -64,7 +64,7 @@ then
# it is either located within ThirdParty, or a central installation # it is either located within ThirdParty, or a central installation
# outside of ThirdParty and must be added to the lib-path. # outside of ThirdParty and must be added to the lib-path.
ending="${FFTW_ARCH_PATH_PATH##*-}" ending="${FFTW_ARCH_PATH##*-}"
if [ "$ending" != none -a "$ending" != system ] if [ "$ending" != none -a "$ending" != system ]
then then
_foamAddLib $FFTW_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH _foamAddLib $FFTW_ARCH_PATH/lib$WM_COMPILER_LIB_ARCH

View File

@ -129,23 +129,29 @@ Foam::plane::plane(const vector& normalVector)
} }
Foam::plane::plane(const point& basePoint, const vector& normalVector) Foam::plane::plane
(
const point& basePoint,
const vector& normalVector,
const bool normalise
)
: :
unitVector_(normalVector), unitVector_(normalVector),
basePoint_(basePoint) basePoint_(basePoint)
{ {
scalar magUnitVector(mag(unitVector_)); scalar magSqrUnitVector(magSqr(unitVector_));
if (magUnitVector > VSMALL) if (magSqrUnitVector < VSMALL)
{
unitVector_ /= magUnitVector;
}
else
{ {
FatalErrorInFunction FatalErrorInFunction
<< "plane normal has zero length. basePoint:" << basePoint_ << "plane normal has zero length. basePoint:" << basePoint_
<< abort(FatalError); << abort(FatalError);
} }
if (normalise)
{
unitVector_ /= Foam::sqrt(magSqrUnitVector);
}
} }

View File

@ -131,7 +131,12 @@ public:
plane(const vector& normalVector); plane(const vector& normalVector);
//- Construct from normal vector and point in plane //- Construct from normal vector and point in plane
plane(const point& basePoint, const vector& normalVector); plane
(
const point& basePoint,
const vector& normalVector,
const bool normalise = true
);
//- Construct from three points in plane //- Construct from three points in plane
plane(const point& point1, const point& point2, const point& point3); plane(const point& point1, const point& point2, const point& point3);

View File

@ -47,26 +47,31 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
label nD = mesh.nGeometricD(); label nD = mesh.nGeometricD();
const cellList& cells = mesh.cells(); const cellList& cells = mesh.cells();
const vectorField& cellC = mesh.cellCentres();
const vectorField& faceC = mesh.faceCentres();
const vectorField faceN(mesh.faceAreas()/mag(mesh.faceAreas()));
scalarField hmax(cells.size()); scalarField hmax(cells.size());
forAll(cells,cellI) forAll(cells, celli)
{ {
scalar deltaMaxTmp = 0.0; scalar deltaMaxTmp = 0.0;
const labelList& cFaces = mesh.cells()[cellI]; const labelList& cFaces = cells[celli];
const point& centrevector = mesh.cellCentres()[cellI]; const point& cc = cellC[celli];
forAll(cFaces, cFaceI) forAll(cFaces, cFacei)
{ {
label faceI = cFaces[cFaceI]; label facei = cFaces[cFacei];
const point& facevector = mesh.faceCentres()[faceI]; const point& fc = faceC[facei];
scalar tmp = mag(facevector - centrevector); const vector& n = faceN[facei];
scalar tmp = magSqr(n*(n & (fc - cc)));
if (tmp > deltaMaxTmp) if (tmp > deltaMaxTmp)
{ {
deltaMaxTmp = tmp; deltaMaxTmp = tmp;
} }
} }
hmax[cellI] = deltaCoeff_*deltaMaxTmp; hmax[celli] = deltaCoeff_*Foam::sqrt(deltaMaxTmp);
} }
if (nD == 3) if (nD == 3)
@ -76,7 +81,7 @@ void Foam::LESModels::maxDeltaxyz::calcDelta()
else if (nD == 2) else if (nD == 2)
{ {
WarningInFunction WarningInFunction
<< "Case is 2D, LES is not strictly applicable\n" << "Case is 2D, LES is not strictly applicable" << nl
<< endl; << endl;
delta_.internalField() = hmax; delta_.internalField() = hmax;

View File

@ -332,6 +332,20 @@ void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs()
areaFraction = 1 - openFraction_; areaFraction = 1 - openFraction_;
} }
if (patch().boundaryMesh().mesh().moving())
{
// All areas have been recalculated already so are consistent
// with the new points. Note we already only do this routine
// once per timestep. The alternative is to use the upToDate
// mechanism for regIOobject + polyMesh::upToDatePoints
initWallSf_ = patch().Sf();
initCyclicSf_ = patch().boundaryMesh()[cyclicPatchLabel_].Sf();
nbrCyclicSf_ = refCast<const cyclicFvPatch>
(
patch().boundaryMesh()[cyclicPatchLabel_]
).neighbFvPatch().Sf();
}
// Update this wall patch // Update this wall patch
vectorField::subField Sfw = patch().patch().faceAreas(); vectorField::subField Sfw = patch().patch().faceAreas();
vectorField newSfw((1 - areaFraction)*initWallSf_); vectorField newSfw((1 - areaFraction)*initWallSf_);

View File

@ -48,7 +48,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
: :
uniformJumpFvPatchField<Type>(p, iF), uniformJumpFvPatchField<Type>(p, iF),
phiName_("phi"), phiName_("phi"),
rhoName_("rho") rhoName_("rho"),
uniformJump_(false)
{} {}
@ -62,7 +63,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
: :
uniformJumpFvPatchField<Type>(p, iF, dict), uniformJumpFvPatchField<Type>(p, iF, dict),
phiName_(dict.lookupOrDefault<word>("phi", "phi")), phiName_(dict.lookupOrDefault<word>("phi", "phi")),
rhoName_(dict.lookupOrDefault<word>("rho", "rho")) rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
uniformJump_(dict.lookupOrDefault<bool>("uniformJump", "false"))
{} {}
@ -77,7 +79,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
: :
uniformJumpFvPatchField<Type>(ptf, p, iF, mapper), uniformJumpFvPatchField<Type>(ptf, p, iF, mapper),
phiName_(ptf.phiName_), phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_) rhoName_(ptf.rhoName_),
uniformJump_(ptf.uniformJump_)
{} {}
@ -89,7 +92,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
: :
uniformJumpFvPatchField<Type>(ptf), uniformJumpFvPatchField<Type>(ptf),
phiName_(ptf.phiName_), phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_) rhoName_(ptf.rhoName_),
uniformJump_(ptf.uniformJump_)
{} {}
@ -102,7 +106,8 @@ Foam::fanFvPatchField<Type>::fanFvPatchField
: :
uniformJumpFvPatchField<Type>(ptf, iF), uniformJumpFvPatchField<Type>(ptf, iF),
phiName_(ptf.phiName_), phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_) rhoName_(ptf.rhoName_),
uniformJump_(ptf.uniformJump_)
{} {}
@ -129,6 +134,10 @@ void Foam::fanFvPatchField<Type>::write(Ostream& os) const
uniformJumpFvPatchField<Type>::write(os); uniformJumpFvPatchField<Type>::write(os);
this->template writeEntryIfDifferent<word>(os, "phi", "phi", phiName_); this->template writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
this->template writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_); this->template writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
this->template writeEntryIfDifferent<bool>
(
os, "uniformJump", false, uniformJump_
);
} }

View File

@ -42,6 +42,8 @@ Description
jumpTable | jump data, e.g. \c csvFile | yes | jumpTable | jump data, e.g. \c csvFile | yes |
phi | flux field name | no | phi phi | flux field name | no | phi
rho | density field name | no | none rho | density field name | no | none
uniformJump | applies a uniform pressure based on the averaged
velocity | no | false
\endtable \endtable
Example of the boundary condition specification: Example of the boundary condition specification:
@ -53,10 +55,11 @@ Description
jumpTable csvFile; jumpTable csvFile;
csvFileCoeffs csvFileCoeffs
{ {
hasHeaderLine 1; nHeaderLine 1;
refColumn 0; refColumn 0;
componentColumns 1(1); componentColumns 1(1);
separator ","; separator ",";
mergeSeparators no;
fileName "$FOAM_CASE/constant/pressureVsU"; fileName "$FOAM_CASE/constant/pressureVsU";
} }
value uniform 0; value uniform 0;
@ -67,7 +70,7 @@ Description
the jump condition. the jump condition.
Note Note
The underlying \c patchType should be set to \c cyclic The underlying \c patchType should be set to \c cyclic
SeeAlso SeeAlso
Foam::Function1Types Foam::Function1Types
@ -109,6 +112,9 @@ class fanFvPatchField
// if neccessary // if neccessary
word rhoName_; word rhoName_;
//- Uniform pressure drop
bool uniformJump_;
// Private Member Functions // Private Member Functions

View File

@ -44,6 +44,11 @@ void Foam::fanFvPatchField<Foam::scalar>::calcFanJump()
patch().patchField<surfaceScalarField, scalar>(phi); patch().patchField<surfaceScalarField, scalar>(phi);
scalarField Un(max(phip/patch().magSf(), scalar(0))); scalarField Un(max(phip/patch().magSf(), scalar(0)));
if (uniformJump_)
{
scalar area = gSum(patch().magSf());
Un = gSum(Un*patch().magSf())/area;
}
if (phi.dimensions() == dimDensity*dimVelocity*dimArea) if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
{ {
@ -67,7 +72,8 @@ Foam::fanFvPatchField<Foam::scalar>::fanFvPatchField
: :
uniformJumpFvPatchField<scalar>(p, iF), uniformJumpFvPatchField<scalar>(p, iF),
phiName_(dict.lookupOrDefault<word>("phi", "phi")), phiName_(dict.lookupOrDefault<word>("phi", "phi")),
rhoName_(dict.lookupOrDefault<word>("rho", "rho")) rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
uniformJump_(dict.lookupOrDefault<bool>("uniformJump", false))
{ {
if (this->cyclicPatch().owner()) if (this->cyclicPatch().owner())
{ {

View File

@ -848,7 +848,7 @@ turbulentDFSEMInletFvPatchVectorField
eddy::debug = debug; eddy::debug = debug;
// Set UMean as patch area average value // Set UMean as patch area average value
UMean_ = gSum(U_*patch().magSf())/gSum(patch().magSf()); UMean_ = gSum(U_*patch().magSf())/(gSum(patch().magSf()) + ROOTVSMALL);
} }

View File

@ -40,6 +40,13 @@ namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
const convectionScheme<Type>& boundedConvectionScheme<Type>::scheme() const
{
return scheme_();
}
template<class Type> template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
boundedConvectionScheme<Type>::interpolate boundedConvectionScheme<Type>::interpolate

View File

@ -106,6 +106,8 @@ public:
// Member Functions // Member Functions
const convectionScheme<Type>& scheme() const;
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
( (
const surfaceScalarField&, const surfaceScalarField&,

View File

@ -39,12 +39,6 @@ namespace Foam
defineTypeNameAndDebug(tetOverlapVolume, 0); defineTypeNameAndDebug(tetOverlapVolume, 0);
} }
// When to consider a tet to be zero volume. We want to avoid doing clipping
// against negative volume tets. Tet volume can be calculated incorrectly
// due to truncation errors. The value below works for single and double
// precision but could probably be tighter for double precision.
Foam::scalar Foam::tetOverlapVolume::minTetVolume_ = SMALL*SMALL;
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

View File

@ -54,12 +54,6 @@ class polyMesh;
class tetOverlapVolume class tetOverlapVolume
{ {
// Private data
//- Minimum tet volume to skip test
static scalar minTetVolume_;
// Private classes // Private classes
//- tetPoints handling : sum resulting volumes //- tetPoints handling : sum resulting volumes

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -45,51 +45,90 @@ void Foam::tetOverlapVolume::tetTetOverlap
tetPointRef::storeOp cutInside(cutInsideTets, nCutInside); tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
tetPointRef::dummyOp outside; tetPointRef::dummyOp outside;
if (tetA.tet().mag() < minTetVolume_ || tetB.tet().mag() < minTetVolume_) // Precompute the tet face areas and exit early if any face area is
// too small
static FixedList<vector, 4> tetAFaceAreas;
static FixedList<scalar, 4> tetAMag2FaceAreas;
tetPointRef tetATet = tetA.tet();
for (label facei = 0; facei < 4; ++facei)
{ {
return; tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
{
return;
}
} }
// face0 static FixedList<vector, 4> tetBFaceAreas;
plane pl0(tetB[1], tetB[3], tetB[2]); static FixedList<scalar, 4> tetBMag2FaceAreas;
tetA.tet().sliceWithPlane(pl0, cutInside, outside); tetPointRef tetBTet = tetB.tet();
if (nCutInside == 0) for (label facei = 0; facei < 4; ++facei)
{ {
return; tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
{
return;
}
} }
// face1
plane pl1(tetB[0], tetB[2], tetB[3]); // Face 0
nInside = 0;
for (label i = 0; i < nCutInside; i++)
{ {
const tetPointRef t = cutInsideTets[i].tet(); vector n = tetBFaceAreas[0]/Foam::sqrt(tetBMag2FaceAreas[0]);
t.sliceWithPlane(pl1, inside, outside); plane pl0(tetBTet.tri(0).a(), n, false);
}
if (nInside == 0) tetA.tet().sliceWithPlane(pl0, cutInside, outside);
{ if (nCutInside == 0)
return; {
return;
}
} }
// face2 // Face 1
plane pl2(tetB[0], tetB[3], tetB[1]);
nCutInside = 0;
for (label i = 0; i < nInside; i++)
{ {
const tetPointRef t = insideTets[i].tet(); vector n = tetBFaceAreas[1]/Foam::sqrt(tetBMag2FaceAreas[1]);
t.sliceWithPlane(pl2, cutInside, outside); plane pl1(tetBTet.tri(1).a(), n, false);
}
if (nCutInside == 0) nInside = 0;
{ for (label i = 0; i < nCutInside; i++)
return; {
const tetPointRef t = cutInsideTets[i].tet();
t.sliceWithPlane(pl1, inside, outside);
}
if (nInside == 0)
{
return;
}
} }
// face3 // Face 2
plane pl3(tetB[0], tetB[1], tetB[2]);
for (label i = 0; i < nCutInside; i++)
{ {
const tetPointRef t = cutInsideTets[i].tet(); vector n = tetBFaceAreas[2]/Foam::sqrt(tetBMag2FaceAreas[2]);
t.sliceWithPlane(pl3, insideOp, outside); plane pl2(tetBTet.tri(2).a(), n, false);
nCutInside = 0;
for (label i = 0; i < nInside; i++)
{
const tetPointRef t = insideTets[i].tet();
t.sliceWithPlane(pl2, cutInside, outside);
}
if (nCutInside == 0)
{
return;
}
}
// Face 3
{
vector n = tetBFaceAreas[3]/Foam::sqrt(tetBMag2FaceAreas[3]);
plane pl3(tetBTet.tri(3).a(), n, false);
for (label i = 0; i < nCutInside; i++)
{
const tetPointRef t = cutInsideTets[i].tet();
t.sliceWithPlane(pl3, insideOp, outside);
}
} }
} }

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -335,7 +335,8 @@ Foam::regionSizeDistribution::regionSizeDistribution
active_(true), active_(true),
alphaName_(dict.lookup("field")), alphaName_(dict.lookup("field")),
patchNames_(dict.lookup("patches")), patchNames_(dict.lookup("patches")),
log_(true) log_(true),
isoPlanes_(dict.lookupOrDefault<bool>("isoPlanes", false))
{ {
// Check if the available mesh is an fvMesh, otherwise deactivate // Check if the available mesh is an fvMesh, otherwise deactivate
if (isA<fvMesh>(obr_)) if (isA<fvMesh>(obr_))
@ -388,6 +389,17 @@ void Foam::regionSizeDistribution::read(const dictionary& dict)
<< "Transforming all vectorFields with coordinate system " << "Transforming all vectorFields with coordinate system "
<< coordSysPtr_().name() << endl; << coordSysPtr_().name() << endl;
} }
if (isoPlanes_)
{
dict.lookup("origin") >> origin_;
dict.lookup("direction") >> direction_;
dict.lookup("maxDiameter") >> maxDiameter_;
dict.lookup("nDownstreamBins") >> nDownstreamBins_;
dict.lookup("maxDownstream") >> maxDownstream_;
direction_ /= mag(direction_);
}
} }
} }
@ -691,6 +703,8 @@ void Foam::regionSizeDistribution::write()
// (patchRegions) and background regions from maps. // (patchRegions) and background regions from maps.
// Note that we have to use mesh volume (allRegionVolume) and not // Note that we have to use mesh volume (allRegionVolume) and not
// allRegionAlphaVolume since background might not have alpha in it. // allRegionAlphaVolume since background might not have alpha in it.
// Deleting regions where the volume-alpha-weighted is lower than
// threshold
forAllIter(Map<scalar>, allRegionVolume, vIter) forAllIter(Map<scalar>, allRegionVolume, vIter)
{ {
label regionI = vIter.key(); label regionI = vIter.key();
@ -698,6 +712,7 @@ void Foam::regionSizeDistribution::write()
( (
patchRegions.found(regionI) patchRegions.found(regionI)
|| vIter() >= maxDropletVol || vIter() >= maxDropletVol
|| (allRegionAlphaVolume[regionI]/vIter() < threshold_)
) )
{ {
allRegionVolume.erase(vIter); allRegionVolume.erase(vIter);
@ -735,6 +750,108 @@ void Foam::regionSizeDistribution::write()
) )
); );
vectorField centroids(sortedVols.size(), vector::zero);
// Check if downstream bins are calculed
if (isoPlanes_)
{
vectorField alphaDistance
(
(alpha.internalField()*mesh.V())
* (mesh.C().internalField() - origin_)()
);
Map<vector> allRegionAlphaDistance
(
regionSum
(
regions,
alphaDistance
)
);
// 2. centroid
vectorField sortedMoment
(
extractData
(
sortedRegions,
allRegionAlphaDistance
)
);
centroids = sortedMoment/sortedVols + origin_;
// Bin according to centroid
scalarField distToPlane((centroids - origin_) & direction_);
vectorField radialDistToOrigin
(
(centroids - origin_) - (distToPlane*direction_)
);
const scalar deltaX = maxDownstream_/nDownstreamBins_;
labelList downstreamIndices(distToPlane.size(), -1);
forAll(distToPlane, i)
{
if
(
(mag(radialDistToOrigin[i]) < maxDiameter_)
&& (distToPlane[i] < maxDownstream_)
)
{
downstreamIndices[i] = distToPlane[i]/deltaX;
}
}
scalarField binDownCount(nDownstreamBins_, 0.0);
forAll(distToPlane, i)
{
if (downstreamIndices[i] != -1)
{
binDownCount[downstreamIndices[i]] += 1.0;
}
}
// Write
if (Pstream::master())
{
// Construct mids of bins for plotting
pointField xBin(nDownstreamBins_);
scalar x = 0.5*deltaX;
forAll(xBin, i)
{
xBin[i] = point(x, 0, 0);
x += deltaX;
}
const coordSet coords("distance", "x", xBin, mag(xBin));
writeGraph(coords, "isoPlanes", binDownCount);
}
// Write to screen
if (log_)
{
Info<< " Iso-planes Bins:" << endl;
Info<< " " << token::TAB << "Bin"
<< token::TAB << "Min distance"
<< token::TAB << "Count:"
<< endl;
scalar delta = 0.0;
forAll(binDownCount, binI)
{
Info<< " " << token::TAB << binI
<< token::TAB << delta
<< token::TAB << binDownCount[binI] << endl;
delta += deltaX;
}
Info<< endl;
}
}
// Calculate the diameters // Calculate the diameters
scalarField sortedDiameters(sortedVols.size()); scalarField sortedDiameters(sortedVols.size());
forAll(sortedDiameters, i) forAll(sortedDiameters, i)

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -35,6 +35,8 @@ Description
based on where the field is below the threshold value. These based on where the field is below the threshold value. These
regions ("droplets") can now be analysed. regions ("droplets") can now be analysed.
Regions: Regions:
- print the regions connected to a user-defined set of patches. - print the regions connected to a user-defined set of patches.
(in spray calculation these form the liquid core) (in spray calculation these form the liquid core)
@ -53,6 +55,9 @@ Description
- write graph of sum, average and deviation of droplet volume 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 - write graph of sum, average and deviation of user-defined fields. For
volVectorFields these are those of the 3 components and the magnitude. volVectorFields these are those of the 3 components and the magnitude.
- (optional) write graph of histogram of centroids on iso planes
downstream of the injector determined by origin, direction and maxDiameter
up to maxDownstream
Example of function object specification: Example of function object specification:
\verbatim \verbatim
@ -76,6 +81,24 @@ Description
e3 (0 1 1); e3 (0 1 1);
e1 (1 0 0); e1 (1 0 0);
} }
// Optional downstream iso-plane bins.
isoPlanes true;
// Plane normal and point definition
direction (1 0 1);
origin (1e-4 0 5e-4);
// Maximum diameter of the cylinder formed by the origin point
// and direction
maxDiameter 3e-4;
// Maximum downstream distance
maxDownstream 6e-4;
// Number of iso-plane bins
nDownstreamBins 20;
} }
\endverbatim \endverbatim
@ -177,6 +200,27 @@ class regionSizeDistribution
//- Optional coordinate system //- Optional coordinate system
autoPtr<coordinateSystem> coordSysPtr_; autoPtr<coordinateSystem> coordSysPtr_;
// Optional extra definition of bins on planes downstream to the origin
// point and maximum diameter
//- Switch to enable iso-planes sampling
bool isoPlanes_;
//- Optional origin point
vector origin_;
//- Optional plane normal direction
vector direction_;
//- Optional maximum diameter on plane
scalar maxDiameter_;
//- Optional number of bins for
scalar nDownstreamBins_;
//- Optional maximum downstream coordinate from origin
scalar maxDownstream_;
// Private Member Functions // Private Member Functions

View File

@ -43,6 +43,7 @@ License
#include "vtkSphereSource.h" #include "vtkSphereSource.h"
#include "vtkTextActor.h" #include "vtkTextActor.h"
#include "vtkTextProperty.h" #include "vtkTextProperty.h"
#include "vtkCellDataToPointData.h"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@ -151,6 +152,9 @@ void Foam::fieldVisualisationBase::addScalarBar
const vector textColour = colours_["text"]->value(position); const vector textColour = colours_["text"]->value(position);
// Work-around to supply our own scalarbar title // Work-around to supply our own scalarbar title
// - Default scalar bar title text is scales by the scalar bar box
// dimensions so if the title is a long string, the text is shrunk to fit
// Instead, suppress title and set the title using a vtkTextActor
vtkSmartPointer<vtkTextActor> titleActor = vtkSmartPointer<vtkTextActor> titleActor =
vtkSmartPointer<vtkTextActor>::New(); vtkSmartPointer<vtkTextActor>::New();
sbar->SetTitle(" "); sbar->SetTitle(" ");
@ -170,19 +174,18 @@ void Foam::fieldVisualisationBase::addScalarBar
titleActor->GetPositionCoordinate()-> titleActor->GetPositionCoordinate()->
SetCoordinateSystemToNormalizedViewport(); SetCoordinateSystemToNormalizedViewport();
/* // How to use the standard scalar bar text
sbar->SetTitle(scalarBar_.title_.c_str()); // sbar->SetTitle(scalarBar_.title_.c_str());
sbar->GetTitleTextProperty()->SetColor // sbar->GetTitleTextProperty()->SetColor
( // (
textColour[0], // textColour[0],
textColour[1], // textColour[1],
textColour[2] // textColour[2]
); // );
sbar->GetTitleTextProperty()->SetFontSize(scalarBar_.fontSize_); // sbar->GetTitleTextProperty()->SetFontSize(scalarBar_.fontSize_);
sbar->GetTitleTextProperty()->ShadowOff(); // sbar->GetTitleTextProperty()->ShadowOff();
sbar->GetTitleTextProperty()->BoldOn(); // sbar->GetTitleTextProperty()->BoldOn();
sbar->GetTitleTextProperty()->ItalicOff(); // sbar->GetTitleTextProperty()->ItalicOff();
*/
sbar->GetLabelTextProperty()->SetColor sbar->GetLabelTextProperty()->SetColor
( (
@ -217,8 +220,8 @@ void Foam::fieldVisualisationBase::addScalarBar
sbar->SetWidth(0.75); sbar->SetWidth(0.75);
sbar->SetHeight(0.07); sbar->SetHeight(0.07);
sbar->SetBarRatio(0.5); sbar->SetBarRatio(0.5);
// sbar->SetHeight(0.1); // sbar->SetHeight(0.1);
// sbar->SetTitleRatio(0.01); // sbar->SetTitleRatio(0.01);
sbar->SetTextPositionToPrecedeScalarBar(); sbar->SetTextPositionToPrecedeScalarBar();
} }
@ -228,10 +231,10 @@ void Foam::fieldVisualisationBase::addScalarBar
scalarBar_.position_.second() + sbar->GetHeight() scalarBar_.position_.second() + sbar->GetHeight()
); );
// sbar->DrawFrameOn(); // sbar->DrawFrameOn();
// sbar->DrawBackgroundOn(); // sbar->DrawBackgroundOn();
// sbar->UseOpacityOff(); // sbar->UseOpacityOff();
// sbar->VisibilityOff(); // sbar->VisibilityOff();
sbar->VisibilityOn(); sbar->VisibilityOn();
renderer->AddActor(sbar); renderer->AddActor(sbar);
@ -266,19 +269,21 @@ void Foam::fieldVisualisationBase::setField
lut->SetVectorMode(vtkScalarsToColors::MAGNITUDE); lut->SetVectorMode(vtkScalarsToColors::MAGNITUDE);
// Configure the mapper // Configure the mapper
mapper->SelectColorArray(colourFieldName.c_str());
mapper->SetScalarRange(range_.first(), range_.second());
// Set to use either cell or point data
const char* fieldName = colourFieldName.c_str(); const char* fieldName = colourFieldName.c_str();
if (pData->GetCellData()->HasArray(fieldName) == 1) mapper->SelectColorArray(fieldName);
{
mapper->SetScalarModeToUseCellFieldData(); // Set to use either point or cell data
} // Note: if both point and cell data exists, preferentially
else if (pData->GetPointData()->HasArray(fieldName) == 1) // choosing point data. This is often the case when using
// glyphs
if (pData->GetPointData()->HasArray(fieldName) == 1)
{ {
mapper->SetScalarModeToUsePointFieldData(); mapper->SetScalarModeToUsePointFieldData();
} }
else if (pData->GetCellData()->HasArray(fieldName) == 1)
{
mapper->SetScalarModeToUseCellFieldData();
}
else else
{ {
WarningInFunction WarningInFunction
@ -286,7 +291,7 @@ void Foam::fieldVisualisationBase::setField
<< "- assuming point data"; << "- assuming point data";
mapper->SetScalarModeToUsePointFieldData(); mapper->SetScalarModeToUsePointFieldData();
} }
mapper->SetScalarRange(range_.first(), range_.second());
mapper->SetColorModeToMapScalars(); mapper->SetColorModeToMapScalars();
mapper->SetLookupTable(lut); mapper->SetLookupTable(lut);
mapper->ScalarVisibilityOn(); mapper->ScalarVisibilityOn();
@ -322,9 +327,37 @@ void Foam::fieldVisualisationBase::addGlyphs
glyph->ScalingOn(); glyph->ScalingOn();
bool ok = true; bool ok = true;
label nComponents = // Determine whether we have scalar or vector data
data->GetPointData()->GetArray(scaleFieldName.c_str()) label nComponents = -1;
->GetNumberOfComponents(); const char* scaleFieldNameChar = scaleFieldName.c_str();
if (data->GetPointData()->HasArray(scaleFieldNameChar) == 1)
{
nComponents =
data->GetPointData()->GetArray(scaleFieldNameChar)
->GetNumberOfComponents();
}
else if (data->GetCellData()->HasArray(scaleFieldNameChar) == 1)
{
// Need to convert cell data to point data
vtkSmartPointer<vtkCellDataToPointData> cellToPoint =
vtkSmartPointer<vtkCellDataToPointData>::New();
cellToPoint->SetInputData(data);
cellToPoint->Update();
vtkDataSet* pds = cellToPoint->GetOutput();
vtkDataArray* pData = pds->GetPointData()->GetArray(scaleFieldNameChar);
// Store in main vtkPolyData
data->GetPointData()->AddArray(pData);
nComponents = pData->GetNumberOfComponents();
}
else
{
WarningInFunction
<< "Glyphs can only be added to scalar or vector data. "
<< "Unable to process field " << scaleFieldName << endl;
return;
}
if (nComponents == 1) if (nComponents == 1)
{ {
@ -332,9 +365,10 @@ void Foam::fieldVisualisationBase::addGlyphs
vtkSmartPointer<vtkSphereSource>::New(); vtkSmartPointer<vtkSphereSource>::New();
sphere->SetCenter(0, 0, 0); sphere->SetCenter(0, 0, 0);
sphere->SetRadius(0.5); sphere->SetRadius(0.5);
// Setting higher resolution slows the rendering significantly
// sphere->SetPhiResolution(20); // Setting higher resolution slows the rendering significantly
// sphere->SetThetaResolution(20); // sphere->SetPhiResolution(20);
// sphere->SetThetaResolution(20);
glyph->SetSourceConnection(sphere->GetOutputPort()); glyph->SetSourceConnection(sphere->GetOutputPort());
@ -342,18 +376,18 @@ void Foam::fieldVisualisationBase::addGlyphs
{ {
double range[2]; double range[2];
// Can use values to find range // Can use values to find range
// vtkDataArray* values = // vtkDataArray* values =
// data->GetPointData()->GetScalars(scaleFieldName.c_str()); // data->GetPointData()->GetScalars(scaleFieldNameChar);
// values->GetRange(range); // values->GetRange(range);
// set range accoding to user-supplied limits // Set range accoding to user-supplied limits
range[0] = range_.first(); range[0] = range_.first();
range[1] = range_.second(); range[1] = range_.second();
glyph->ClampingOn(); glyph->ClampingOn();
glyph->SetRange(range); glyph->SetRange(range);
// if range[0] != min(value), maxGlyphLength behaviour will not // If range[0] != min(value), maxGlyphLength behaviour will not
// be correct... // be correct...
glyph->SetScaleFactor(maxGlyphLength); glyph->SetScaleFactor(maxGlyphLength);
} }
@ -370,7 +404,7 @@ void Foam::fieldVisualisationBase::addGlyphs
0, 0,
0, 0,
vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataObject::FIELD_ASSOCIATION_POINTS,
scaleFieldName.c_str() scaleFieldNameChar
); );
} }
else if (nComponents == 3) else if (nComponents == 3)
@ -388,21 +422,21 @@ void Foam::fieldVisualisationBase::addGlyphs
if (maxGlyphLength > 0) if (maxGlyphLength > 0)
{ {
vtkDataArray* values = vtkDataArray* values =
data->GetPointData()->GetVectors(scaleFieldName.c_str()); data->GetPointData()->GetVectors(scaleFieldNameChar);
double range[6]; double range[6];
values->GetRange(range); values->GetRange(range);
/*
// Attempt to set range for vectors... // Attempt to set range for vectors...
scalar x0 = sqrt(sqr(range_.first())/3.0); // scalar x0 = sqrt(sqr(range_.first())/3.0);
scalar x1 = sqrt(sqr(range_.second())/3.0); // scalar x1 = sqrt(sqr(range_.second())/3.0);
range[0] = x0; // range[0] = x0;
range[1] = x0; // range[1] = x0;
range[2] = x0; // range[2] = x0;
range[3] = x1; // range[3] = x1;
range[4] = x1; // range[4] = x1;
range[5] = x1; // range[5] = x1;
*/
glyph->ClampingOn(); glyph->ClampingOn();
glyph->SetRange(range); glyph->SetRange(range);
glyph->SetScaleFactor(maxGlyphLength); glyph->SetScaleFactor(maxGlyphLength);
@ -421,7 +455,7 @@ void Foam::fieldVisualisationBase::addGlyphs
0, 0,
0, 0,
vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataObject::FIELD_ASSOCIATION_POINTS,
scaleFieldName.c_str() scaleFieldNameChar
); );
} }
else else

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -88,6 +88,8 @@ SourceFiles
#include "OFstream.H" #include "OFstream.H"
#include "Switch.H" #include "Switch.H"
#include "convectionScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam namespace Foam
@ -140,6 +142,14 @@ class blendingFactor
//- Disallow default bitwise assignment //- Disallow default bitwise assignment
void operator=(const blendingFactor&); void operator=(const blendingFactor&);
//- Helper function to calculate the blending factor for the scheme
template<class Type>
void calcScheme
(
const GeometricField<Type, fvPatchField, volMesh>& field,
const typename fv::convectionScheme<Type>& cs
);
//- Calculate the blending factor //- Calculate the blending factor
template<class Type> template<class Type>
void calc(); void calc();

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. \\/ M anipulation | Copyright (C) 2015-2016 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "gaussConvectionScheme.H" #include "gaussConvectionScheme.H"
#include "boundedConvectionScheme.H"
#include "blendedSchemeBase.H" #include "blendedSchemeBase.H"
#include "fvcCellReduce.H" #include "fvcCellReduce.H"
#include "zeroGradientFvPatchFields.H" #include "zeroGradientFvPatchFields.H"
@ -31,39 +32,35 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type> template<class Type>
void Foam::blendingFactor::calc() void Foam::blendingFactor::calcScheme
(
const GeometricField<Type, fvPatchField, volMesh>& field,
const typename fv::convectionScheme<Type>& cs
)
{ {
typedef GeometricField<Type, fvPatchField, volMesh> fieldType; if (!isA<fv::gaussConvectionScheme<Type>>(cs))
if (!obr_.foundObject<fieldType>(fieldName_))
{ {
WarningInFunction
<< "Scheme for field " << field.name() << " is not a "
<< fv::gaussConvectionScheme<Type>::typeName
<< " scheme. Not calculating " << resultName_ << endl;
return; return;
} }
const fvMesh& mesh = refCast<const fvMesh>(obr_);
const fieldType& field = mesh.lookupObject<fieldType>(fieldName_);
const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
ITstream& its = mesh.divScheme(divScheme);
const surfaceScalarField& phi =
mesh.lookupObject<surfaceScalarField>(phiName_);
tmp<fv::convectionScheme<Type>> cs =
fv::convectionScheme<Type>::New(mesh, phi, its);
const fv::gaussConvectionScheme<Type>& gcs = const fv::gaussConvectionScheme<Type>& gcs =
refCast<const fv::gaussConvectionScheme<Type>>(cs()); refCast<const fv::gaussConvectionScheme<Type>>(cs);
const surfaceInterpolationScheme<Type>& interpScheme =
gcs.interpScheme(); const surfaceInterpolationScheme<Type>& interpScheme = gcs.interpScheme();
if (!isA<blendedSchemeBase<Type>>(interpScheme)) if (!isA<blendedSchemeBase<Type>>(interpScheme))
{ {
FatalErrorInFunction WarningInFunction
<< interpScheme.typeName << " is not a blended scheme" << interpScheme.type() << " is not a blended scheme"
<< exit(FatalError); << ". Not calculating " << resultName_ << endl;
return;
} }
// Retrieve the face-based blending factor // Retrieve the face-based blending factor
@ -123,4 +120,46 @@ void Foam::blendingFactor::calc()
} }
template<class Type>
void Foam::blendingFactor::calc()
{
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
if (!obr_.foundObject<fieldType>(fieldName_))
{
return;
}
const fvMesh& mesh = refCast<const fvMesh>(obr_);
const fieldType& field = mesh.lookupObject<fieldType>(fieldName_);
const word divScheme("div(" + phiName_ + ',' + fieldName_ + ')');
ITstream& its = mesh.divScheme(divScheme);
const surfaceScalarField& phi =
mesh.lookupObject<surfaceScalarField>(phiName_);
tmp<fv::convectionScheme<Type>> tcs =
fv::convectionScheme<Type>::New(mesh, phi, its);
const fv::convectionScheme<Type>& cs = tcs();
if (isA<fv::boundedConvectionScheme<Type>>(cs))
{
const fv::boundedConvectionScheme<Type>& bcs =
refCast<const fv::boundedConvectionScheme<Type>>(cs);
calcScheme(field, bcs.scheme());
}
else
{
const fv::gaussConvectionScheme<Type>& gcs =
refCast<const fv::gaussConvectionScheme<Type>>(cs);
calcScheme(field, gcs);
}
}
// ************************************************************************* // // ************************************************************************* //

View File

@ -122,6 +122,13 @@ class mapFieldsFO
//- Helper function to create the mesh-to-mesh interpolation //- Helper function to create the mesh-to-mesh interpolation
void createInterpolation(const dictionary& dict); void createInterpolation(const dictionary& dict);
//- Helper function to evaluate constraint patches after mapping
template<class Type>
void evaluateConstraintTypes
(
GeometricField<Type, fvPatchField, volMesh>& fld
) const;
//- Helper function to interpolate and write the fied //- Helper function to interpolate and write the fied
template<class Type> template<class Type>
bool writeFieldType() const; bool writeFieldType() const;

View File

@ -25,6 +25,91 @@ License
#include "meshToMesh.H" #include "meshToMesh.H"
template<class Type>
void Foam::mapFieldsFO::evaluateConstraintTypes
(
GeometricField<Type, fvPatchField, volMesh>& fld
) const
{
typename GeometricField<Type, fvPatchField, volMesh>::
GeometricBoundaryField& fldBf = fld.boundaryField();
if
(
Pstream::defaultCommsType == Pstream::blocking
|| Pstream::defaultCommsType == Pstream::nonBlocking
)
{
label nReq = Pstream::nRequests();
forAll(fldBf, patchi)
{
fvPatchField<Type>& tgtField = fldBf[patchi];
if
(
tgtField.type() == tgtField.patch().patch().type()
&& polyPatch::constraintType(tgtField.patch().patch().type())
)
{
tgtField.initEvaluate(Pstream::defaultCommsType);
}
}
// Block for any outstanding requests
if
(
Pstream::parRun()
&& Pstream::defaultCommsType == Pstream::nonBlocking
)
{
Pstream::waitRequests(nReq);
}
forAll(fldBf, patchi)
{
fvPatchField<Type>& tgtField = fldBf[patchi];
if
(
tgtField.type() == tgtField.patch().patch().type()
&& polyPatch::constraintType(tgtField.patch().patch().type())
)
{
tgtField.evaluate(Pstream::defaultCommsType);
}
}
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
const lduSchedule& patchSchedule =
fld.mesh().globalData().patchSchedule();
forAll(patchSchedule, patchEvali)
{
label patchi = patchSchedule[patchEvali].patch;
fvPatchField<Type>& tgtField = fldBf[patchi];
if
(
tgtField.type() == tgtField.patch().patch().type()
&& polyPatch::constraintType(tgtField.patch().patch().type())
)
{
if (patchSchedule[patchEvali].init)
{
tgtField.initEvaluate(Pstream::scheduled);
}
else
{
tgtField.evaluate(Pstream::scheduled);
}
}
}
}
}
template<class Type> template<class Type>
bool Foam::mapFieldsFO::writeFieldType() const bool Foam::mapFieldsFO::writeFieldType() const
{ {
@ -53,6 +138,9 @@ bool Foam::mapFieldsFO::writeFieldType() const
if (log_) Info<< ": interpolated"; if (log_) Info<< ": interpolated";
FieldType fieldMapRegion(mapRegionIO, tfieldMapRegion); FieldType fieldMapRegion(mapRegionIO, tfieldMapRegion);
evaluateConstraintTypes(fieldMapRegion);
fieldMapRegion.write(); fieldMapRegion.write();
if (log_) Info<< " and written" << nl; if (log_) Info<< " and written" << nl;

View File

@ -127,10 +127,13 @@ void Foam::noiseFFT::octaveBandInfo
fBandIDs = bandIDs.toc(); fBandIDs = bandIDs.toc();
// Remove the last centre frequency (beyond upper frequency limit) if (fc.size())
fc.remove(); {
// Remove the last centre frequency (beyond upper frequency limit)
fc.remove();
fCentre.transfer(fc); fCentre.transfer(fc);
}
} }

View File

@ -494,8 +494,21 @@ void surfaceNoise::calculate()
octave13FreqCentre octave13FreqCentre
); );
List<scalarField> surfPSD13f(octave13BandIDs.size() - 1); label bandSize = 0;
List<scalarField> surfPrms13f2(octave13BandIDs.size() - 1); if (octave13BandIDs.empty())
{
WarningInFunction
<< "Ocatve band calculation failed (zero sized). "
<< "please check your input data"
<< endl;
}
else
{
bandSize = octave13BandIDs.size() - 1;
}
List<scalarField> surfPSD13f(bandSize);
List<scalarField> surfPrms13f2(bandSize);
forAll(surfPSD13f, freqI) forAll(surfPSD13f, freqI)
{ {
surfPSD13f[freqI].setSize(nLocalFace); surfPSD13f[freqI].setSize(nLocalFace);

View File

@ -335,6 +335,10 @@ void reactingOneDim::solveEnergy()
hEqn += fvc::div(phiQr); hEqn += fvc::div(phiQr);
} }
/* NOTE: The moving mesh option is only correct for reaction such as
Solid -> Gas, thus the ddt term is compesated exaclty by chemistrySh and
the mesh flux is not necessary.
if (regionMesh().moving()) if (regionMesh().moving())
{ {
surfaceScalarField phihMesh surfaceScalarField phihMesh
@ -344,7 +348,7 @@ void reactingOneDim::solveEnergy()
hEqn -= fvc::div(phihMesh); hEqn -= fvc::div(phihMesh);
} }
*/
hEqn.relax(); hEqn.relax();
hEqn.solve(); hEqn.solve();
} }

View File

@ -34,7 +34,7 @@ namespace Foam
} }
void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const void Foam::ensightSurfaceReader::skip(const label n, Istream& is) const
{ {
label i = 0; label i = 0;
token t; token t;
@ -59,6 +59,40 @@ void Foam::ensightSurfaceReader::skip(const label n, IFstream& is) const
} }
void Foam::ensightSurfaceReader::readLine(IFstream& is, string& buffer) const
{
buffer = "";
while (is.good() && buffer == "")
{
is.getLine(buffer);
}
}
void Foam::ensightSurfaceReader::debugSection
(
const word& expected,
IFstream& is
) const
{
string actual = "";
readLine(is, actual);
if (expected != actual)
{
FatalIOErrorInFunction(is)
<< "Expected section header '" << expected
<< "' but read the word '" << actual << "'"
<< exit(FatalIOError);
}
if (debug)
{
Info<< "Read section header: " << expected << endl;
}
}
void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const
{ {
// Binary flag string if applicable // Binary flag string if applicable
@ -101,29 +135,6 @@ void Foam::ensightSurfaceReader::readGeometryHeader(ensightReadFile& is) const
} }
void Foam::ensightSurfaceReader::debugSection
(
const word& expected,
IFstream& is
) const
{
word actual(is);
if (expected != actual)
{
FatalIOErrorInFunction(is)
<< "Expected section header '" << expected
<< "' but read the word '" << actual << "'"
<< exit(FatalIOError);
}
if (debug)
{
Info<< "Read section header: " << expected << endl;
}
}
void Foam::ensightSurfaceReader::readCase(IFstream& is) void Foam::ensightSurfaceReader::readCase(IFstream& is)
{ {
if (debug) if (debug)
@ -141,31 +152,54 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
string buffer; string buffer;
// Read the file // Read the file
debugSection("FORMAT", is); debugSection("FORMAT", is);
skip(3, is); // type: ensight gold readLine(is, buffer); // type: ensight gold
debugSection("GEOMETRY", is); debugSection("GEOMETRY", is);
readSkip(is, 2, meshFileName_); readLine(is, buffer);
readFromLine(2, buffer, meshFileName_); // model: 1 xxx.0000.mesh
debugSection("VARIABLE", is); debugSection("VARIABLE", is);
// Read the field description
DynamicList<word> fieldNames(10); DynamicList<word> fieldNames(10);
DynamicList<word> fieldFileNames(10); DynamicList<string> fieldFileNames(10);
word fieldName; word fieldName;
word fieldFileName; string fieldFileName;
while (is.good()) while (is.good())
{ {
word primitiveType(is); // scalar, vector readLine(is, buffer);
if (primitiveType == "TIME") if (buffer == "TIME")
{ {
break; break;
} }
readSkip(is, 3, fieldName); // p, U etc IStringStream iss(buffer);
// Read the field name, e.g. p U etc
readFromLine(4, iss, fieldName);
fieldNames.append(fieldName); fieldNames.append(fieldName);
is >> fieldFileName; // surfaceName.****.fieldName // Field file name may contain /'s e.g.
// surfaceName.****.fieldName
// This is not parser friendly - simply take remainder of buffer
label iPos = iss.stdStream().tellg();
fieldFileName = buffer(iPos, buffer.size() - iPos);
size_t p0 = fieldFileName.find_first_not_of(' ');
if (p0 == string::npos)
{
WarningInFunction
<< "Error reading field file name. "
<< "Current buffer: " << buffer
<< endl;
}
else
{
size_t p1 = fieldFileName.find_last_not_of(' ');
fieldFileName = fieldFileName.substr(p0, p1 - p0 + 1);
}
fieldFileNames.append(fieldFileName); fieldFileNames.append(fieldFileName);
} }
fieldNames_.transfer(fieldNames); fieldNames_.transfer(fieldNames);
@ -178,10 +212,14 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
} }
// Start reading time information // Start reading time information
skip(3, is); // time set: 1 readLine(is, buffer); // time set: 1
readSkip(is, 3, nTimeSteps_);
readSkip(is, 3, timeStartIndex_); readLine(is, buffer);
readSkip(is, 2, timeIncrement_); readFromLine(3, buffer, nTimeSteps_);
readLine(is, buffer);
readFromLine(3, buffer, timeStartIndex_);
readLine(is, buffer);
readFromLine(2, buffer, timeIncrement_);
if (debug) if (debug)
{ {
@ -191,7 +229,7 @@ void Foam::ensightSurfaceReader::readCase(IFstream& is)
} }
// Read the time values // Read the time values
skip(2, is); readLine(is, buffer); // time values:
timeValues_.setSize(nTimeSteps_); timeValues_.setSize(nTimeSteps_);
for (label i = 0; i < nTimeSteps_; i++) for (label i = 0; i < nTimeSteps_; i++)
{ {

View File

@ -69,7 +69,7 @@ protected:
List<word> fieldNames_; List<word> fieldNames_;
//- Field file names //- Field file names
List<word> fieldFileNames_; List<string> fieldFileNames_;
//- Number of time steps //- Number of time steps
label nTimeSteps_; label nTimeSteps_;
@ -91,21 +91,38 @@ protected:
// Protected Member Functions // Protected Member Functions
//- Helper function to skip forward n steps in stream
void skip(const label n, Istream& is) const;
//- Helper function to read an ascii line from file
void readLine(IFstream& is, string& buffer) const;
//- Read and check a section header //- Read and check a section header
void debugSection(const word& expected, IFstream& is) const; void debugSection(const word& expected, IFstream& is) const;
//- Read the case file
void readCase(IFstream& is);
//- Helper function to skip forward n steps in stream
void skip(const label n, IFstream& is) const;
//- Read (and throw away) geometry file header //- Read (and throw away) geometry file header
void readGeometryHeader(ensightReadFile& is) const; void readGeometryHeader(ensightReadFile& is) const;
//- Read the case file
void readCase(IFstream& is);
//- Helper function to return Type after skipping n tokens //- Helper function to return Type after skipping n tokens
template<class Type> template<class Type>
void readSkip(IFstream& is, const label nSkip, Type& value) const; void readFromLine
(
const label nSkip,
IStringStream& is,
Type& value
) const;
//- Helper function to return Type after skipping n tokens
template<class Type>
void readFromLine
(
const label nSkip,
const string& buffer,
Type& value
) const;
//- Helper function to return a field //- Helper function to return a field
template<class Type> template<class Type>

View File

@ -29,10 +29,10 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type> template<class Type>
void Foam::ensightSurfaceReader::readSkip void Foam::ensightSurfaceReader::readFromLine
( (
IFstream& is,
const label nSkip, const label nSkip,
IStringStream& is,
Type& value Type& value
) const ) const
{ {
@ -42,6 +42,20 @@ void Foam::ensightSurfaceReader::readSkip
} }
template<class Type>
void Foam::ensightSurfaceReader::readFromLine
(
const label nSkip,
const string& buffer,
Type& value
) const
{
IStringStream is(buffer);
readFromLine(nSkip, is, value);
}
template<class Type> template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
( (
@ -60,9 +74,19 @@ Foam::tmp<Foam::Field<Type> > Foam::ensightSurfaceReader::readField
fileName fieldFileName(fieldFileNames_[fieldIndex]); fileName fieldFileName(fieldFileNames_[fieldIndex]);
std::ostringstream oss; std::ostringstream oss;
oss << std::setfill('0') << std::setw(4) << fileIndex; label nMask = 0;
for (size_t chari = 0; chari < fieldFileName.size(); chari++)
{
if (fieldFileName[chari] == '*')
{
nMask++;
}
}
const std::string maskStr(nMask, '*');
oss << std::setfill('0') << std::setw(nMask) << fileIndex;
const word indexStr = oss.str(); const word indexStr = oss.str();
fieldFileName.replace("****", indexStr); fieldFileName.replace(maskStr, indexStr);
ensightReadFile is(baseDir_/fieldFileName, streamFormat_); ensightReadFile is(baseDir_/fieldFileName, streamFormat_);

View File

@ -23,7 +23,7 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
\page pageThermophsyicalModels Thermophsyical Models \page pageThermophsyicalModels Thermophysical Models
\section secSchemes Overview \section secSchemes Overview
The available thermophysical models are grouped into the following categories: The available thermophysical models are grouped into the following categories:

View File

@ -4,6 +4,8 @@ cd ${0%/*} || exit 1 # Run from this directory
# Source tutorial run functions # Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions . $WM_PROJECT_DIR/bin/tools/RunFunctions
application=$(getApplication)
# Create mesh # Create mesh
runApplication blockMesh runApplication blockMesh
@ -16,7 +18,7 @@ runApplication topoSet
# Create baffles and fields # Create baffles and fields
runApplication createBaffles -overwrite runApplication createBaffles -overwrite
runApplication $(getApplication) runApplication $application
#- RedistributePar to do decomposition #- RedistributePar to do decomposition
runParallel redistributePar -decompose -cellDist runParallel redistributePar -decompose -cellDist