zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2010-04-07 11:42:32 +01:00
parent 149a190e22 90189b35d4
commit 357ff92e8d
35 changed files with 897 additions and 1017 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
COPYING
View File

@ -3,7 +3,7 @@
OpenFOAM(R) is Copyright (C) 1991-2010 OpenCFD Ltd. OpenFOAM(R) is Copyright (C) 1991-2010 OpenCFD Ltd.
Contact: OpenCFD (enquiries@OpenCFD.co.uk) 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 of GNU General Public License version 3, which is displayed below, or
(at your option) any later version. (at your option) any later version.

1
applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
View File

@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U) + turbulence->divDevRhoReff(U)
); );

1
applications/solvers/compressible/rhoSimpleFoam/UEqn.H
View File

@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U) + turbulence->divDevRhoReff(U)
); );

9
applications/solvers/compressible/rhoSimpleFoam/createFields.H
View File

@ -44,9 +44,14 @@
scalar pRefValue = 0.0; scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue); 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; Info<< "Creating turbulence model\n" << endl;

4
applications/solvers/compressible/rhoSimpleFoam/hEqn.H
View File

@ -5,8 +5,8 @@
- fvm::Sp(fvc::div(phi), h) - fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h) - fvm::laplacian(turbulence->alphaEff(), h)
== ==
fvc::div(phi/fvc::interpolate(rho), p, "div(U,p)") fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
- p*fvc::div(phi/fvc::interpolate(rho)) - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
); );
hEqn.relax(); hEqn.relax();

15
applications/solvers/compressible/rhoSimpleFoam/pEqn.H
View File

@ -1,4 +1,7 @@
rho = thermo.rho(); rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rUA = 1.0/UEqn().A(); volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H(); U = rUA*UEqn().H();
@ -82,15 +85,9 @@ else
// Explicitly relax pressure for momentum corrector // Explicitly relax pressure for momentum corrector
p.relax(); p.relax();
rho = thermo.rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
U -= rUA*fvc::grad(p); U -= rUA*fvc::grad(p);
U.correctBoundaryConditions(); U.correctBoundaryConditions();
bound(p, pMin);
// For closed-volume cases adjust the pressure and density levels // For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity // to obey overall mass continuity
if (closedVolume) if (closedVolume)
@ -98,3 +95,9 @@ if (closedVolume)
p += (initialMass - fvc::domainIntegrate(psi*p)) p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi); /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;

1
applications/solvers/compressible/rhoSimplecFoam/UEqn.H
View File

@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U) + turbulence->divDevRhoReff(U)
); );

1
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
View File

@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U) + turbulence->divDevReff(U)
); );

1
applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
View File

@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U) + turbulence->divDevRhoReff(U)
); );

1
applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
View File

@ -1,7 +1,6 @@
regionProperties/regionProperties.C regionProperties/regionProperties.C
derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
fluid/compressibleCourantNo.C fluid/compressibleCourantNo.C
solid/solidRegionDiffNo.C solid/solidRegionDiffNo.C

381
applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
View File

@ -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
// ************************************************************************* //

192
applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H
View File

@ -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
// ************************************************************************* //

4
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
View File

@ -1,6 +1,5 @@
scalar CoNum = -GREAT; scalar CoNum = -GREAT;
if (fluidRegions.size())
{
forAll(fluidRegions, regionI) forAll(fluidRegions, regionI)
{ {
CoNum = max CoNum = max
@ -15,4 +14,3 @@
CoNum CoNum
); );
} }
}

4
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H
View File

@ -1,6 +1,5 @@
scalar DiNum = -GREAT; scalar DiNum = -GREAT;
if (solidRegions.size())
{
forAll(solidRegions, regionI) forAll(solidRegions, regionI)
{ {
DiNum = max DiNum = max
@ -15,4 +14,3 @@
DiNum DiNum
); );
} }
}

29
applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
View File

@ -37,7 +37,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF const DimensionedField<scalar, volMesh>& iF
) )
: :
fixedValueFvPatchScalarField(p, iF), fixedGradientFvPatchScalarField(p, iF),
q_(p.size(), 0.0), q_(p.size(), 0.0),
KName_("undefined-K") KName_("undefined-K")
{} {}
@ -52,7 +52,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const fvPatchFieldMapper& mapper const fvPatchFieldMapper& mapper
) )
: :
fixedValueFvPatchScalarField(ptf, p, iF, mapper), fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
q_(ptf.q_, mapper), q_(ptf.q_, mapper),
KName_(ptf.KName_) KName_(ptf.KName_)
{} {}
@ -66,7 +66,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const dictionary& dict const dictionary& dict
) )
: :
fixedValueFvPatchScalarField(p, iF, dict), fixedGradientFvPatchScalarField(p, iF, dict),
q_("q", dict, p.size()), q_("q", dict, p.size()),
KName_(dict.lookup("K")) KName_(dict.lookup("K"))
{} {}
@ -78,7 +78,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf
) )
: :
fixedValueFvPatchScalarField(tppsf), fixedGradientFvPatchScalarField(tppsf),
q_(tppsf.q_), q_(tppsf.q_),
KName_(tppsf.KName_) KName_(tppsf.KName_)
{} {}
@ -91,7 +91,7 @@ solidWallHeatFluxTemperatureFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF const DimensionedField<scalar, volMesh>& iF
) )
: :
fixedValueFvPatchScalarField(tppsf, iF), fixedGradientFvPatchScalarField(tppsf, iF),
q_(tppsf.q_), q_(tppsf.q_),
KName_(tppsf.KName_) KName_(tppsf.KName_)
{} {}
@ -104,7 +104,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::autoMap
const fvPatchFieldMapper& m const fvPatchFieldMapper& m
) )
{ {
fixedValueFvPatchScalarField::autoMap(m); fixedGradientFvPatchScalarField::autoMap(m);
q_.autoMap(m); q_.autoMap(m);
} }
@ -115,7 +115,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
const labelList& addr const labelList& addr
) )
{ {
fixedValueFvPatchScalarField::rmap(ptf, addr); fixedGradientFvPatchScalarField::rmap(ptf, addr);
const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf = const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf =
refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf); refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf);
@ -131,14 +131,14 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
return; return;
} }
const scalarField& Kw = const scalarField& Kw = patch().lookupPatchField<volScalarField, scalar>
patch().lookupPatchField<volScalarField, scalar>(KName_); (
KName_
);
const fvPatchScalarField& Tw = *this; gradient() = q_/Kw;
operator==(q_/(patch().deltaCoeffs()*Kw) + Tw.patchInternalField()); fixedGradientFvPatchScalarField::updateCoeffs();
fixedValueFvPatchScalarField::updateCoeffs();
} }
@ -147,9 +147,10 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
Ostream& os Ostream& os
) const ) const
{ {
fixedValueFvPatchScalarField::write(os); fixedGradientFvPatchScalarField::write(os);
q_.writeEntry("q", os); q_.writeEntry("q", os);
os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl; os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
this->writeEntry("value", os);
} }

4
applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H
View File

@ -45,7 +45,7 @@ SourceFiles
#ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H #ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
#define solidWallHeatFluxTemperatureFvPatchScalarField_H #define solidWallHeatFluxTemperatureFvPatchScalarField_H
#include "fixedValueFvPatchFields.H" #include "fixedGradientFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -58,7 +58,7 @@ namespace Foam
class solidWallHeatFluxTemperatureFvPatchScalarField class solidWallHeatFluxTemperatureFvPatchScalarField
: :
public fixedValueFvPatchScalarField public fixedGradientFvPatchScalarField
{ {
// Private data // Private data

1
applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H
View File

@ -2,7 +2,6 @@
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turb.divDevRhoReff(U) + turb.divDevRhoReff(U)
); );

1
applications/solvers/incompressible/porousSimpleFoam/UEqn.H
View File

@ -3,7 +3,6 @@
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U) + turbulence->divDevReff(U)
); );

16
applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
View File

@ -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; return noFailedChecks;
} }

662
applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
View File

@ -39,17 +39,29 @@ Description
- mesh with cells put into cellZones (-makeCellZones) - mesh with cells put into cellZones (-makeCellZones)
Note: 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. - Should work in parallel.
cellZones can differ on either side of processor boundaries in which case cellZones can differ on either side of processor boundaries in which case
the faces get moved from processor patch to directMapped patch. Not 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 - If a cell zone gets split into more than one region it can detect
the largest matching region (-sloppyCellZones). This will accept any the largest matching region (-sloppyCellZones). This will accept any
region that covers more than 50% of the zone. It has to be a subset region that covers more than 50% of the zone. It has to be a subset
so cannot have any cells in any other zone. 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: - writes maps like decomposePar back to original mesh:
- pointRegionAddressing : for every point in this region the point in - pointRegionAddressing : for every point in this region the point in
the original mesh 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 // Find region that covers most of cell zone
label findCorrespondingRegion label findCorrespondingRegion
( (
@ -1314,62 +1269,20 @@ label findCorrespondingRegion
//} //}
// Main program: // Get zone per cell
int main(int argc, char *argv[])
{
# 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 // - non-unique zoning
// - coupled zones // - coupled zones
void getZoneID
(
const polyMesh& mesh,
const cellZoneMesh& cellZones,
labelList& zoneID,
labelList& neiZoneID
)
{
// Existing zoneID // Existing zoneID
labelList zoneID(mesh.nCells(), -1); zoneID.setSize(mesh.nCells());
zoneID = -1;
forAll(cellZones, zoneI) forAll(cellZones, zoneI)
{ {
@ -1384,7 +1297,7 @@ int main(int argc, char *argv[])
} }
else else
{ {
FatalErrorIn(args.executable()) FatalErrorIn("getZoneID(..)")
<< "Cell " << cellI << " with cell centre " << "Cell " << cellI << " with cell centre "
<< mesh.cellCentres()[cellI] << mesh.cellCentres()[cellI]
<< " is multiple zones. This is not allowed." << endl << " is multiple zones. This is not allowed." << endl
@ -1396,175 +1309,42 @@ int main(int argc, char *argv[])
} }
// Neighbour zoneID. // Neighbour zoneID.
labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces()); neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
forAll(neiZoneID, i) forAll(neiZoneID, i)
{ {
neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]]; neiZoneID[i] = zoneID[mesh.faceOwner()[i+mesh.nInternalFaces()]];
} }
syncTools::swapBoundaryFaceList(mesh, neiZoneID, false); syncTools::swapBoundaryFaceList(mesh, neiZoneID, false);
// Determine connected regions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mark additional faces that are blocked
boolList blockedFace;
// Read from faceSet
if (blockedFacesName.size())
{
faceSet blockedFaceSet(mesh, blockedFacesName);
Info<< "Read " << returnReduce(blockedFaceSet.size(), sumOp<label>())
<< " blocked faces from set " << blockedFacesName << nl << endl;
blockedFace.setSize(mesh.nFaces(), false);
forAllConstIter(faceSet, blockedFaceSet, iter)
{
blockedFace[iter.key()] = true;
}
}
// Imply from differing cellZones
if (useCellZones)
{
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;
}
}
} }
// Determine per cell the region it belongs to void matchRegions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// cellRegion is the labelList with the region per cell.
labelList cellRegion;
label nCellRegions = 0;
if (useCellZonesOnly)
{
label unzonedCellI = findIndex(zoneID, -1);
if (unzonedCellI != -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);
}
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;
// Write to manual decomposition option
{
labelIOList cellToRegion
( (
IOobject const bool sloppyCellZones,
( const polyMesh& mesh,
"cellToRegion",
mesh.facesInstance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
cellRegion
);
cellToRegion.write();
Info<< "Writing region per cell file (for manual decomposition) to " const label nCellRegions,
<< cellToRegion.objectPath() << nl << endl; const labelList& cellRegion,
}
// Write for postprocessing labelList& regionToZone,
wordList& regionNames,
labelList& zoneToRegion
)
{ {
volScalarField cellToRegion const cellZoneMesh& cellZones = mesh.cellZones();
(
IOobject
(
"cellToRegion",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
zeroGradientFvPatchScalarField::typeName
);
forAll(cellRegion, cellI)
{
cellToRegion[cellI] = cellRegion[cellI];
}
cellToRegion.write();
Info<< "Writing region per cell as volScalarField to " regionToZone.setSize(nCellRegions, -1);
<< cellToRegion.objectPath() << nl << endl; regionNames.setSize(nCellRegions);
} zoneToRegion.setSize(cellZones.size(), -1);
// 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;
// Get current per cell zoneID
labelList zoneID(mesh.nCells(), -1);
labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
getZoneID(mesh, cellZones, zoneID, neiZoneID);
// Sizes per cellzone // Sizes per cellzone
// ~~~~~~~~~~~~~~~~~~
labelList zoneSizes(cellZones.size(), 0); labelList zoneSizes(cellZones.size(), 0);
if (useCellZones || makeCellZones || sloppyCellZones)
{ {
List<wordList> zoneNames(Pstream::nProcs()); List<wordList> zoneNames(Pstream::nProcs());
zoneNames[Pstream::myProcNo()] = cellZones.names(); zoneNames[Pstream::myProcNo()] = cellZones.names();
@ -1575,7 +1355,7 @@ int main(int argc, char *argv[])
{ {
if (zoneNames[procI] != zoneNames[0]) if (zoneNames[procI] != zoneNames[0])
{ {
FatalErrorIn(args.executable()) FatalErrorIn("matchRegions(..)")
<< "cellZones not synchronised across processors." << endl << "cellZones not synchronised across processors." << endl
<< "Master has cellZones " << zoneNames[0] << endl << "Master has cellZones " << zoneNames[0] << endl
<< "Processor " << procI << "Processor " << procI
@ -1595,16 +1375,6 @@ int main(int argc, char *argv[])
} }
// Whether region corresponds to a cellzone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Region per zone
labelList regionToZone(nCellRegions, -1);
// Name of region
wordList regionNames(nCellRegions);
// Zone to region
labelList zoneToRegion(cellZones.size(), -1);
if (sloppyCellZones) if (sloppyCellZones)
{ {
Info<< "Trying to match regions to existing cell zones;" Info<< "Trying to match regions to existing cell zones;"
@ -1624,7 +1394,7 @@ int main(int argc, char *argv[])
if (regionI != -1) if (regionI != -1)
{ {
Info<< "Sloppily matched region " << regionI Info<< "Sloppily matched region " << regionI
<< " size " << regionSizes[regionI] //<< " size " << regionSizes[regionI]
<< " to zone " << zoneI << " size " << zoneSizes[zoneI] << " to zone " << zoneI << " size " << zoneSizes[zoneI]
<< endl; << endl;
zoneToRegion[zoneI] = regionI; zoneToRegion[zoneI] = regionI;
@ -1664,6 +1434,330 @@ int main(int argc, char *argv[])
regionNames[regionI] = "domain" + Foam::name(regionI); regionNames[regionI] = "domain" + Foam::name(regionI);
} }
} }
}
void writeCellToRegion(const fvMesh& mesh, const labelList& cellRegion)
{
// Write to manual decomposition option
{
labelIOList cellToRegion
(
IOobject
(
"cellToRegion",
mesh.facesInstance(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
cellRegion
);
cellToRegion.write();
Info<< "Writing region per cell file (for manual decomposition) to "
<< cellToRegion.objectPath() << nl << endl;
}
// Write for postprocessing
{
volScalarField cellToRegion
(
IOobject
(
"cellToRegion",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
zeroGradientFvPatchScalarField::typeName
);
forAll(cellRegion, cellI)
{
cellToRegion[cellI] = cellRegion[cellI];
}
cellToRegion.write();
Info<< "Writing region per cell as volScalarField to "
<< cellToRegion.objectPath() << nl << endl;
}
}
// Main program:
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))
{
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 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
)
)
{
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);
}
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();
// Existing zoneID
labelList zoneID(mesh.nCells(), -1);
// Neighbour zoneID.
labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
getZoneID(mesh, cellZones, zoneID, neiZoneID);
// Determine per cell the region it belongs to
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// cellRegion is the labelList with the region per cell.
labelList cellRegion;
// Region per zone
labelList regionToZone;
// Name of region
wordList regionNames;
// Zone to region
labelList zoneToRegion;
label nCellRegions = 0;
if (useCellZonesOnly)
{
Info<< "Using current cellZones to split mesh into regions."
<< " This requires all"
<< " cells to be in one and only one cellZone." << nl << endl;
label unzonedCellI = findIndex(zoneID, -1);
if (unzonedCellI != -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);
}
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
(
zoneFile,
mesh.facesInstance(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
mesh
);
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
{
// Determine connected regions
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Mark additional faces that are blocked
boolList blockedFace;
// Read from faceSet
if (blockedFacesName.size())
{
faceSet blockedFaceSet(mesh, blockedFacesName);
Info<< "Read "
<< returnReduce(blockedFaceSet.size(), sumOp<label>())
<< " blocked faces from set " << blockedFacesName << nl << endl;
blockedFace.setSize(mesh.nFaces(), false);
forAllConstIter(faceSet, blockedFaceSet, iter)
{
blockedFace[iter.key()] = true;
}
}
// Imply from differing cellZones
if (useCellZones)
{
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 // Print region to zone
@ -1676,16 +1770,6 @@ int main(int argc, char *argv[])
} }
Info<< endl; 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 // Since we're going to mess with patches make sure all non-processor ones

17
src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -37,6 +37,9 @@ ListType Foam::renumber
// Create copy // Create copy
ListType newLst(lst.size()); ListType newLst(lst.size());
// ensure consistent addressable size (eg, DynamicList)
newLst.setSize(lst.size());
forAll(lst, elemI) forAll(lst, elemI)
{ {
if (lst[elemI] >= 0) if (lst[elemI] >= 0)
@ -76,6 +79,9 @@ ListType Foam::reorder
// Create copy // Create copy
ListType newLst(lst.size()); ListType newLst(lst.size());
// ensure consistent addressable size (eg, DynamicList)
newLst.setSize(lst.size());
forAll(lst, elemI) forAll(lst, elemI)
{ {
if (oldToNew[elemI] >= 0) if (oldToNew[elemI] >= 0)
@ -101,6 +107,9 @@ void Foam::inplaceReorder
// Create copy // Create copy
ListType newLst(lst.size()); ListType newLst(lst.size());
// ensure consistent addressable size (eg, DynamicList)
newLst.setSize(lst.size());
forAll(lst, elemI) forAll(lst, elemI)
{ {
if (oldToNew[elemI] >= 0) if (oldToNew[elemI] >= 0)
@ -258,6 +267,9 @@ ListType Foam::subset
ListType newLst(lst.size()); ListType newLst(lst.size());
// ensure consistent addressable size (eg, DynamicList)
newLst.setSize(lst.size());
label nElem = 0; label nElem = 0;
forAll(lst, elemI) forAll(lst, elemI)
{ {
@ -318,6 +330,9 @@ ListType Foam::subset
ListType newLst(lst.size()); ListType newLst(lst.size());
// ensure consistent addressable size (eg, DynamicList)
newLst.setSize(lst.size());
label nElem = 0; label nElem = 0;
forAll(lst, elemI) forAll(lst, elemI)
{ {

8
src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
View File

@ -70,6 +70,7 @@ namespace Foam
class mapPolyMesh; class mapPolyMesh;
class globalIndex; class globalIndex;
class PstreamBuffers;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class mapDistribute Declaration 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. //- Correct for topo change.
void updateMesh(const mapPolyMesh&) void updateMesh(const mapPolyMesh&)
{ {

78
src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
View File

@ -185,7 +185,7 @@ void Foam::mapDistribute::distribute
{ {
if (!contiguous<T>()) if (!contiguous<T>())
{ {
PstreamBuffers pBuffs(Pstream::nonBlocking); PstreamBuffers pBufs(Pstream::nonBlocking);
// Stream data into buffer // Stream data into buffer
for (label domain = 0; domain < Pstream::nProcs(); domain++) for (label domain = 0; domain < Pstream::nProcs(); domain++)
@ -195,13 +195,13 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size()) if (domain != Pstream::myProcNo() && map.size())
{ {
// Put data into send buffer // Put data into send buffer
UOPstream toDomain(domain, pBuffs); UOPstream toDomain(domain, pBufs);
toDomain << UIndirectList<T>(field, map); toDomain << UIndirectList<T>(field, map);
} }
} }
// Start receiving // Start receiving
pBuffs.finishedSends(); pBufs.finishedSends();
{ {
// Set up 'send' to myself // Set up 'send' to myself
@ -231,7 +231,7 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size()) if (domain != Pstream::myProcNo() && map.size())
{ {
UIPstream str(domain, pBuffs); UIPstream str(domain, pBufs);
List<T> recvField(str); List<T> recvField(str);
if (recvField.size() != map.size()) if (recvField.size() != map.size())
@ -551,7 +551,7 @@ void Foam::mapDistribute::distribute
{ {
if (!contiguous<T>()) if (!contiguous<T>())
{ {
PstreamBuffers pBuffs(Pstream::nonBlocking); PstreamBuffers pBufs(Pstream::nonBlocking);
// Stream data into buffer // Stream data into buffer
for (label domain = 0; domain < Pstream::nProcs(); domain++) for (label domain = 0; domain < Pstream::nProcs(); domain++)
@ -561,13 +561,13 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size()) if (domain != Pstream::myProcNo() && map.size())
{ {
// Put data into send buffer // Put data into send buffer
UOPstream toDomain(domain, pBuffs); UOPstream toDomain(domain, pBufs);
toDomain << UIndirectList<T>(field, map); toDomain << UIndirectList<T>(field, map);
} }
} }
// Start receiving // Start receiving
pBuffs.finishedSends(); pBufs.finishedSends();
{ {
// Set up 'send' to myself // Set up 'send' to myself
@ -597,7 +597,7 @@ void Foam::mapDistribute::distribute
if (domain != Pstream::myProcNo() && map.size()) if (domain != Pstream::myProcNo() && map.size())
{ {
UIPstream str(domain, pBuffs); UIPstream str(domain, pBufs);
List<T> recvField(str); List<T> recvField(str);
if (recvField.size() != map.size()) 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];
}
}
}
}
// ************************************************************************* // // ************************************************************************* //

10
src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
View File

@ -291,6 +291,9 @@ class primitiveMesh
//- Skewness warning threshold //- Skewness warning threshold
static scalar skewThreshold_; static scalar skewThreshold_;
//- Threshold where faces are considered coplanar
static scalar planarCosAngle_;
protected: protected:
@ -646,6 +649,13 @@ public:
labelHashSet* setPtr = NULL labelHashSet* setPtr = NULL
) const; ) 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. //- Check mesh topology for correctness.
// Returns false for no error. // Returns false for no error.

116
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
View File

@ -37,6 +37,7 @@ Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000; Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg
Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4; Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4;
Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * 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 * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::primitiveMesh::checkTopology(const bool report) const bool Foam::primitiveMesh::checkTopology(const bool report) const

6
src/OpenFOAM/primitives/Lists/stringListOpsTemplates.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -57,8 +57,12 @@ StringListType Foam::subsetMatchingStrings
const bool invert const bool invert
) )
{ {
// Create copy
StringListType newLst(lst.size()); StringListType newLst(lst.size());
// ensure consistent addressable size (eg, DynamicList)
newLst.setSize(lst.size());
label nElem = 0; label nElem = 0;
forAll(lst, elemI) forAll(lst, elemI)
{ {

35
src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
View File

@ -202,12 +202,45 @@ Foam::label Foam::removePoints::countPointUsage
pointCanBeDeleted[pointI] = true; pointCanBeDeleted[pointI] = true;
nDeleted++; nDeleted++;
} }
} }
edge0.clear(); edge0.clear();
edge1.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 // Point can be deleted only if all processors want to delete it
syncTools::syncPointList syncTools::syncPointList
( (

51
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View File

@ -1375,7 +1375,7 @@ void Foam::autoLayerDriver::growNoExtrusion
pointField& patchDisp, pointField& patchDisp,
labelList& patchNLayers, labelList& patchNLayers,
List<extrudeMode>& extrudeStatus List<extrudeMode>& extrudeStatus
) ) const
{ {
Info<< nl << "Growing non-extrusion points by one layer ..." << endl; Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
@ -1406,7 +1406,7 @@ void Foam::autoLayerDriver::growNoExtrusion
{ {
if if
( (
extrudeStatus[f[fp]] == NOEXTRUDE extrudeStatus[f[fp]] == EXTRUDE
&& grownExtrudeStatus[f[fp]] != NOEXTRUDE && grownExtrudeStatus[f[fp]] != NOEXTRUDE
) )
{ {
@ -1419,6 +1419,31 @@ void Foam::autoLayerDriver::growNoExtrusion
extrudeStatus.transfer(grownExtrudeStatus); 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) forAll(extrudeStatus, patchPointI)
{ {
if (extrudeStatus[patchPointI] == NOEXTRUDE) if (extrudeStatus[patchPointI] == NOEXTRUDE)
@ -2711,8 +2736,6 @@ void Foam::autoLayerDriver::addLayers
const labelList& meshPoints = patches[patchI].meshPoints(); const labelList& meshPoints = patches[patchI].meshPoints();
//scalar maxThickness = -VGREAT;
//scalar minThickness = VGREAT;
scalar sumThickness = 0; scalar sumThickness = 0;
scalar sumNearWallThickness = 0; scalar sumNearWallThickness = 0;
@ -2720,8 +2743,6 @@ void Foam::autoLayerDriver::addLayers
{ {
label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]]; label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
//maxThickness = max(maxThickness, thickness[ppPointI]);
//minThickness = min(minThickness, thickness[ppPointI]);
sumThickness += thickness[ppPointI]; sumThickness += thickness[ppPointI];
label nLay = patchNLayers[ppPointI]; label nLay = patchNLayers[ppPointI];
@ -2749,8 +2770,6 @@ void Foam::autoLayerDriver::addLayers
if (totNPoints > 0) if (totNPoints > 0)
{ {
//reduce(maxThickness, maxOp<scalar>());
//reduce(minThickness, minOp<scalar>());
avgThickness = avgThickness =
returnReduce(sumThickness, sumOp<scalar>()) returnReduce(sumThickness, sumOp<scalar>())
/ totNPoints; / totNPoints;
@ -2766,8 +2785,6 @@ void Foam::autoLayerDriver::addLayers
<< " " << setw(6) << layerParams.numLayers()[patchI] << " " << setw(6) << layerParams.numLayers()[patchI]
<< " " << setw(8) << avgNearWallThickness << " " << setw(8) << avgNearWallThickness
<< " " << setw(8) << avgThickness << " " << setw(8) << avgThickness
//<< " " << setw(8) << minThickness
//<< " " << setw(8) << maxThickness
<< endl; << endl;
} }
Info<< endl; Info<< endl;
@ -3147,6 +3164,19 @@ void Foam::autoLayerDriver::addLayers
meshMover().movePoints(oldPoints); meshMover().movePoints(oldPoints);
meshMover().correct(); meshMover().correct();
// Grow out region of non-extrusion
for (label i = 0; i < layerParams.nGrow(); i++)
{
growNoExtrusion
(
pp,
patchDisp,
patchNLayers,
extrudeStatus
);
}
Info<< endl; Info<< endl;
} }
@ -3293,7 +3323,6 @@ void Foam::autoLayerDriver::doLayers
// Balance // Balance
if (Pstream::parRun())
if (Pstream::parRun() && preBalance) if (Pstream::parRun() && preBalance)
{ {
Info<< nl Info<< nl

4
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
View File

@ -232,13 +232,13 @@ class autoLayerDriver
) const; ) const;
//- Grow no-extrusion layer. //- Grow no-extrusion layer.
static void growNoExtrusion void growNoExtrusion
( (
const indirectPrimitivePatch& pp, const indirectPrimitivePatch& pp,
pointField& patchDisp, pointField& patchDisp,
labelList& patchNLayers, labelList& patchNLayers,
List<extrudeMode>& extrudeStatus List<extrudeMode>& extrudeStatus
); ) const;
//- Calculate pointwise wanted and minimum thickness. //- Calculate pointwise wanted and minimum thickness.
// thickness: wanted thickness // thickness: wanted thickness

79
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View File

@ -245,27 +245,8 @@ void Foam::meshRefinement::getBafflePatches
const pointField& cellCentres = mesh_.cellCentres(); const pointField& cellCentres = mesh_.cellCentres();
// Build list of surfaces that are not to be baffled. // Surfaces that need to be baffled
const wordList& cellZoneNames = surfaces_.cellZoneNames(); const labelList surfacesToBaffle(surfaces_.getUnnamedSurfaces());
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);
ownPatch.setSize(mesh_.nFaces()); ownPatch.setSize(mesh_.nFaces());
ownPatch = -1; 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()); labelList namedSurfaces(surfaces_.getNamedSurfaces());
boolList isNamedSurface(cellZoneNames.size(), false);
forAll(namedSurfaces, i) forAll(namedSurfaces, i)
{ {
label surfI = namedSurfaces[i]; label surfI = namedSurfaces[i];
isNamedSurface[surfI] = true;
Info<< "Surface : " << surfaces_.names()[surfI] << nl Info<< "Surface : " << surfaces_.names()[surfI] << nl
<< " faceZone : " << faceZoneNames[surfI] << nl << " faceZone : " << faceZoneNames[surfI] << nl
<< " cellZone : " << cellZoneNames[surfI] << endl; << " cellZone : " << cellZoneNames[surfI] << endl;
@ -2125,6 +2111,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
{ {
label surfI = namedSurfaces[i]; label surfI = namedSurfaces[i];
if (cellZoneNames[surfI] != word::null)
{
label zoneI = cellZones.findZoneID(cellZoneNames[surfI]); label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
if (zoneI == -1) if (zoneI == -1)
@ -2151,6 +2139,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
} }
surfaceToCellZone[surfI] = zoneI; surfaceToCellZone[surfI] = zoneI;
} }
}
// Check they are synced // Check they are synced
List<wordList> allCellZones(Pstream::nProcs()); List<wordList> allCellZones(Pstream::nProcs());
@ -2321,6 +2310,11 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
if (closedNamedSurfaces.size()) if (closedNamedSurfaces.size())
{ {
Info<< "Found " << closedNamedSurfaces.size()
<< " closed, named surfaces. Assigning cells in/outside"
<< " these surfaces to the corresponding cellZone."
<< nl << endl;
findCellZoneGeometric findCellZoneGeometric
( (
closedNamedSurfaces, // indices of closed surfaces closedNamedSurfaces, // indices of closed surfaces
@ -2333,8 +2327,12 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Set using walking // 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 // Topological walk
findCellZoneTopo findCellZoneTopo
( (
@ -2349,6 +2347,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Make sure namedSurfaceIndex is unset inbetween same cell cell zones. // Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
if (!allowFreeStandingZoneFaces) if (!allowFreeStandingZoneFaces)
{ {
Info<< "Only keeping zone faces inbetween different cellZones."
<< nl << endl;
makeConsistentFaceIndex(cellToZone, namedSurfaceIndex); makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
} }
@ -2521,6 +2522,32 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
mesh_.setInstance(oldInstance()); 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... // None of the faces has changed, only the zones. Still...
updateMesh(map, labelList()); updateMesh(map, labelList());

72
src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
View File

@ -77,8 +77,29 @@ Foam::refinementSurfaces::refinementSurfaces
if (dict.found("faceZone")) if (dict.found("faceZone"))
{ {
dict.lookup("faceZone") >> faceZoneNames_[surfI]; dict.lookup("faceZone") >> faceZoneNames_[surfI];
dict.lookup("cellZone") >> cellZoneNames_[surfI]; bool hasSide = dict.readIfPresent("zoneInside", zoneInside_[surfI]);
dict.lookup("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 // Global perpendicular angle
@ -314,8 +335,40 @@ Foam::refinementSurfaces::refinementSurfaces
if (dict.found("faceZone")) if (dict.found("faceZone"))
{ {
dict.lookup("faceZone") >> faceZoneNames_[surfI]; dict.lookup("faceZone") >> faceZoneNames_[surfI];
dict.lookup("cellZone") >> cellZoneNames_[surfI]; bool hasSide = dict.readIfPresent
dict.lookup("zoneInside") >> zoneInside_[surfI]; (
"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 // Global perpendicular angle
@ -475,20 +528,19 @@ Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
// Get indices of closed named surfaces // Get indices of closed named surfaces
Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
{ {
labelList named(getNamedSurfaces()); labelList closed(cellZoneNames_.size());
labelList closed(named.size());
label closedI = 0; label closedI = 0;
forAll(cellZoneNames_, surfI)
forAll(named, i) {
if (cellZoneNames_[surfI].size())
{ {
label surfI = named[i];
if (allGeometry_[surfaces_[surfI]].hasVolumeType()) if (allGeometry_[surfaces_[surfI]].hasVolumeType())
{ {
closed[closedI++] = surfI; closed[closedI++] = surfI;
} }
} }
}
closed.setSize(closedI); closed.setSize(closedI);
return closed; return closed;

5
src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
View File

@ -146,7 +146,8 @@ public:
return faceZoneNames_; 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 const wordList& cellZoneNames() const
{ {
return cellZoneNames_; return cellZoneNames_;
@ -158,7 +159,7 @@ public:
//- Get indices of named surfaces (surfaces with faceZoneName) //- Get indices of named surfaces (surfaces with faceZoneName)
labelList getNamedSurfaces() const; labelList getNamedSurfaces() const;
//- Get indices of closed named surfaces //- Get indices of closed surfaces with a cellZone
labelList getClosedNamedSurfaces() const; labelList getClosedNamedSurfaces() const;
//- From local region number to global region number //- From local region number to global region number

4
src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C
View File

@ -637,7 +637,7 @@ Foam::directMappedPatchBase::directMappedPatchBase
samplePatch_(samplePatch), samplePatch_(samplePatch),
uniformOffset_(true), uniformOffset_(true),
offset_(offset), offset_(offset),
offsets_(0), offsets_(pp.size(), offset_),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()), sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mapPtr_(NULL) mapPtr_(NULL)
{} {}
@ -670,7 +670,7 @@ Foam::directMappedPatchBase::directMappedPatchBase
offsets_ offsets_
( (
uniformOffset_ uniformOffset_
? pointField(patch_.size(), offset_) ? pointField(pp.size(), offset_)
: dict.lookup("offsets") : dict.lookup("offsets")
), ),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()), sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),

2
src/turbulenceModels/incompressible/RAS/kOmega/kOmega.H
View File

@ -33,7 +33,7 @@ Description
"Turbulence Modeling for CFD" "Turbulence Modeling for CFD"
D. C. Wilcox, D. C. Wilcox,
DCW Industries, Inc., La Canada, DCW Industries, Inc., La Canada,
California, 1998. California, 1988.
See also: See also:
http://www.cfd-online.com/Wiki/Wilcox's_k-omega_model http://www.cfd-online.com/Wiki/Wilcox's_k-omega_model

1
tutorials/incompressible/MRFSimpleFoam/MRFSimpleFoam/MRFSimpleFoam.C
View File

@ -64,7 +64,6 @@ int main(int argc, char *argv[])
tmp<fvVectorMatrix> UEqn tmp<fvVectorMatrix> UEqn
( (
fvm::div(phi, U) fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevReff(U) + turbulence->divDevReff(U)
); );
mrfZones.addCoriolis(UEqn()); mrfZones.addCoriolis(UEqn());