diff --git a/COPYING b/COPYING
index 81ee282d31c6594ffc4d3ceace3523bf5364f3c3..7e02780f0b5e88e02abfb88633829987f065d70d 100644
--- a/COPYING
+++ b/COPYING
@@ -3,7 +3,7 @@
OpenFOAM(R) is Copyright (C) 1991-2010 OpenCFD Ltd.
Contact: OpenCFD (enquiries@OpenCFD.co.uk)
- You may use, distribute and copy the OpenCFD CFD Toolbox under the terms
+ You may use, distribute and copy the OpenFOAM CFD Toolbox under the terms
of GNU General Public License version 3, which is displayed below, or
(at your option) any later version.
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
index 4150cd7503f41d781ce066ae50ee9acfd915e00c..5460944288567d39edaf3f1c3d94ab7051b6bae5 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
index c41bc9b6c7b93cab65751d97b441f519e4b3d05f..401e1203d2617de63e107506389ac1beaea3f4ec 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index 690b53760d75b22618fc9cd0c1358c015d867854..122280cfac48ecd0264c6463fa711928868d7c52 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -44,9 +44,14 @@
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
- dimensionedScalar pMin
+ dimensionedScalar rhoMax
(
- mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
+ mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
+ );
+
+ dimensionedScalar rhoMin
+ (
+ mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
);
Info<< "Creating turbulence model\n" << endl;
diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
index e299d99f83c9b1e0c0ab95b35d015300f4f18d7b..57395e977ed3520730cb0efaea428e2bae2b5da9 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
@@ -5,8 +5,8 @@
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
- fvc::div(phi/fvc::interpolate(rho), p, "div(U,p)")
- - p*fvc::div(phi/fvc::interpolate(rho))
+ fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
+ - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
);
hEqn.relax();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 38922c99c9c5700518251a347c167b6848f3de62..329ebbd27fab99c901e1f727669ed2ce21ed5ebb 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,4 +1,7 @@
rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
@@ -82,15 +85,9 @@ else
// Explicitly relax pressure for momentum corrector
p.relax();
-rho = thermo.rho();
-rho.relax();
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
-
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
-bound(p, pMin);
-
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
@@ -98,3 +95,9 @@ if (closedVolume)
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
+
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/compressible/rhoSimplecFoam/UEqn.H b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
index c41bc9b6c7b93cab65751d97b441f519e4b3d05f..401e1203d2617de63e107506389ac1beaea3f4ec 100644
--- a/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimplecFoam/UEqn.H
@@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index a813a2d8b7aab17f8a45dfc9295e9ee6f161a683..bcf99e166809c915d3c5763dc6abd0366d761507 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U)
);
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
index aa476b4dc808f46619dd5c94fd9d8defaa10c55d..0170d6e38446f8831bb4c170c7190a1742708e10 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
index dcec73b26841787fb72ab1a6395e530c826a4dcf..68d70ae7187555261c2768c20fbb586f74fdf3da 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
@@ -1,7 +1,6 @@
regionProperties/regionProperties.C
derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
-derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
fluid/compressibleCourantNo.C
solid/solidRegionDiffNo.C
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
deleted file mode 100644
index 42a73343ccae04725b61e3b6b39a59adec180c63..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
+++ /dev/null
@@ -1,381 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2009 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 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 "solidWallMixedTemperatureCoupledFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "directMappedPatchBase.H"
-#include "regionProperties.H"
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-bool Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::interfaceOwner
-(
- const polyMesh& nbrRegion
-) const
-{
- const fvMesh& myRegion = patch().boundaryMesh().mesh();
-
- const regionProperties& props =
- myRegion.objectRegistry::parent().lookupObject<regionProperties>
- (
- "regionProperties"
- );
-
- label myIndex = findIndex(props.fluidRegionNames(), myRegion.name());
- if (myIndex == -1)
- {
- label i = findIndex(props.solidRegionNames(), myRegion.name());
-
- if (i == -1)
- {
- FatalErrorIn
- (
- "solidWallMixedTemperatureCoupledFvPatchScalarField"
- "::interfaceOwner(const polyMesh&) const"
- ) << "Cannot find region " << myRegion.name()
- << " neither in fluids " << props.fluidRegionNames()
- << " nor in solids " << props.solidRegionNames()
- << exit(FatalError);
- }
- myIndex = props.fluidRegionNames().size() + i;
- }
- label nbrIndex = findIndex(props.fluidRegionNames(), nbrRegion.name());
- if (nbrIndex == -1)
- {
- label i = findIndex(props.solidRegionNames(), nbrRegion.name());
-
- if (i == -1)
- {
- FatalErrorIn("coupleManager::interfaceOwner(const polyMesh&) const")
- << "Cannot find region " << nbrRegion.name()
- << " neither in fluids " << props.fluidRegionNames()
- << " nor in solids " << props.solidRegionNames()
- << exit(FatalError);
- }
- nbrIndex = props.fluidRegionNames().size() + i;
- }
-
- return myIndex < nbrIndex;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
- const fvPatch& p,
- const DimensionedField<scalar, volMesh>& iF
-)
-:
- mixedFvPatchScalarField(p, iF),
- neighbourFieldName_("undefined-neighbourFieldName"),
- KName_("undefined-K")
-{
- this->refValue() = 0.0;
- this->refGrad() = 0.0;
- this->valueFraction() = 1.0;
- this->fixesValue_ = true;
-}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
- const solidWallMixedTemperatureCoupledFvPatchScalarField& ptf,
- const fvPatch& p,
- const DimensionedField<scalar, volMesh>& iF,
- const fvPatchFieldMapper& mapper
-)
-:
- mixedFvPatchScalarField(ptf, p, iF, mapper),
- neighbourFieldName_(ptf.neighbourFieldName_),
- KName_(ptf.KName_),
- fixesValue_(ptf.fixesValue_)
-{}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
- const fvPatch& p,
- const DimensionedField<scalar, volMesh>& iF,
- const dictionary& dict
-)
-:
- mixedFvPatchScalarField(p, iF),
- neighbourFieldName_(dict.lookup("neighbourFieldName")),
- KName_(dict.lookup("K"))
-{
- if (!isA<directMappedPatchBase>(this->patch().patch()))
- {
- FatalErrorIn
- (
- "solidWallMixedTemperatureCoupledFvPatchScalarField::"
- "solidWallMixedTemperatureCoupledFvPatchScalarField\n"
- "(\n"
- " const fvPatch& p,\n"
- " const DimensionedField<scalar, volMesh>& iF,\n"
- " const dictionary& dict\n"
- ")\n"
- ) << "\n patch type '" << p.type()
- << "' not type '" << directMappedPatchBase::typeName << "'"
- << "\n for patch " << p.name()
- << " of field " << dimensionedInternalField().name()
- << " in file " << dimensionedInternalField().objectPath()
- << exit(FatalError);
- }
-
- fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
-
- if (dict.found("refValue"))
- {
- // Full restart
- refValue() = scalarField("refValue", dict, p.size());
- refGrad() = scalarField("refGradient", dict, p.size());
- valueFraction() = scalarField("valueFraction", dict, p.size());
- fixesValue_ = readBool(dict.lookup("fixesValue"));
- }
- else
- {
- // Start from user entered data. Assume fixedValue.
- refValue() = *this;
- refGrad() = 0.0;
- valueFraction() = 1.0;
- fixesValue_ = true;
- }
-}
-
-
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::
-solidWallMixedTemperatureCoupledFvPatchScalarField
-(
- const solidWallMixedTemperatureCoupledFvPatchScalarField& wtcsf,
- const DimensionedField<scalar, volMesh>& iF
-)
-:
- mixedFvPatchScalarField(wtcsf, iF),
- neighbourFieldName_(wtcsf.neighbourFieldName_),
- KName_(wtcsf.KName_),
- fixesValue_(wtcsf.fixesValue_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-const Foam::fvPatchScalarField&
-Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::K() const
-{
- return this->patch().lookupPatchField<volScalarField, scalar>(KName_);
-}
-
-
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::updateCoeffs()
-{
- if (updated())
- {
- return;
- }
-
- // Get the coupling information from the directMappedPatchBase
- const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
- (
- patch().patch()
- );
- const polyMesh& nbrMesh = mpp.sampleMesh();
- // Force recalculation of mapping and schedule
- const mapDistribute& distMap = mpp.map();
- (void)distMap.schedule();
-
- tmp<scalarField> intFld = patchInternalField();
-
- if (interfaceOwner(nbrMesh))
- {
- // Note: other side information could be cached - it only needs
- // to be updated the first time round the iteration (i.e. when
- // switching regions) but unfortunately we don't have this information.
-
- const fvPatch& nbrPatch = refCast<const fvMesh>
- (
- nbrMesh
- ).boundary()[mpp.samplePolyPatch().index()];
-
-
- // Calculate the temperature by harmonic averaging
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- const solidWallMixedTemperatureCoupledFvPatchScalarField& nbrField =
- refCast<const solidWallMixedTemperatureCoupledFvPatchScalarField>
- (
- nbrPatch.lookupPatchField<volScalarField, scalar>
- (
- neighbourFieldName_
- )
- );
-
- // Swap to obtain full local values of neighbour internal field
- scalarField nbrIntFld = nbrField.patchInternalField();
- mapDistribute::distribute
- (
- Pstream::defaultCommsType,
- distMap.schedule(),
- distMap.constructSize(),
- distMap.subMap(), // what to send
- distMap.constructMap(), // what to receive
- nbrIntFld
- );
-
- // Swap to obtain full local values of neighbour K*delta
- scalarField nbrKDelta = nbrField.K()*nbrPatch.deltaCoeffs();
- mapDistribute::distribute
- (
- Pstream::defaultCommsType,
- distMap.schedule(),
- distMap.constructSize(),
- distMap.subMap(), // what to send
- distMap.constructMap(), // what to receive
- nbrKDelta
- );
-
-
- tmp<scalarField> myKDelta = K()*patch().deltaCoeffs();
-
- // Calculate common wall temperature. Reuse *this to store common value.
- scalarField Twall
- (
- (myKDelta()*intFld() + nbrKDelta*nbrIntFld)
- / (myKDelta() + nbrKDelta)
- );
- // Assign to me
- fvPatchScalarField::operator=(Twall);
- // Distribute back and assign to neighbour
- mapDistribute::distribute
- (
- Pstream::defaultCommsType,
- distMap.schedule(),
- nbrField.size(),
- distMap.constructMap(), // reverse : what to send
- distMap.subMap(),
- Twall
- );
- const_cast<solidWallMixedTemperatureCoupledFvPatchScalarField&>
- (
- nbrField
- ).fvPatchScalarField::operator=(Twall);
- }
-
-
- // Switch between fixed value (of harmonic avg) or gradient
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- label nFixed = 0;
-
- // Like snGrad but bypass switching on refValue/refGrad.
- tmp<scalarField> normalGradient = (*this-intFld())*patch().deltaCoeffs();
-
- if (debug)
- {
- scalar Q = gSum(K()*patch().magSf()*normalGradient());
-
- Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
- << "updateCoeffs() :"
- << " patch:" << patch().name()
- << " heatFlux:" << Q
- << " walltemperature "
- << " min:" << gMin(*this)
- << " max:" << gMax(*this)
- << " avg:" << gAverage(*this)
- << endl;
- }
-
- forAll(*this, i)
- {
- // if outgoing flux use fixed value.
- if (normalGradient()[i] < 0.0)
- {
- this->refValue()[i] = operator[](i);
- this->refGrad()[i] = 0.0; // not used by me
- this->valueFraction()[i] = 1.0;
- nFixed++;
- }
- else
- {
- // Fixed gradient. Make sure to have valid refValue (even though
- // I am not using it - other boundary conditions might)
- this->refValue()[i] = operator[](i);
- this->refGrad()[i] = normalGradient()[i];
- this->valueFraction()[i] = 0.0;
- }
- }
-
- reduce(nFixed, sumOp<label>());
-
- fixesValue_ = (nFixed > 0);
-
- if (debug)
- {
- label nTotSize = returnReduce(this->size(), sumOp<label>());
-
- Info<< "solidWallMixedTemperatureCoupledFvPatchScalarField::"
- << "updateCoeffs() :"
- << " patch:" << patch().name()
- << " out of:" << nTotSize
- << " fixedBC:" << nFixed
- << " gradient:" << nTotSize-nFixed << endl;
- }
-
- mixedFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::solidWallMixedTemperatureCoupledFvPatchScalarField::write
-(
- Ostream& os
-) const
-{
- mixedFvPatchScalarField::write(os);
- os.writeKeyword("neighbourFieldName")<< neighbourFieldName_
- << token::END_STATEMENT << nl;
- os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
- os.writeKeyword("fixesValue") << fixesValue_ << token::END_STATEMENT << nl;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-makePatchTypeField
-(
- fvPatchScalarField,
- solidWallMixedTemperatureCoupledFvPatchScalarField
-);
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
deleted file mode 100644
index f4766e5dc337baf2654a1386323806d8441d906a..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
+++ /dev/null
@@ -1,192 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2009 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 3 of the License, or
- (at your option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
-
-Class
- solidWallMixedTemperatureCoupledFvPatchScalarField
-
-Description
- Mixed boundary condition for temperature, to be used by the
- conjugate heat transfer solver.
- If my temperature is T1, neighbour is T2:
-
- T1 > T2: my side becomes fixedValue T2 bc, other side becomes fixedGradient.
-
-
- Example usage:
- myInterfacePatchName
- {
- type solidWallMixedTemperatureCoupled;
- neighbourFieldName T;
- K K;
- value uniform 300;
- }
-
- Needs to be on underlying directMapped(Wall)FvPatch.
-
- Note: runs in parallel with arbitrary decomposition. Uses directMapped
- functionality to calculate exchange.
-
- Note: lags interface data so both sides use same data.
- - problem: schedule to calculate average would interfere
- with standard processor swaps.
- - so: updateCoeffs sets both to same Twall. Only need to do
- this for last outer iteration but don't have access to this.
-
-SourceFiles
- solidWallMixedTemperatureCoupledFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef solidWallMixedTemperatureCoupledFvPatchScalarField_H
-#define solidWallMixedTemperatureCoupledFvPatchScalarField_H
-
-#include "fvPatchFields.H"
-#include "mixedFvPatchFields.H"
-#include "fvPatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
- Class solidWallMixedTemperatureCoupledFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class solidWallMixedTemperatureCoupledFvPatchScalarField
-:
- public mixedFvPatchScalarField
-{
- // Private data
-
- //- Name of field on the neighbour region
- const word neighbourFieldName_;
-
- //- Name of thermal conductivity field
- const word KName_;
-
- bool fixesValue_;
-
-
- // Private Member Functions
-
- //- Am I or neighbour owner of interface
- bool interfaceOwner(const polyMesh& nbrRegion) const;
-
-
-public:
-
- //- Runtime type information
- TypeName("solidWallMixedTemperatureCoupled");
-
-
- // Constructors
-
- //- Construct from patch and internal field
- solidWallMixedTemperatureCoupledFvPatchScalarField
- (
- const fvPatch&,
- const DimensionedField<scalar, volMesh>&
- );
-
- //- Construct from patch, internal field and dictionary
- solidWallMixedTemperatureCoupledFvPatchScalarField
- (
- const fvPatch&,
- const DimensionedField<scalar, volMesh>&,
- const dictionary&
- );
-
- //- Construct by mapping given
- // solidWallMixedTemperatureCoupledFvPatchScalarField onto a new patch
- solidWallMixedTemperatureCoupledFvPatchScalarField
- (
- const solidWallMixedTemperatureCoupledFvPatchScalarField&,
- const fvPatch&,
- const DimensionedField<scalar, volMesh>&,
- const fvPatchFieldMapper&
- );
-
- //- Construct and return a clone
- virtual tmp<fvPatchScalarField> clone() const
- {
- return tmp<fvPatchScalarField>
- (
- new solidWallMixedTemperatureCoupledFvPatchScalarField(*this)
- );
- }
-
- //- Construct as copy setting internal field reference
- solidWallMixedTemperatureCoupledFvPatchScalarField
- (
- const solidWallMixedTemperatureCoupledFvPatchScalarField&,
- const DimensionedField<scalar, volMesh>&
- );
-
- //- Construct and return a clone setting internal field reference
- virtual tmp<fvPatchScalarField> clone
- (
- const DimensionedField<scalar, volMesh>& iF
- ) const
- {
- return tmp<fvPatchScalarField>
- (
- new solidWallMixedTemperatureCoupledFvPatchScalarField
- (
- *this,
- iF
- )
- );
- }
-
-
- // Member functions
-
- //- Get corresponding K field
- const fvPatchScalarField& K() const;
-
- //- Return true if this patch field fixes a value.
- // Needed to check if a level has to be specified while solving
- // Poissons equations.
- virtual bool fixesValue() const
- {
- return fixesValue_;
- }
-
- //- Update the coefficients associated with the patch field
- virtual void updateCoeffs();
-
- //- Write
- virtual void write(Ostream&) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
index 8f1f654d2b6d133a3d7adabcf77dfcf83955f59c..b5db5078f1ab81134f8efe86fbc6d67a8fab3984 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
@@ -1,18 +1,16 @@
scalar CoNum = -GREAT;
- if (fluidRegions.size())
+
+ forAll(fluidRegions, regionI)
{
- forAll(fluidRegions, regionI)
- {
- CoNum = max
+ CoNum = max
+ (
+ compressibleCourantNo
(
- compressibleCourantNo
- (
- fluidRegions[regionI],
- runTime,
- rhoFluid[regionI],
- phiFluid[regionI]
- ),
- CoNum
- );
- }
+ fluidRegions[regionI],
+ runTime,
+ rhoFluid[regionI],
+ phiFluid[regionI]
+ ),
+ CoNum
+ );
}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
index 20ca225676f734a10382746e9533badfd6b12d78..c2e5668a684d16cf09256550d3da8e67a271c587 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
@@ -1,18 +1,16 @@
scalar DiNum = -GREAT;
- if (solidRegions.size())
+
+ forAll(solidRegions, regionI)
{
- forAll(solidRegions, regionI)
- {
- DiNum = max
+ DiNum = max
+ (
+ solidRegionDiffNo
(
- solidRegionDiffNo
- (
- solidRegions[regionI],
- runTime,
- rhosCps[regionI],
- Ks[regionI]
- ),
- DiNum
- );
- }
+ solidRegions[regionI],
+ runTime,
+ rhosCps[regionI],
+ Ks[regionI]
+ ),
+ DiNum
+ );
}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
index bf48d472925603c82b5f45b4f3e2f37b17f240db..e2b1d536d29814d07781d20e86761c6fd409fac4 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
@@ -37,7 +37,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
- fixedValueFvPatchScalarField(p, iF),
+ fixedGradientFvPatchScalarField(p, iF),
q_(p.size(), 0.0),
KName_("undefined-K")
{}
@@ -52,7 +52,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
- fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+ fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
q_(ptf.q_, mapper),
KName_(ptf.KName_)
{}
@@ -66,7 +66,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const dictionary& dict
)
:
- fixedValueFvPatchScalarField(p, iF, dict),
+ fixedGradientFvPatchScalarField(p, iF, dict),
q_("q", dict, p.size()),
KName_(dict.lookup("K"))
{}
@@ -78,7 +78,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf
)
:
- fixedValueFvPatchScalarField(tppsf),
+ fixedGradientFvPatchScalarField(tppsf),
q_(tppsf.q_),
KName_(tppsf.KName_)
{}
@@ -91,7 +91,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
- fixedValueFvPatchScalarField(tppsf, iF),
+ fixedGradientFvPatchScalarField(tppsf, iF),
q_(tppsf.q_),
KName_(tppsf.KName_)
{}
@@ -104,7 +104,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::autoMap
const fvPatchFieldMapper& m
)
{
- fixedValueFvPatchScalarField::autoMap(m);
+ fixedGradientFvPatchScalarField::autoMap(m);
q_.autoMap(m);
}
@@ -115,7 +115,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
const labelList& addr
)
{
- fixedValueFvPatchScalarField::rmap(ptf, addr);
+ fixedGradientFvPatchScalarField::rmap(ptf, addr);
const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf =
refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf);
@@ -131,14 +131,14 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
return;
}
- const scalarField& Kw =
- patch().lookupPatchField<volScalarField, scalar>(KName_);
-
- const fvPatchScalarField& Tw = *this;
+ const scalarField& Kw = patch().lookupPatchField<volScalarField, scalar>
+ (
+ KName_
+ );
- operator==(q_/(patch().deltaCoeffs()*Kw) + Tw.patchInternalField());
+ gradient() = q_/Kw;
- fixedValueFvPatchScalarField::updateCoeffs();
+ fixedGradientFvPatchScalarField::updateCoeffs();
}
@@ -147,9 +147,10 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
Ostream& os
) const
{
- fixedValueFvPatchScalarField::write(os);
+ fixedGradientFvPatchScalarField::write(os);
q_.writeEntry("q", os);
os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+ this->writeEntry("value", os);
}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
index 7558c4df7eb5d088900fb37634fc3906ea8278b1..ffbcf52ed5c01bb77010a1bdce5e5b062895953a 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
@@ -45,7 +45,7 @@ SourceFiles
#ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
#define solidWallHeatFluxTemperatureFvPatchScalarField_H
-#include "fixedValueFvPatchFields.H"
+#include "fixedGradientFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,7 +58,7 @@ namespace Foam
class solidWallHeatFluxTemperatureFvPatchScalarField
:
- public fixedValueFvPatchScalarField
+ public fixedGradientFvPatchScalarField
{
// Private data
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
index 1cbafe368aa4b6d2c4f9dddaa63fbe8db0ca1175..7f317cc2e6bb727b81534d45bdc91bbf79df96bc 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
@@ -2,7 +2,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turb.divDevRhoReff(U)
);
diff --git a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
index c1925c7708c8fb821515092efcc420b22f816ad1..688120809a1dfb9821d74ff7d7274e4fb50183d0 100644
--- a/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
+++ b/applications/solvers/incompressible/porousSimpleFoam/UEqn.H
@@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U)
);
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index e16d90bbbd1adbd1be557973de6b8107776f9d97..50ff6918b089b2b67e97828b3a06cb5589ca5981 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -291,6 +291,22 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
}
}
+ if (allGeometry)
+ {
+ cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
+ if (mesh.checkConcaveCells(true, &cells))
+ {
+ noFailedChecks++;
+
+ label nCells = returnReduce(cells.size(), sumOp<label>());
+
+ Info<< " <<Writing " << nCells
+ << " concave cells to set " << cells.name() << endl;
+ cells.instance() = mesh.pointsInstance();
+ cells.write();
+ }
+ }
+
return noFailedChecks;
}
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index 88ca441e9874c7ce3e8a753ca33ef96144c1b163..7da3c2226908e908245cd1dcd1a2b6bcd986bb81 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -39,17 +39,29 @@ Description
- mesh with cells put into cellZones (-makeCellZones)
Note:
+ - cellZonesOnly does not do a walk and uses the cellZones only. Use
+ this if you don't mind having disconnected domains in a single region.
+ This option requires all cells to be in one (and one only) cellZone.
+
+ - cellZonesFileOnly behaves like -cellZonesOnly but reads the cellZones
+ from the specified file. This allows one to explicitly specify the region
+ distribution and still have multiple cellZones per region.
+
+ - useCellZonesOnly does not do a walk and uses the cellZones only. Use
+ this if you don't mind having disconnected domains in a single region.
+ This option requires all cells to be in one (and one only) cellZone.
+
+
- Should work in parallel.
cellZones can differ on either side of processor boundaries in which case
the faces get moved from processor patch to directMapped patch. Not
- ery well tested.
+ very well tested.
+
- If a cell zone gets split into more than one region it can detect
the largest matching region (-sloppyCellZones). This will accept any
region that covers more than 50% of the zone. It has to be a subset
so cannot have any cells in any other zone.
- - useCellZonesOnly does not do a walk and uses the cellZones only. Use
- this if you don't mind having disconnected domains in a single region.
- This option requires all cells to be in one (and one only) cellZone.
+
- writes maps like decomposePar back to original mesh:
- pointRegionAddressing : for every point in this region the point in
the original mesh
@@ -1137,63 +1149,6 @@ EdgeMap<label> addRegionPatches
}
-//// Checks if regionI in cellRegion is subset of existing cellZone. Returns -1
-//// if no zone found, zone otherwise
-//label findCorrespondingSubZone
-//(
-// const cellZoneMesh& cellZones,
-// const labelList& existingZoneID,
-// const labelList& cellRegion,
-// const label regionI
-//)
-//{
-// // Zone corresponding to region. No corresponding zone.
-// label zoneI = labelMax;
-//
-// labelList regionCells = findIndices(cellRegion, regionI);
-//
-// if (regionCells.empty())
-// {
-// // My local portion is empty. Maps to any empty cellZone. Mark with
-// // special value which can get overwritten by other processors.
-// zoneI = -1;
-// }
-// else
-// {
-// // Get zone for first element.
-// zoneI = existingZoneID[regionCells[0]];
-//
-// if (zoneI == -1)
-// {
-// zoneI = labelMax;
-// }
-// else
-// {
-// // 1. All regionCells in zoneI?
-// forAll(regionCells, i)
-// {
-// if (existingZoneID[regionCells[i]] != zoneI)
-// {
-// zoneI = labelMax;
-// break;
-// }
-// }
-// }
-// }
-//
-// // Determine same zone over all processors.
-// reduce(zoneI, maxOp<label>());
-//
-// if (zoneI == labelMax)
-// {
-// // Cells in region that are not in zoneI
-// zoneI = -1;
-// }
-//
-// return zoneI;
-//}
-
-
// Find region that covers most of cell zone
label findCorrespondingRegion
(
@@ -1314,62 +1269,20 @@ label findCorrespondingRegion
//}
-// Main program:
-
-int main(int argc, char *argv[])
+// Get zone per cell
+// - non-unique zoning
+// - coupled zones
+void getZoneID
+(
+ const polyMesh& mesh,
+ const cellZoneMesh& cellZones,
+ labelList& zoneID,
+ labelList& neiZoneID
+)
{
-# include "addOverwriteOption.H"
- argList::addBoolOption("cellZones");
- argList::addBoolOption("cellZonesOnly");
- argList::addOption("blockedFaces", "faceSet");
- argList::addBoolOption("makeCellZones");
- argList::addBoolOption("largestOnly");
- argList::addOption("insidePoint", "point");
- argList::addBoolOption("detectOnly");
- argList::addBoolOption("sloppyCellZones");
-
-# include "setRootCase.H"
-# include "createTime.H"
- runTime.functionObjects().off();
-# include "createMesh.H"
- const word oldInstance = mesh.pointsInstance();
-
- word blockedFacesName;
- if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
- {
- Info<< "Reading blocked internal faces from faceSet "
- << blockedFacesName << nl << endl;
- }
-
- const bool makeCellZones = args.optionFound("makeCellZones");
- const bool largestOnly = args.optionFound("largestOnly");
- const bool insidePoint = args.optionFound("insidePoint");
- const bool useCellZones = args.optionFound("cellZones");
- const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
- const bool overwrite = args.optionFound("overwrite");
- const bool detectOnly = args.optionFound("detectOnly");
- const bool sloppyCellZones = args.optionFound("sloppyCellZones");
-
- if (insidePoint && largestOnly)
- {
- FatalErrorIn(args.executable())
- << "You cannot specify both -largestOnly"
- << " (keep region with most cells)"
- << " and -insidePoint (keep region containing point)"
- << exit(FatalError);
- }
-
-
- const cellZoneMesh& cellZones = mesh.cellZones();
-
-
- // Collect zone per cell
- // ~~~~~~~~~~~~~~~~~~~~~
- // - non-unique zoning
- // - coupled zones
-
// Existing zoneID
- labelList zoneID(mesh.nCells(), -1);
+ zoneID.setSize(mesh.nCells());
+ zoneID = -1;
forAll(cellZones, zoneI)
{
@@ -1384,7 +1297,7 @@ int main(int argc, char *argv[])
}
else
{
- FatalErrorIn(args.executable())
+ FatalErrorIn("getZoneID(..)")
<< "Cell " << cellI << " with cell centre "
<< mesh.cellCentres()[cellI]
<< " is multiple zones. This is not allowed." << endl
@@ -1396,98 +1309,136 @@ int main(int argc, char *argv[])
}
// Neighbour zoneID.
- labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+ neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
forAll(neiZoneID, i)
{
neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
}
syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
+}
+
+
+void matchRegions
+(
+ const bool sloppyCellZones,
+ const polyMesh& mesh,
+ const label nCellRegions,
+ const labelList& cellRegion,
+
+ labelList& regionToZone,
+ wordList& regionNames,
+ labelList& zoneToRegion
+)
+{
+ const cellZoneMesh& cellZones = mesh.cellZones();
- // Determine connected regions
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ regionToZone.setSize(nCellRegions, -1);
+ regionNames.setSize(nCellRegions);
+ zoneToRegion.setSize(cellZones.size(), -1);
- // Mark additional faces that are blocked
- boolList blockedFace;
+ // Get current per cell zoneID
+ labelList zoneID(mesh.nCells(), -1);
+ labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+ getZoneID(mesh, cellZones, zoneID, neiZoneID);
- // Read from faceSet
- if (blockedFacesName.size())
+ // Sizes per cellzone
+ labelList zoneSizes(cellZones.size(), 0);
{
- faceSet blockedFaceSet(mesh, blockedFacesName);
- Info<< "Read " << returnReduce(blockedFaceSet.size(), sumOp<label>())
- << " blocked faces from set " << blockedFacesName << nl << endl;
+ List<wordList> zoneNames(Pstream::nProcs());
+ zoneNames[Pstream::myProcNo()] = cellZones.names();
+ Pstream::gatherList(zoneNames);
+ Pstream::scatterList(zoneNames);
- blockedFace.setSize(mesh.nFaces(), false);
+ forAll(zoneNames, procI)
+ {
+ if (zoneNames[procI] != zoneNames[0])
+ {
+ FatalErrorIn("matchRegions(..)")
+ << "cellZones not synchronised across processors." << endl
+ << "Master has cellZones " << zoneNames[0] << endl
+ << "Processor " << procI
+ << " has cellZones " << zoneNames[procI]
+ << exit(FatalError);
+ }
+ }
- forAllConstIter(faceSet, blockedFaceSet, iter)
+ forAll(cellZones, zoneI)
{
- blockedFace[iter.key()] = true;
+ zoneSizes[zoneI] = returnReduce
+ (
+ cellZones[zoneI].size(),
+ sumOp<label>()
+ );
}
}
- // Imply from differing cellZones
- if (useCellZones)
+
+ if (sloppyCellZones)
{
- blockedFace.setSize(mesh.nFaces(), false);
+ Info<< "Trying to match regions to existing cell zones;"
+ << " region can be subset of cell zone." << nl << endl;
- for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+ forAll(cellZones, zoneI)
{
- label own = mesh.faceOwner()[faceI];
- label nei = mesh.faceNeighbour()[faceI];
+ label regionI = findCorrespondingRegion
+ (
+ zoneID,
+ cellRegion,
+ nCellRegions,
+ zoneI,
+ label(0.5*zoneSizes[zoneI]) // minimum overlap
+ );
- if (zoneID[own] != zoneID[nei])
+ if (regionI != -1)
{
- blockedFace[faceI] = true;
+ Info<< "Sloppily matched region " << regionI
+ //<< " size " << regionSizes[regionI]
+ << " to zone " << zoneI << " size " << zoneSizes[zoneI]
+ << endl;
+ zoneToRegion[zoneI] = regionI;
+ regionToZone[regionI] = zoneI;
+ regionNames[regionI] = cellZones[zoneI].name();
}
}
+ }
+ else
+ {
+ Info<< "Trying to match regions to existing cell zones." << nl << endl;
- // Different cellZones on either side of processor patch.
- forAll(neiZoneID, i)
+ forAll(cellZones, zoneI)
{
- label faceI = i+mesh.nInternalFaces();
+ label regionI = findCorrespondingRegion
+ (
+ zoneID,
+ cellRegion,
+ nCellRegions,
+ zoneI,
+ 1 // minimum overlap
+ );
- if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
+ if (regionI != -1)
{
- blockedFace[faceI] = true;
+ zoneToRegion[zoneI] = regionI;
+ regionToZone[regionI] = zoneI;
+ regionNames[regionI] = cellZones[zoneI].name();
}
}
}
-
-
- // Determine per cell the region it belongs to
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // cellRegion is the labelList with the region per cell.
- labelList cellRegion;
- label nCellRegions = 0;
- if (useCellZonesOnly)
+ // Allocate region names for unmatched regions.
+ forAll(regionToZone, regionI)
{
- label unzonedCellI = findIndex(zoneID, -1);
- if (unzonedCellI != -1)
+ if (regionToZone[regionI] == -1)
{
- FatalErrorIn(args.executable())
- << "For the cellZonesOnly option all cells "
- << "have to be in a cellZone." << endl
- << "Cell " << unzonedCellI
- << " at" << mesh.cellCentres()[unzonedCellI]
- << " is not in a cellZone. There might be more unzoned cells."
- << exit(FatalError);
+ regionNames[regionI] = "domain" + Foam::name(regionI);
}
- cellRegion = zoneID;
- nCellRegions = gMax(cellRegion)+1;
- }
- else
- {
- // Do a topological walk to determine regions
- regionSplit regions(mesh, blockedFace);
- nCellRegions = regions.nRegions();
- cellRegion.transfer(regions);
}
-
- Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
+}
+void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
+{
// Write to manual decomposition option
{
labelIOList cellToRegion
@@ -1515,7 +1466,7 @@ int main(int argc, char *argv[])
IOobject
(
"cellToRegion",
- runTime.timeName(),
+ mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
@@ -1534,137 +1485,280 @@ int main(int argc, char *argv[])
Info<< "Writing region per cell as volScalarField to "
<< cellToRegion.objectPath() << nl << endl;
}
+}
- // Sizes per region
- // ~~~~~~~~~~~~~~~~
- labelList regionSizes(nCellRegions, 0);
+// Main program:
- forAll(cellRegion, cellI)
+int main(int argc, char *argv[])
+{
+# include "addOverwriteOption.H"
+ argList::addBoolOption("cellZones");
+ argList::addBoolOption("cellZonesOnly");
+ argList::addOption("cellZonesFileOnly", "cellZonesName");
+ argList::addOption("blockedFaces", "faceSet");
+ argList::addBoolOption("makeCellZones");
+ argList::addBoolOption("largestOnly");
+ argList::addOption("insidePoint", "point");
+ argList::addBoolOption("detectOnly");
+ argList::addBoolOption("sloppyCellZones");
+
+# include "setRootCase.H"
+# include "createTime.H"
+ runTime.functionObjects().off();
+# include "createMesh.H"
+ const word oldInstance = mesh.pointsInstance();
+
+ word blockedFacesName;
+ if (args.optionReadIfPresent("blockedFaces", blockedFacesName))
{
- regionSizes[cellRegion[cellI]]++;
+ Info<< "Reading blocked internal faces from faceSet "
+ << blockedFacesName << nl << endl;
}
- forAll(regionSizes, regionI)
+
+ const bool makeCellZones = args.optionFound("makeCellZones");
+ const bool largestOnly = args.optionFound("largestOnly");
+ const bool insidePoint = args.optionFound("insidePoint");
+ const bool useCellZones = args.optionFound("cellZones");
+ const bool useCellZonesOnly = args.optionFound("cellZonesOnly");
+ const bool useCellZonesFile = args.optionFound("cellZonesFileOnly");
+ const bool overwrite = args.optionFound("overwrite");
+ const bool detectOnly = args.optionFound("detectOnly");
+ const bool sloppyCellZones = args.optionFound("sloppyCellZones");
+ if
+ (
+ (useCellZonesOnly || useCellZonesFile)
+ && (
+ blockedFacesName != word::null
+ || useCellZones
+ )
+ )
{
- reduce(regionSizes[regionI], sumOp<label>());
+ FatalErrorIn(args.executable())
+ << "You cannot specify both -cellZonesOnly or -cellZonesFileOnly"
+ << " (which specify complete split)"
+ << " in combination with -blockedFaces or -cellZones"
+ << " (which imply a split based on topology)"
+ << exit(FatalError);
}
- Info<< "Region\tCells" << nl
- << "------\t-----" << endl;
- forAll(regionSizes, regionI)
+
+ if (insidePoint && largestOnly)
{
- Info<< regionI << '\t' << regionSizes[regionI] << nl;
+ FatalErrorIn(args.executable())
+ << "You cannot specify both -largestOnly"
+ << " (keep region with most cells)"
+ << " and -insidePoint (keep region containing point)"
+ << exit(FatalError);
}
- Info<< endl;
- // Sizes per cellzone
- // ~~~~~~~~~~~~~~~~~~
-
- labelList zoneSizes(cellZones.size(), 0);
- if (useCellZones || makeCellZones || sloppyCellZones)
- {
- List<wordList> zoneNames(Pstream::nProcs());
- zoneNames[Pstream::myProcNo()] = cellZones.names();
- Pstream::gatherList(zoneNames);
- Pstream::scatterList(zoneNames);
+ const cellZoneMesh& cellZones = mesh.cellZones();
- forAll(zoneNames, procI)
- {
- if (zoneNames[procI] != zoneNames[0])
- {
- FatalErrorIn(args.executable())
- << "cellZones not synchronised across processors." << endl
- << "Master has cellZones " << zoneNames[0] << endl
- << "Processor " << procI
- << " has cellZones " << zoneNames[procI]
- << exit(FatalError);
- }
- }
+ // Existing zoneID
+ labelList zoneID(mesh.nCells(), -1);
+ // Neighbour zoneID.
+ labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+ getZoneID(mesh, cellZones, zoneID, neiZoneID);
- forAll(cellZones, zoneI)
- {
- zoneSizes[zoneI] = returnReduce
- (
- cellZones[zoneI].size(),
- sumOp<label>()
- );
- }
- }
- // Whether region corresponds to a cellzone
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Determine per cell the region it belongs to
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // cellRegion is the labelList with the region per cell.
+ labelList cellRegion;
// Region per zone
- labelList regionToZone(nCellRegions, -1);
+ labelList regionToZone;
// Name of region
- wordList regionNames(nCellRegions);
+ wordList regionNames;
// Zone to region
- labelList zoneToRegion(cellZones.size(), -1);
+ labelList zoneToRegion;
- if (sloppyCellZones)
+ label nCellRegions = 0;
+ if (useCellZonesOnly)
{
- Info<< "Trying to match regions to existing cell zones;"
- << " region can be subset of cell zone." << nl << endl;
+ Info<< "Using current cellZones to split mesh into regions."
+ << " This requires all"
+ << " cells to be in one and only one cellZone." << nl << endl;
- forAll(cellZones, zoneI)
+ label unzonedCellI = findIndex(zoneID, -1);
+ if (unzonedCellI != -1)
{
- label regionI = findCorrespondingRegion
+ FatalErrorIn(args.executable())
+ << "For the cellZonesOnly option all cells "
+ << "have to be in a cellZone." << endl
+ << "Cell " << unzonedCellI
+ << " at" << mesh.cellCentres()[unzonedCellI]
+ << " is not in a cellZone. There might be more unzoned cells."
+ << exit(FatalError);
+ }
+ cellRegion = zoneID;
+ nCellRegions = gMax(cellRegion)+1;
+ regionToZone.setSize(nCellRegions);
+ regionNames.setSize(nCellRegions);
+ zoneToRegion.setSize(cellZones.size(), -1);
+ for (label regionI = 0; regionI < nCellRegions; regionI++)
+ {
+ regionToZone[regionI] = regionI;
+ zoneToRegion[regionI] = regionI;
+ regionNames[regionI] = cellZones[regionI].name();
+ }
+ }
+ else if (useCellZonesFile)
+ {
+ const word zoneFile = args.option("cellZonesFileOnly");
+ Info<< "Reading split from cellZones file " << zoneFile << endl
+ << "This requires all"
+ << " cells to be in one and only one cellZone." << nl << endl;
+
+ cellZoneMesh newCellZones
+ (
+ IOobject
(
- zoneID,
- cellRegion,
- nCellRegions,
- zoneI,
- label(0.5*zoneSizes[zoneI]) // minimum overlap
- );
+ zoneFile,
+ mesh.facesInstance(),
+ polyMesh::meshSubDir,
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh
+ );
- if (regionI != -1)
- {
- Info<< "Sloppily matched region " << regionI
- << " size " << regionSizes[regionI]
- << " to zone " << zoneI << " size " << zoneSizes[zoneI]
- << endl;
- zoneToRegion[zoneI] = regionI;
- regionToZone[regionI] = zoneI;
- regionNames[regionI] = cellZones[zoneI].name();
- }
+ labelList newZoneID(mesh.nCells(), -1);
+ labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+ getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
+
+ label unzonedCellI = findIndex(newZoneID, -1);
+ if (unzonedCellI != -1)
+ {
+ FatalErrorIn(args.executable())
+ << "For the cellZonesFileOnly option all cells "
+ << "have to be in a cellZone." << endl
+ << "Cell " << unzonedCellI
+ << " at" << mesh.cellCentres()[unzonedCellI]
+ << " is not in a cellZone. There might be more unzoned cells."
+ << exit(FatalError);
+ }
+ cellRegion = newZoneID;
+ nCellRegions = gMax(cellRegion)+1;
+ zoneToRegion.setSize(newCellZones.size(), -1);
+ regionToZone.setSize(nCellRegions);
+ regionNames.setSize(nCellRegions);
+ for (label regionI = 0; regionI < nCellRegions; regionI++)
+ {
+ regionToZone[regionI] = regionI;
+ zoneToRegion[regionI] = regionI;
+ regionNames[regionI] = newCellZones[regionI].name();
}
}
else
{
- Info<< "Trying to match regions to existing cell zones." << nl << endl;
+ // Determine connected regions
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
- forAll(cellZones, zoneI)
+ // Mark additional faces that are blocked
+ boolList blockedFace;
+
+ // Read from faceSet
+ if (blockedFacesName.size())
{
- label regionI = findCorrespondingRegion
- (
- zoneID,
- cellRegion,
- nCellRegions,
- zoneI,
- 1 // minimum overlap
- );
+ faceSet blockedFaceSet(mesh, blockedFacesName);
+ Info<< "Read "
+ << returnReduce(blockedFaceSet.size(), sumOp<label>())
+ << " blocked faces from set " << blockedFacesName << nl << endl;
- if (regionI != -1)
+ blockedFace.setSize(mesh.nFaces(), false);
+
+ forAllConstIter(faceSet, blockedFaceSet, iter)
{
- zoneToRegion[zoneI] = regionI;
- regionToZone[regionI] = zoneI;
- regionNames[regionI] = cellZones[zoneI].name();
+ blockedFace[iter.key()] = true;
}
}
- }
- // Allocate region names for unmatched regions.
- forAll(regionToZone, regionI)
- {
- if (regionToZone[regionI] == -1)
+
+ // Imply from differing cellZones
+ if (useCellZones)
{
- regionNames[regionI] = "domain" + Foam::name(regionI);
+ blockedFace.setSize(mesh.nFaces(), false);
+
+ for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+ {
+ label own = mesh.faceOwner()[faceI];
+ label nei = mesh.faceNeighbour()[faceI];
+
+ if (zoneID[own] != zoneID[nei])
+ {
+ blockedFace[faceI] = true;
+ }
+ }
+
+ // Different cellZones on either side of processor patch.
+ forAll(neiZoneID, i)
+ {
+ label faceI = i+mesh.nInternalFaces();
+
+ if (zoneID[mesh.faceOwner()[faceI]] != neiZoneID[i])
+ {
+ blockedFace[faceI] = true;
+ }
+ }
}
+
+ // Do a topological walk to determine regions
+ regionSplit regions(mesh, blockedFace);
+ nCellRegions = regions.nRegions();
+ cellRegion.transfer(regions);
+
+ // Make up region names. If possible match them to existing zones.
+ matchRegions
+ (
+ sloppyCellZones,
+ mesh,
+ nCellRegions,
+ cellRegion,
+
+ regionToZone,
+ regionNames,
+ zoneToRegion
+ );
+ }
+
+ Info<< endl << "Number of regions:" << nCellRegions << nl << endl;
+
+
+ // Write decomposition to file
+ writeCellToRegion(mesh, cellRegion);
+
+
+
+ // Sizes per region
+ // ~~~~~~~~~~~~~~~~
+
+ labelList regionSizes(nCellRegions, 0);
+
+ forAll(cellRegion, cellI)
+ {
+ regionSizes[cellRegion[cellI]]++;
+ }
+ forAll(regionSizes, regionI)
+ {
+ reduce(regionSizes[regionI], sumOp<label>());
}
+ Info<< "Region\tCells" << nl
+ << "------\t-----" << endl;
+
+ forAll(regionSizes, regionI)
+ {
+ Info<< regionI << '\t' << regionSizes[regionI] << nl;
+ }
+ Info<< endl;
+
+
// Print region to zone
Info<< "Region\tZone\tName" << nl
@@ -1676,16 +1770,6 @@ int main(int argc, char *argv[])
}
Info<< endl;
- //// Print zone to region
- //Info<< "Zone\tName\tRegion" << nl
- // << "----\t----\t------" << endl;
- //forAll(zoneToRegion, zoneI)
- //{
- // Info<< zoneI << '\t' << cellZones[zoneI].name() << '\t'
- // << zoneToRegion[zoneI] << nl;
- //}
- //Info<< endl;
-
// Since we're going to mess with patches make sure all non-processor ones
diff --git a/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C b/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
index d50fbcf0f8b15f52185f1af8c8c41a338d98a3b5..65bd73563dee93ded8b7aeda9f89c1e9527840a2 100644
--- a/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
+++ b/src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -37,6 +37,9 @@ ListType Foam::renumber
// Create copy
ListType newLst(lst.size());
+ // ensure consistent addressable size (eg, DynamicList)
+ newLst.setSize(lst.size());
+
forAll(lst, elemI)
{
if (lst[elemI] >= 0)
@@ -76,6 +79,9 @@ ListType Foam::reorder
// Create copy
ListType newLst(lst.size());
+ // ensure consistent addressable size (eg, DynamicList)
+ newLst.setSize(lst.size());
+
forAll(lst, elemI)
{
if (oldToNew[elemI] >= 0)
@@ -101,6 +107,9 @@ void Foam::inplaceReorder
// Create copy
ListType newLst(lst.size());
+ // ensure consistent addressable size (eg, DynamicList)
+ newLst.setSize(lst.size());
+
forAll(lst, elemI)
{
if (oldToNew[elemI] >= 0)
@@ -258,6 +267,9 @@ ListType Foam::subset
ListType newLst(lst.size());
+ // ensure consistent addressable size (eg, DynamicList)
+ newLst.setSize(lst.size());
+
label nElem = 0;
forAll(lst, elemI)
{
@@ -318,6 +330,9 @@ ListType Foam::subset
ListType newLst(lst.size());
+ // ensure consistent addressable size (eg, DynamicList)
+ newLst.setSize(lst.size());
+
label nElem = 0;
forAll(lst, elemI)
{
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index ab7193a682998d02a45d6a6af7b6faeee3d328aa..0d47c7db9d4aea9f7811e58bbc5887ffc1d946cd 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -70,6 +70,7 @@ namespace Foam
class mapPolyMesh;
class globalIndex;
+class PstreamBuffers;
/*---------------------------------------------------------------------------*\
Class mapDistribute Declaration
@@ -318,6 +319,13 @@ public:
}
}
+ //- Do all sends using PstreamBuffers
+ template<class T>
+ void send(PstreamBuffers&, const List<T>&) const;
+ //- Do all receives using PstreamBuffers
+ template<class T>
+ void receive(PstreamBuffers&, List<T>&) const;
+
//- Correct for topo change.
void updateMesh(const mapPolyMesh&)
{
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index e219bc95abd78046382169ef4346b893c6ac6cc0..8797519d4dfc239511dd3b5c5351d86e0071e543 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -185,7 +185,7 @@ void Foam::mapDistribute::distribute
{
if (!contiguous<T>())
{
- PstreamBuffers pBuffs(Pstream::nonBlocking);
+ PstreamBuffers pBufs(Pstream::nonBlocking);
// Stream data into buffer
for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -195,13 +195,13 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size())
{
// Put data into send buffer
- UOPstream toDomain(domain, pBuffs);
+ UOPstream toDomain(domain, pBufs);
toDomain << UIndirectList<T>(field, map);
}
}
// Start receiving
- pBuffs.finishedSends();
+ pBufs.finishedSends();
{
// Set up 'send' to myself
@@ -231,7 +231,7 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size())
{
- UIPstream str(domain, pBuffs);
+ UIPstream str(domain, pBufs);
List<T> recvField(str);
if (recvField.size() != map.size())
@@ -551,7 +551,7 @@ void Foam::mapDistribute::distribute
{
if (!contiguous<T>())
{
- PstreamBuffers pBuffs(Pstream::nonBlocking);
+ PstreamBuffers pBufs(Pstream::nonBlocking);
// Stream data into buffer
for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -561,13 +561,13 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size())
{
// Put data into send buffer
- UOPstream toDomain(domain, pBuffs);
+ UOPstream toDomain(domain, pBufs);
toDomain << UIndirectList<T>(field, map);
}
}
// Start receiving
- pBuffs.finishedSends();
+ pBufs.finishedSends();
{
// Set up 'send' to myself
@@ -597,7 +597,7 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size())
{
- UIPstream str(domain, pBuffs);
+ UIPstream str(domain, pBufs);
List<T> recvField(str);
if (recvField.size() != map.size())
@@ -757,4 +757,66 @@ void Foam::mapDistribute::distribute
}
+template<class T>
+void Foam::mapDistribute::send(PstreamBuffers& pBufs, const List<T>& field)
+const
+{
+ // Stream data into buffer
+ for (label domain = 0; domain < Pstream::nProcs(); domain++)
+ {
+ const labelList& map = subMap_[domain];
+
+ if (map.size())
+ {
+ // Put data into send buffer
+ UOPstream toDomain(domain, pBufs);
+ toDomain << UIndirectList<T>(field, map);
+ }
+ }
+
+ // Start sending and receiving but do not block.
+ pBufs.finishedSends(false);
+}
+
+
+template<class T>
+void Foam::mapDistribute::receive(PstreamBuffers& pBufs, List<T>& field) const
+{
+ // Consume
+ field.setSize(constructSize_);
+
+ for (label domain = 0; domain < Pstream::nProcs(); domain++)
+ {
+ const labelList& map = constructMap_[domain];
+
+ if (map.size())
+ {
+ UIPstream str(domain, pBufs);
+ List<T> recvField(str);
+
+ if (recvField.size() != map.size())
+ {
+ FatalErrorIn
+ (
+ "template<class T>\n"
+ "void mapDistribute::receive\n"
+ "(\n"
+ " PstreamBuffers&,\n"
+ " List<T>&\n"
+ ")\n"
+ ) << "Expected from processor " << domain
+ << " " << map.size() << " but received "
+ << recvField.size() << " elements."
+ << abort(FatalError);
+ }
+
+ forAll(map, i)
+ {
+ field[map[i]] = recvField[i];
+ }
+ }
+ }
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 7e7272b03236f8512b9457708fc1bd4fd87a15f0..718df480e00713c8991728542e35faa1bb3c4f22 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -291,6 +291,9 @@ class primitiveMesh
//- Skewness warning threshold
static scalar skewThreshold_;
+ //- Threshold where faces are considered coplanar
+ static scalar planarCosAngle_;
+
protected:
@@ -646,6 +649,13 @@ public:
labelHashSet* setPtr = NULL
) const;
+ //- Check for concave cells by the planes of faces
+ bool checkConcaveCells
+ (
+ const bool report = false,
+ labelHashSet* setPtr = NULL
+ ) const;
+
//- Check mesh topology for correctness.
// Returns false for no error.
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 30a851e237913f7d97e1ee3066d3b118160d2563..f28e9affd09efcd72ac43f738b3de241c4669536 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -37,6 +37,7 @@ Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg
Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4;
+Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -2093,6 +2094,121 @@ bool Foam::primitiveMesh::checkCellDeterminant
}
+bool Foam::primitiveMesh::checkConcaveCells
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ if (debug)
+ {
+ Info<< "bool primitiveMesh::checkConcaveCells(const bool"
+ << ", labelHashSet*) const: "
+ << "checking for concave cells" << endl;
+ }
+
+ const cellList& c = cells();
+ const labelList& fOwner = faceOwner();
+ const vectorField& fAreas = faceAreas();
+ const pointField& fCentres = faceCentres();
+
+ label nConcaveCells = 0;
+
+ forAll (c, cellI)
+ {
+ const cell& cFaces = c[cellI];
+
+ bool concave = false;
+
+ forAll(cFaces, i)
+ {
+ if (concave)
+ {
+ break;
+ }
+
+ label fI = cFaces[i];
+
+ const point& fC = fCentres[fI];
+
+ vector fN = fAreas[fI];
+
+ fN /= max(mag(fN), VSMALL);
+
+ // Flip normal if required so that it is always pointing out of
+ // the cell
+ if (fOwner[fI] != cellI)
+ {
+ fN *= -1;
+ }
+
+ // Is the centre of any other face of the cell on the
+ // wrong side of the plane of this face?
+
+ forAll(cFaces, j)
+ {
+ if (j != i)
+ {
+ label fJ = cFaces[j];
+
+ const point& pt = fCentres[fJ];
+
+ // If the cell is concave, the point will be on the
+ // positive normal side of the plane of f, defined by
+ // its centre and normal, and the angle between (pt -
+ // fC) and fN will be less than 90 degrees, so the dot
+ // product will be positive.
+
+ vector pC = (pt - fC);
+
+ pC /= max(mag(pC), VSMALL);
+
+ if ((pC & fN) > -planarCosAngle_)
+ {
+ // Concave or planar face
+
+ concave = true;
+
+ if (setPtr)
+ {
+ setPtr->insert(cellI);
+ }
+
+ nConcaveCells++;
+
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ reduce(nConcaveCells, sumOp<label>());
+
+ if (nConcaveCells > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Concave cells (using face planes) found,"
+ << " number of cells: " << nConcaveCells << endl;
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " Concave cell check OK." << endl;
+ }
+
+ return false;
+ }
+
+ return false;
+}
+
+
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::primitiveMesh::checkTopology(const bool report) const
diff --git a/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C b/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C
index 9fdfee1b0a5943d724a855226c3ae829790bdceb..695d211d61e71c63a1b20a75d419d7f780135ef6 100644
--- a/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C
+++ b/src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -57,8 +57,12 @@ StringListType Foam::subsetMatchingStrings
const bool invert
)
{
+ // Create copy
StringListType newLst(lst.size());
+ // ensure consistent addressable size (eg, DynamicList)
+ newLst.setSize(lst.size());
+
label nElem = 0;
forAll(lst, elemI)
{
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
index 07ceb696108defb4890df7a81f9a538c16220517..92ba164bb82de00df1319b2d2f31f5a869585d7c 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
@@ -202,12 +202,45 @@ Foam::label Foam::removePoints::countPointUsage
pointCanBeDeleted[pointI] = true;
nDeleted++;
}
-
}
edge0.clear();
edge1.clear();
+ // Protect any points on faces that would collapse down to nothing
+ // No particular intelligence so might protect too many points
+ forAll(mesh_.faces(), faceI)
+ {
+ const face& f = mesh_.faces()[faceI];
+
+ label nCollapse = 0;
+ forAll(f, fp)
+ {
+ if (pointCanBeDeleted[f[fp]])
+ {
+ nCollapse++;
+ }
+ }
+
+ if ((f.size() - nCollapse) < 3)
+ {
+ // Just unmark enough points
+ forAll(f, fp)
+ {
+ if (pointCanBeDeleted[f[fp]])
+ {
+ pointCanBeDeleted[f[fp]] = false;
+ --nCollapse;
+ if (nCollapse == 0)
+ {
+ break;
+ }
+ }
+ }
+ }
+ }
+
+
// Point can be deleted only if all processors want to delete it
syncTools::syncPointList
(
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 828e7e65ff10b7f3071b5974af321f8c08e3bc12..50942e39c72629a59aa17dd957a4f6c4fa0f487b 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -1375,7 +1375,7 @@ void Foam::autoLayerDriver::growNoExtrusion
pointField& patchDisp,
labelList& patchNLayers,
List<extrudeMode>& extrudeStatus
-)
+) const
{
Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
@@ -1406,7 +1406,7 @@ void Foam::autoLayerDriver::growNoExtrusion
{
if
(
- extrudeStatus[f[fp]] == NOEXTRUDE
+ extrudeStatus[f[fp]] == EXTRUDE
&& grownExtrudeStatus[f[fp]] != NOEXTRUDE
)
{
@@ -1419,6 +1419,31 @@ void Foam::autoLayerDriver::growNoExtrusion
extrudeStatus.transfer(grownExtrudeStatus);
+
+ // Synchronise since might get called multiple times.
+ // Use the fact that NOEXTRUDE is the minimum value.
+ {
+ labelList status(extrudeStatus.size());
+ forAll(status, i)
+ {
+ status[i] = extrudeStatus[i];
+ }
+ syncTools::syncPointList
+ (
+ meshRefiner_.mesh(),
+ pp.meshPoints(),
+ status,
+ minEqOp<label>(),
+ labelMax, // null value
+ false // no separation
+ );
+ forAll(status, i)
+ {
+ extrudeStatus[i] = extrudeMode(status[i]);
+ }
+ }
+
+
forAll(extrudeStatus, patchPointI)
{
if (extrudeStatus[patchPointI] == NOEXTRUDE)
@@ -2711,8 +2736,6 @@ void Foam::autoLayerDriver::addLayers
const labelList& meshPoints = patches[patchI].meshPoints();
- //scalar maxThickness = -VGREAT;
- //scalar minThickness = VGREAT;
scalar sumThickness = 0;
scalar sumNearWallThickness = 0;
@@ -2720,8 +2743,6 @@ void Foam::autoLayerDriver::addLayers
{
label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
- //maxThickness = max(maxThickness, thickness[ppPointI]);
- //minThickness = min(minThickness, thickness[ppPointI]);
sumThickness += thickness[ppPointI];
label nLay = patchNLayers[ppPointI];
@@ -2749,8 +2770,6 @@ void Foam::autoLayerDriver::addLayers
if (totNPoints > 0)
{
- //reduce(maxThickness, maxOp<scalar>());
- //reduce(minThickness, minOp<scalar>());
avgThickness =
returnReduce(sumThickness, sumOp<scalar>())
/ totNPoints;
@@ -2766,8 +2785,6 @@ void Foam::autoLayerDriver::addLayers
<< " " << setw(6) << layerParams.numLayers()[patchI]
<< " " << setw(8) << avgNearWallThickness
<< " " << setw(8) << avgThickness
- //<< " " << setw(8) << minThickness
- //<< " " << setw(8) << maxThickness
<< endl;
}
Info<< endl;
@@ -3147,6 +3164,19 @@ void Foam::autoLayerDriver::addLayers
meshMover().movePoints(oldPoints);
meshMover().correct();
+
+ // Grow out region of non-extrusion
+ for (label i = 0; i < layerParams.nGrow(); i++)
+ {
+ growNoExtrusion
+ (
+ pp,
+ patchDisp,
+ patchNLayers,
+ extrudeStatus
+ );
+ }
+
Info<< endl;
}
@@ -3293,7 +3323,6 @@ void Foam::autoLayerDriver::doLayers
// Balance
- if (Pstream::parRun())
if (Pstream::parRun() && preBalance)
{
Info<< nl
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index b49a6934b87f9332f168bf6bac543b52944b540e..915614aad91113fd88289bea54c0186f81c85182 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -232,13 +232,13 @@ class autoLayerDriver
) const;
//- Grow no-extrusion layer.
- static void growNoExtrusion
+ void growNoExtrusion
(
const indirectPrimitivePatch& pp,
pointField& patchDisp,
labelList& patchNLayers,
List<extrudeMode>& extrudeStatus
- );
+ ) const;
//- Calculate pointwise wanted and minimum thickness.
// thickness: wanted thickness
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 882da504f9f4a3c1f86773b7471856fa1498dcd6..0d9cecd039d5c852a972049b07b1091571e18367 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -245,27 +245,8 @@ void Foam::meshRefinement::getBafflePatches
const pointField& cellCentres = mesh_.cellCentres();
- // Build list of surfaces that are not to be baffled.
- const wordList& cellZoneNames = surfaces_.cellZoneNames();
-
- labelList surfacesToBaffle(cellZoneNames.size());
- label baffleI = 0;
- forAll(cellZoneNames, surfI)
- {
- if (cellZoneNames[surfI].size())
- {
- if (debug)
- {
- Pout<< "getBafflePatches : Not baffling surface "
- << surfaces_.names()[surfI] << endl;
- }
- }
- else
- {
- surfacesToBaffle[baffleI++] = surfI;
- }
- }
- surfacesToBaffle.setSize(baffleI);
+ // Surfaces that need to be baffled
+ const labelList surfacesToBaffle(surfaces_.getUnnamedSurfaces());
ownPatch.setSize(mesh_.nFaces());
ownPatch = -1;
@@ -1439,6 +1420,15 @@ void Foam::meshRefinement::makeConsistentFaceIndex
}
}
}
+ else
+ {
+ // Unzonify boundary faces
+ forAll(pp, i)
+ {
+ label faceI = pp.start()+i;
+ namedSurfaceIndex[faceI] = -1;
+ }
+ }
}
}
@@ -2042,14 +2032,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
labelList namedSurfaces(surfaces_.getNamedSurfaces());
- boolList isNamedSurface(cellZoneNames.size(), false);
-
forAll(namedSurfaces, i)
{
label surfI = namedSurfaces[i];
- isNamedSurface[surfI] = true;
-
Info<< "Surface : " << surfaces_.names()[surfI] << nl
<< " faceZone : " << faceZoneNames[surfI] << nl
<< " cellZone : " << cellZoneNames[surfI] << endl;
@@ -2125,31 +2111,34 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
{
label surfI = namedSurfaces[i];
- label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
-
- if (zoneI == -1)
+ if (cellZoneNames[surfI] != word::null)
{
- zoneI = cellZones.size();
- cellZones.setSize(zoneI+1);
- cellZones.set
- (
- zoneI,
- new cellZone
+ label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
+
+ if (zoneI == -1)
+ {
+ zoneI = cellZones.size();
+ cellZones.setSize(zoneI+1);
+ cellZones.set
(
- cellZoneNames[surfI], //name
- labelList(0), //addressing
- zoneI, //index
- cellZones //cellZoneMesh
- )
- );
- }
+ zoneI,
+ new cellZone
+ (
+ cellZoneNames[surfI], //name
+ labelList(0), //addressing
+ zoneI, //index
+ cellZones //cellZoneMesh
+ )
+ );
+ }
- if (debug)
- {
- Pout<< "Cells inside " << surfaces_.names()[surfI]
- << " will go into cellZone " << zoneI << endl;
+ if (debug)
+ {
+ Pout<< "Cells inside " << surfaces_.names()[surfI]
+ << " will go into cellZone " << zoneI << endl;
+ }
+ surfaceToCellZone[surfI] = zoneI;
}
- surfaceToCellZone[surfI] = zoneI;
}
// Check they are synced
@@ -2321,6 +2310,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
if (closedNamedSurfaces.size())
{
+ Info<< "Found " << closedNamedSurfaces.size()
+ << " closed, named surfaces. Assigning cells in/outside"
+ << " these surfaces to the corresponding cellZone."
+ << nl << endl;
+
findCellZoneGeometric
(
closedNamedSurfaces, // indices of closed surfaces
@@ -2333,8 +2327,12 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Set using walking
// ~~~~~~~~~~~~~~~~~
- //if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
+ //if (!allowFreeStandingZoneFaces)
{
+ Info<< "Walking from location-in-mesh " << keepPoint
+ << " to assign cellZones "
+ << "- crossing a faceZone face changes cellZone" << nl << endl;
+
// Topological walk
findCellZoneTopo
(
@@ -2349,6 +2347,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
if (!allowFreeStandingZoneFaces)
{
+ Info<< "Only keeping zone faces inbetween different cellZones."
+ << nl << endl;
+
makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
}
@@ -2521,6 +2522,32 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
mesh_.setInstance(oldInstance());
}
+ // Print some stats (note: zones are synchronised)
+ if (mesh_.cellZones().size() > 0)
+ {
+ Info<< "CellZones:" << endl;
+ forAll(mesh_.cellZones(), zoneI)
+ {
+ const cellZone& cz = mesh_.cellZones()[zoneI];
+ Info<< " " << cz.name()
+ << "\tsize:" << returnReduce(cz.size(), sumOp<label>())
+ << endl;
+ }
+ Info<< endl;
+ }
+ if (mesh_.faceZones().size() > 0)
+ {
+ Info<< "FaceZones:" << endl;
+ forAll(mesh_.faceZones(), zoneI)
+ {
+ const faceZone& fz = mesh_.faceZones()[zoneI];
+ Info<< " " << fz.name()
+ << "\tsize:" << returnReduce(fz.size(), sumOp<label>())
+ << endl;
+ }
+ Info<< endl;
+ }
+
// None of the faces has changed, only the zones. Still...
updateMesh(map, labelList());
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index d3687b7bdbfcb09b3b46f5142437dd3335ffa5c4..ac01e052cdb6b4a8062505b32f9a7c8ac7dfa16d 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -77,8 +77,29 @@ Foam::refinementSurfaces::refinementSurfaces
if (dict.found("faceZone"))
{
dict.lookup("faceZone") >> faceZoneNames_[surfI];
- dict.lookup("cellZone") >> cellZoneNames_[surfI];
- dict.lookup("zoneInside") >> zoneInside_[surfI];
+ bool hasSide = dict.readIfPresent("zoneInside", zoneInside_[surfI]);
+ if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
+ {
+ if (hasSide && !allGeometry_[surfaces_[surfI]].hasVolumeType())
+ {
+ IOWarningIn
+ (
+ "refinementSurfaces::refinementSurfaces(..)",
+ dict
+ ) << "Unused entry zoneInside for faceZone "
+ << faceZoneNames_[surfI]
+ << " since surface " << names_[surfI]
+ << " is not closed." << endl;
+ }
+ }
+ else if (hasSide)
+ {
+ IOWarningIn("refinementSurfaces::refinementSurfaces(..)", dict)
+ << "Unused entry zoneInside for faceZone "
+ << faceZoneNames_[surfI]
+ << " since no cellZone specified."
+ << endl;
+ }
}
// Global perpendicular angle
@@ -314,8 +335,40 @@ Foam::refinementSurfaces::refinementSurfaces
if (dict.found("faceZone"))
{
dict.lookup("faceZone") >> faceZoneNames_[surfI];
- dict.lookup("cellZone") >> cellZoneNames_[surfI];
- dict.lookup("zoneInside") >> zoneInside_[surfI];
+ bool hasSide = dict.readIfPresent
+ (
+ "zoneInside",
+ zoneInside_[surfI]
+ );
+ if (dict.readIfPresent("cellZone", cellZoneNames_[surfI]))
+ {
+ if
+ (
+ hasSide
+ && !allGeometry_[surfaces_[surfI]].hasVolumeType()
+ )
+ {
+ IOWarningIn
+ (
+ "refinementSurfaces::refinementSurfaces(..)",
+ dict
+ ) << "Unused entry zoneInside for faceZone "
+ << faceZoneNames_[surfI]
+ << " since surface " << names_[surfI]
+ << " is not closed." << endl;
+ }
+ }
+ else if (hasSide)
+ {
+ IOWarningIn
+ (
+ "refinementSurfaces::refinementSurfaces(..)",
+ dict
+ ) << "Unused entry zoneInside for faceZone "
+ << faceZoneNames_[surfI]
+ << " since no cellZone specified."
+ << endl;
+ }
}
// Global perpendicular angle
@@ -475,18 +528,17 @@ Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
// Get indices of closed named surfaces
Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
{
- labelList named(getNamedSurfaces());
+ labelList closed(cellZoneNames_.size());
- labelList closed(named.size());
label closedI = 0;
-
- forAll(named, i)
+ forAll(cellZoneNames_, surfI)
{
- label surfI = named[i];
-
- if (allGeometry_[surfaces_[surfI]].hasVolumeType())
+ if (cellZoneNames_[surfI].size())
{
- closed[closedI++] = surfI;
+ if (allGeometry_[surfaces_[surfI]].hasVolumeType())
+ {
+ closed[closedI++] = surfI;
+ }
}
}
closed.setSize(closedI);
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index 4a7b6f96c9992d47c603099cc67145db504c121c..bafe949a08adc469ac611a8d6060410b9c8096b5 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -146,7 +146,8 @@ public:
return faceZoneNames_;
}
- //- Per 'interface' surface : name of cellZone to put cells into
+ //- Per 'interface' surface : empty or name of cellZone to put
+ // cells into
const wordList& cellZoneNames() const
{
return cellZoneNames_;
@@ -158,7 +159,7 @@ public:
//- Get indices of named surfaces (surfaces with faceZoneName)
labelList getNamedSurfaces() const;
- //- Get indices of closed named surfaces
+ //- Get indices of closed surfaces with a cellZone
labelList getClosedNamedSurfaces() const;
//- From local region number to global region number
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
index 166d74ae2a77b03eb7065f793184e5c15de74477..91c24ecf28f07ca09c8d3538b406afdf95a34777 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
@@ -637,7 +637,7 @@ Foam::directMappedPatchBase::directMappedPatchBase
samplePatch_(samplePatch),
uniformOffset_(true),
offset_(offset),
- offsets_(0),
+ offsets_(pp.size(), offset_),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mapPtr_(NULL)
{}
@@ -670,7 +670,7 @@ Foam::directMappedPatchBase::directMappedPatchBase
offsets_
(
uniformOffset_
- ? pointField(patch_.size(), offset_)
+ ? pointField(pp.size(), offset_)
: dict.lookup("offsets")
),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
diff --git a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
index cf8f8a07ce65b679bef3184432cdb2c7c66167b3..d2759d29274673be2ce3fa1f94f3655bdb2699ba 100644
--- a/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
+++ b/src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
@@ -33,7 +33,7 @@ Description
"Turbulence Modeling for CFD"
D. C. Wilcox,
DCW Industries, Inc., La Canada,
- California, 1998.
+ California, 1988.
See also:
http://www.cfd-online.com/Wiki/Wilcox's_k-omega_model
diff --git a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
index c727a74d65c6623604098287863c4f8b159a1fec..0941f0352873a68295be41f4b44baf60fba1a797 100644
--- a/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
+++ b/tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
@@ -64,7 +64,6 @@ int main(int argc, char *argv[])
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- - fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U)
);
mrfZones.addCoriolis(UEqn());