Commit 6640acc6 authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

parents 012cbab0 be3d701b
......@@ -41,23 +41,23 @@ namespace Foam
void block::blockPoints()
{
// set local variables for mesh specification
label ni = blockDef_.n().x();
label nj = blockDef_.n().y();
label nk = blockDef_.n().z();
const label ni = blockDef_.n().x();
const label nj = blockDef_.n().y();
const label nk = blockDef_.n().z();
point start = blockDef_.points()[blockDef_.blockShape()[0]];
point xEnd = blockDef_.points()[blockDef_.blockShape()[1]];
point xyEnd = blockDef_.points()[blockDef_.blockShape()[2]];
point yEnd = blockDef_.points()[blockDef_.blockShape()[3]];
const point p000 = blockDef_.points()[blockDef_.blockShape()[0]];
const point p100 = blockDef_.points()[blockDef_.blockShape()[1]];
const point p110 = blockDef_.points()[blockDef_.blockShape()[2]];
const point p010 = blockDef_.points()[blockDef_.blockShape()[3]];
point zEnd = blockDef_.points()[blockDef_.blockShape()[4]];
point xzEnd = blockDef_.points()[blockDef_.blockShape()[5]];
point xyzEnd = blockDef_.points()[blockDef_.blockShape()[6]];
point yzEnd = blockDef_.points()[blockDef_.blockShape()[7]];
const point p001 = blockDef_.points()[blockDef_.blockShape()[4]];
const point p101 = blockDef_.points()[blockDef_.blockShape()[5]];
const point p111 = blockDef_.points()[blockDef_.blockShape()[6]];
const point p011 = blockDef_.points()[blockDef_.blockShape()[7]];
// set reference to the list of edge point and weighting factors
const List<List<point> >& edgePoints = blockDef_.blockEdgePoints();
const scalarListList& edgeWeights = blockDef_.blockEdgeWeights();
// list of edge point and weighting factors
const List<List<point> >& p = blockDef_.blockEdgePoints();
const scalarListList& w = blockDef_.blockEdgeWeights();
// generate vertices
......@@ -69,193 +69,138 @@ void block::blockPoints()
{
label vertexNo = vtxLabel(i, j, k);
vector edgex1 = start*(1.0 - edgeWeights[0][i])
+ xEnd*edgeWeights[0][i];
// points on edges
vector edgex1 = p000 + (p100 - p000)*w[0][i];
vector edgex2 = p010 + (p110 - p010)*w[1][i];
vector edgex3 = p011 + (p111 - p011)*w[2][i];
vector edgex4 = p001 + (p101 - p001)*w[3][i];
vector edgex2 = yEnd*(1.0 - edgeWeights[1][i])
+ xyEnd*edgeWeights[1][i];
vector edgey1 = p000 + (p010 - p000)*w[4][j];
vector edgey2 = p100 + (p110 - p100)*w[5][j];
vector edgey3 = p101 + (p111 - p101)*w[6][j];
vector edgey4 = p001 + (p011 - p001)*w[7][j];
vector edgex3 = yzEnd*(1.0 - edgeWeights[2][i])
+ xyzEnd*edgeWeights[2][i];
vector edgex4 = zEnd*(1.0 - edgeWeights[3][i])
+ xzEnd*edgeWeights[3][i];
vector edgey1 = start*(1.0 - edgeWeights[4][j])
+ yEnd*edgeWeights[4][j];
vector edgey2 = xEnd*(1.0 - edgeWeights[5][j])
+ xyEnd*edgeWeights[5][j];
vector edgey3 = xzEnd*(1.0 - edgeWeights[6][j])
+ xyzEnd*edgeWeights[6][j];
vector edgey4 = zEnd*(1.0 - edgeWeights[7][j])
+ yzEnd*edgeWeights[7][j];
vector edgez1 = start*(1.0 - edgeWeights[8][k])
+ zEnd*edgeWeights[8][k];
vector edgez2 = xEnd*(1.0 - edgeWeights[9][k])
+ xzEnd*edgeWeights[9][k];
vector edgez3 = xyEnd*(1.0 - edgeWeights[10][k])
+ xyzEnd*edgeWeights[10][k];
vector edgez4 = yEnd*(1.0 - edgeWeights[11][k])
+ yzEnd*edgeWeights[11][k];
vector edgez1 = p000 + (p001 - p000)*w[8][k];
vector edgez2 = p100 + (p101 - p100)*w[9][k];
vector edgez3 = p110 + (p111 - p110)*w[10][k];
vector edgez4 = p010 + (p011 - p010)*w[11][k];
// calculate the importance factors for all edges
// x - direction
scalar impx1 =
(
(1.0 - edgeWeights[0][i])
*(1.0 - edgeWeights[4][j])
*(1.0 - edgeWeights[8][k])
+ edgeWeights[0][i]
*(1.0 - edgeWeights[5][j])
*(1.0 - edgeWeights[9][k])
(1.0 - w[0][i])*(1.0 - w[4][j])*(1.0 - w[8][k])
+ w[0][i]*(1.0 - w[5][j])*(1.0 - w[9][k])
);
scalar impx2 =
(
(1.0 - edgeWeights[1][i])
*edgeWeights[4][j]
*(1.0 - edgeWeights[11][k])
+ edgeWeights[1][i]
*edgeWeights[5][j]
*(1.0 - edgeWeights[10][k])
(1.0 - w[1][i])*w[4][j]*(1.0 - w[11][k])
+ w[1][i]*w[5][j]*(1.0 - w[10][k])
);
scalar impx3 =
(
(1.0 - edgeWeights[2][i])
*edgeWeights[7][j]
*edgeWeights[11][k]
+ edgeWeights[2][i]
*edgeWeights[6][j]
*edgeWeights[10][k]
);
scalar impx3 =
(
(1.0 - w[2][i])*w[7][j]*w[11][k]
+ w[2][i]*w[6][j]*w[10][k]
);
scalar impx4 =
(
(1.0 - edgeWeights[3][i])
*(1.0 - edgeWeights[7][j])
*edgeWeights[8][k]
+ edgeWeights[3][i]
*(1.0 - edgeWeights[6][j])
*edgeWeights[9][k]
);
scalar impx4 =
(
(1.0 - w[3][i])*(1.0 - w[7][j])*w[8][k]
+ w[3][i]*(1.0 - w[6][j])*w[9][k]
);
scalar magImpx = impx1 + impx2 + impx3 + impx4;
impx1 /= magImpx;
impx2 /= magImpx;
impx3 /= magImpx;
impx4 /= magImpx;
// y - direction
scalar impy1 =
(
(1.0 - edgeWeights[4][j])
*(1.0 - edgeWeights[0][i])
*(1.0 - edgeWeights[8][k])
+ edgeWeights[4][j]
*(1.0 - edgeWeights[1][i])
*(1.0 - edgeWeights[11][k])
(1.0 - w[4][j])*(1.0 - w[0][i])*(1.0 - w[8][k])
+ w[4][j]*(1.0 - w[1][i])*(1.0 - w[11][k])
);
scalar impy2 =
(
(1.0 - edgeWeights[5][j])
*edgeWeights[0][i]
*(1.0 - edgeWeights[9][k])
+ edgeWeights[5][j]
*edgeWeights[1][i]
*(1.0 - edgeWeights[10][k])
(1.0 - w[5][j])*w[0][i]*(1.0 - w[9][k])
+ w[5][j]*w[1][i]*(1.0 - w[10][k])
);
scalar impy3 =
(
(1.0 - edgeWeights[6][j])
*edgeWeights[3][i]
*edgeWeights[9][k]
+ edgeWeights[6][j]
*edgeWeights[2][i]
*edgeWeights[10][k]
(1.0 - w[6][j])*w[3][i]*w[9][k]
+ w[6][j]*w[2][i]*w[10][k]
);
scalar impy4 =
(
(1.0 - edgeWeights[7][j])
*(1.0 - edgeWeights[3][i])
*edgeWeights[8][k]
+ edgeWeights[7][j]
*(1.0 - edgeWeights[2][i])
*edgeWeights[11][k]
(1.0 - w[7][j])*(1.0 - w[3][i])*w[8][k]
+ w[7][j]*(1.0 - w[2][i])*w[11][k]
);
scalar magImpy = impy1 + impy2 + impy3 + impy4;
impy1 /= magImpy;
impy2 /= magImpy;
impy3 /= magImpy;
impy4 /= magImpy;
// z - direction
scalar impz1 =
(
(1.0 - edgeWeights[8][k])
*(1.0 - edgeWeights[0][i])
*(1.0 - edgeWeights[4][j])
+ edgeWeights[8][k]
*(1.0 - edgeWeights[3][i])
*(1.0 - edgeWeights[7][j])
(1.0 - w[8][k])*(1.0 - w[0][i])*(1.0 - w[4][j])
+ w[8][k]*(1.0 - w[3][i])*(1.0 - w[7][j])
);
scalar impz2 =
(
(1.0 - edgeWeights[9][k])
*edgeWeights[0][i]
*(1.0 - edgeWeights[5][j])
+ edgeWeights[9][k]
*edgeWeights[3][i]
*(1.0 - edgeWeights[6][j])
(1.0 - w[9][k])*w[0][i]*(1.0 - w[5][j])
+ w[9][k]*w[3][i]*(1.0 - w[6][j])
);
scalar impz3 =
(
(1.0 - edgeWeights[10][k])
*edgeWeights[1][i]
*edgeWeights[5][j]
+ edgeWeights[10][k]
*edgeWeights[2][i]
*edgeWeights[6][j]
(1.0 - w[10][k])*w[1][i]*w[5][j]
+ w[10][k]*w[2][i]*w[6][j]
);
scalar impz4 =
(
(1.0 - edgeWeights[11][k])
*(1.0 - edgeWeights[1][i])
*edgeWeights[4][j]
+ edgeWeights[11][k]
*(1.0 - edgeWeights[2][i])
*edgeWeights[7][j]
(1.0 - w[11][k])*(1.0 - w[1][i])*w[4][j]
+ w[11][k]*(1.0 - w[2][i])*w[7][j]
);
scalar magImpz = impz1 + impz2 + impz3 + impz4;
impz1 /= magImpz;
impz2 /= magImpz;
impz3 /= magImpz;
impz4 /= magImpz;
// calculate the correction vectors
vector corx1 = impx1*(edgePoints[0][i] - edgex1);
vector corx2 = impx2*(edgePoints[1][i] - edgex2);
vector corx3 = impx3*(edgePoints[2][i] - edgex3);
vector corx4 = impx4*(edgePoints[3][i] - edgex4);
vector corx1 = impx1*(p[0][i] - edgex1);
vector corx2 = impx2*(p[1][i] - edgex2);
vector corx3 = impx3*(p[2][i] - edgex3);
vector corx4 = impx4*(p[3][i] - edgex4);
vector cory1 = impy1*(edgePoints[4][j] - edgey1);
vector cory2 = impy2*(edgePoints[5][j] - edgey2);
vector cory3 = impy3*(edgePoints[6][j] - edgey3);
vector cory4 = impy4*(edgePoints[7][j] - edgey4);
vector cory1 = impy1*(p[4][j] - edgey1);
vector cory2 = impy2*(p[5][j] - edgey2);
vector cory3 = impy3*(p[6][j] - edgey3);
vector cory4 = impy4*(p[7][j] - edgey4);
vector corz1 = impz1*(edgePoints[8][k] - edgez1);
vector corz2 = impz2*(edgePoints[9][k] - edgez2);
vector corz3 = impz3*(edgePoints[10][k] - edgez3);
vector corz4 = impz4*(edgePoints[11][k] - edgez4);
vector corz1 = impz1*(p[8][k] - edgez1);
vector corz2 = impz2*(p[9][k] - edgez2);
vector corz3 = impz3*(p[10][k] - edgez3);
vector corz4 = impz4*(p[11][k] - edgez4);
// multiply by the importance factor
// x - direction
edgex1 *= impx1;
edgex2 *= impx2;
......@@ -285,7 +230,6 @@ void block::blockPoints()
vertices_[vertexNo] += corx1 + corx2 + corx3 + corx4;
vertices_[vertexNo] += cory1 + cory2 + cory3 + cory4;
vertices_[vertexNo] += corz1 + corz2 + corz3 + corz4;
}
}
}
......
......@@ -79,12 +79,12 @@ int main(int argc, char *argv[])
Info<< " Reading volScalarField " << fieldName << endl;
volScalarField field(fieldHeader, mesh);
scalar area = sum(mesh.magSf().boundaryField()[patchi]);
scalar area = gSum(mesh.magSf().boundaryField()[patchi]);
scalar sumField = 0;
if (area > 0)
{
sumField = sum
sumField = gSum
(
mesh.magSf().boundaryField()[patchi]
* field.boundaryField()[patchi]
......
......@@ -75,14 +75,14 @@ int main(int argc, char *argv[])
}
// Give patch area
Info<< " Patch area = " << sum(mesh.Sf().boundaryField()[patchi]) << endl;
Info<< " Patch area = " << gSum(mesh.Sf().boundaryField()[patchi]) << endl;
if (fieldHeader.headerClassName() == "volScalarField")
{
Info<< " Reading volScalarField " << fieldName << endl;
volScalarField field(fieldHeader, mesh);
vector sumField = sum
vector sumField = gSum
(
mesh.Sf().boundaryField()[patchi]
* field.boundaryField()[patchi]
......@@ -97,7 +97,7 @@ int main(int argc, char *argv[])
Info<< " Reading surfaceScalarField " << fieldName << endl;
surfaceScalarField field(fieldHeader, mesh);
scalar sumField = sum(field.boundaryField()[patchi]);
scalar sumField = gSum(field.boundaryField()[patchi]);
Info<< " Integral of " << fieldName << " over patch "
<< patchName << '[' << patchi << ']' << " = "
......
......@@ -31,7 +31,7 @@
#
#------------------------------------------------------------------------------
export CINTSYSDIR=~/pub/CINT/cint
export CINTSYSDIR=~/pub/CINT/cint7
export PATH=$PATH:$CINTSYSDIR
export MANPATH=$MANPATH:$CINTSYSDIR/doc
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:$CINTSYSDIR
......
......@@ -33,15 +33,68 @@ void Foam::setRefCell
const volScalarField& field,
const dictionary& dict,
label& refCelli,
scalar& refValue
scalar& refValue,
bool forceReference
)
{
if (field.needReference())
if (field.needReference() || forceReference)
{
word refCellName = field.name() + "RefCell";
word refPointName = field.name() + "RefPoint";
word refValueName = field.name() + "RefValue";
refCelli = readLabel(dict.lookup(refCellName));
if (dict.found(refCellName))
{
if (Pstream::master())
{
refCelli = readLabel(dict.lookup(refCellName));
}
else
{
refCelli = -1;
}
}
else if (dict.found(refPointName))
{
point refPointi(dict.lookup(refPointName));
refCelli = field.mesh().findCell(refPointi);
label hasRef = (refCelli >= 0 ? 1 : 0);
label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());
if (sumHasRef != 1)
{
FatalErrorIn
(
"void Foam::setRefCell"
"("
" const volScalarField&,"
" const dictionary&,"
" label& scalar&,"
" bool"
")"
)
<< "Unable to set reference cell for field " << field.name()
<< nl << " Reference point " << refPointName
<< " found on multiple domains" << nl << abort(FatalError);
}
}
else
{
FatalErrorIn
(
"void Foam::setRefCell"
"("
" const volScalarField&,"
" const dictionary&,"
" label& scalar&,"
" bool"
")"
)
<< "Unable to set reference cell for field" << field.name() << nl
<< " Please supply either " << refCellName
<< " or " << refPointName << nl << abort(FatalError);
}
refValue = readScalar(dict.lookup(refValueName));
}
}
......
......@@ -52,7 +52,8 @@ void setRefCell
const volScalarField& field,
const dictionary& dict,
label& refCelli,
scalar& refValue
scalar& refValue,
bool forceReference = false
);
}
......
......@@ -234,6 +234,16 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
(
mpp.samplePatch()
);
if (patchID < 0)
{
FatalErrorIn
(
"void directMappedFixedValueFvPatchField<Type>::"
"updateCoeffs()"
)<< "Unable to find sample patch " << mpp.samplePatch()
<< " for patch " << this->patch().name() << nl
<< abort(FatalError);
}
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
const word& fieldName = this->dimensionedInternalField().name();
const fieldType& sendField =
......
......@@ -483,7 +483,7 @@ void Foam::fvMatrix<Type>::setReference
{
if (psi_.needReference() || forceReference)
{
if (Pstream::master())
if (cell >= 0)
{
source()[cell] += diag()[cell]*value;
diag()[cell] += diag()[cell];
......
......@@ -164,7 +164,7 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::hc() const
(
(((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd
+ a[0]
)*Tstd
)*Tstd + a[5]
);
}
......
......@@ -27,6 +27,7 @@ dynMixedSmagorinsky/dynMixedSmagorinsky.C
/*Smagorinsky2/Smagorinsky2.C*/
kOmegaSSTSAS/kOmegaSSTSAS.C
/* Wall functions */
wallFunctions=derivedFvPatchFields/wallFunctions
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "kOmegaSSTSAS.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace LESModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(kOmegaSSTSAS, 0);
addToRunTimeSelectionTable(LESModel, kOmegaSSTSAS, dictionary);
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<volScalarField> kOmegaSSTSAS::F1(const volScalarField& CDkOmega) const
{
volScalarField CDkOmegaPlus = max
(
CDkOmega,
dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
);
volScalarField arg1 = min
(
min
(
max
(
(scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*nu()/(sqr(y_)*omega_)
),
(4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
),
scalar(10)
);
return tanh(pow4(arg1));
}
tmp<volScalarField> kOmegaSSTSAS::F2() const