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 /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2011-11-02 12:46:31 +00:00
parent 2f980abf24 8f1404db7f
commit b1c932c144
80 changed files with 2187 additions and 727 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/combustion/XiFoam/bEqn.H
View File

@ -63,6 +63,7 @@ if (ign.ignited())
// Solve for b
// ~~~~~~~~~~~
bEqn.relax();
bEqn.solve();
Info<< "min(b) = " << min(b).value() << endl;
@ -201,7 +202,7 @@ if (ign.ignited())
volScalarField XiEq
(
scalar(1.001)
+ (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
+ (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - min(max(b, 0.0), 1.0)))
*(XiEqStar - scalar(1.001))
);

190
applications/solvers/stressAnalysis/solidDisplacementFoam/readMechanicalProperties.H
View File

@ -13,70 +13,176 @@
);
const dictionary& rhoDict(mechanicalProperties.subDict("rho"));
word rhoType(rhoDict.lookup("rho"));
word rhoType(rhoDict.lookup("type"));
volScalarField rho
autoPtr<volScalarField> rhoPtr;
IOobject rhoIO
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
"rho",
runTime.timeName(0),
mesh,
dimensionedScalar("zero", dimMass/dimVolume, 0.0)
IOobject::NO_READ,
IOobject::NO_WRITE
);
if (rhoType == "rhoInf")
if (rhoType == "uniform")
{
rho = rhoDict.lookup("rhoInf");
scalar rhoValue(readScalar(rhoDict.lookup("value")));
rhoPtr.reset
(
new volScalarField
(
rhoIO,
mesh,
dimensionedScalar
(
"rho",
dimMass/dimVolume,
rhoValue
),
zeroGradientFvPatchField<scalar>::typeName
)
);
}
else if (rhoType == "field")
{
rhoIO.readOpt() = IOobject::MUST_READ;
rhoPtr.reset
(
new volScalarField
(
rhoIO,
mesh
)
);
}
else
{
FatalErrorIn
(
"readMechanicalProperties.H"
) << "Valid type entries are uniform or field for rho"
<< abort(FatalError);
}
volScalarField rhoE
(
IOobject
(
"E",
runTime.timeName(0),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("0", dimMass/dimLength/sqr(dimTime), 0.0)
);
volScalarField& rho = rhoPtr();
const dictionary& EDict(mechanicalProperties.subDict("E"));
word EType(EDict.lookup("E"));
if (EType == "EInf")
word EType(EDict.lookup("type"));
autoPtr<volScalarField> EPtr;
IOobject EIO
(
"E",
runTime.timeName(0),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
);
if (EType == "uniform")
{
rhoE = EDict.lookup("EInf");
scalar rhoEValue(readScalar(EDict.lookup("value")));
EPtr.reset
(
new volScalarField
(
EIO,
mesh,
dimensionedScalar
(
"Erho",
dimMass/dimLength/sqr(dimTime),
rhoEValue
),
zeroGradientFvPatchField<scalar>::typeName
)
);
}
else if (EType == "field")
{
EIO.readOpt() = IOobject::MUST_READ;
EPtr.reset
(
new volScalarField
(
EIO,
mesh
)
);
}
else
{
FatalErrorIn
(
"readMechanicalProperties.H"
) << "Valid type entries are uniform or field for E"
<< abort(FatalError);
}
volScalarField& rhoE = EPtr();
volScalarField nu
autoPtr<volScalarField> nuPtr;
IOobject nuIO
(
IOobject
(
"nu",
runTime.timeName(0),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
"nu",
runTime.timeName(0),
mesh,
dimensionedScalar("0", dimless, 0.0)
IOobject::NO_READ,
IOobject::NO_WRITE
);
const dictionary& nuDict(mechanicalProperties.subDict("nu"));
word nuType(nuDict.lookup("nu"));
word nuType(nuDict.lookup("type"));
if (nuType == "nuInf")
if (nuType == "uniform")
{
nu = nuDict.lookup("nuInf");
scalar nuValue(readScalar(nuDict.lookup("value")));
nuPtr.reset
(
new volScalarField
(
nuIO,
mesh,
dimensionedScalar
(
"nu",
dimless,
nuValue
),
zeroGradientFvPatchField<scalar>::typeName
)
);
}
else if(nuType == "field")
{
nuIO.readOpt() = IOobject::MUST_READ;
nuPtr.reset
(
new volScalarField
(
nuIO,
mesh
)
);
}
else
{
FatalErrorIn
(
"readMechanicalProperties.H"
) << "Valid type entries are uniform or field for nu"
<< abort(FatalError);
}
volScalarField& nu = nuPtr();
Info<< "Normalising E : E/rho\n" << endl;
volScalarField E = rhoE/rho;

190
applications/solvers/stressAnalysis/solidDisplacementFoam/readThermalProperties.H
View File

@ -46,70 +46,180 @@ volScalarField DT
if (thermalStress)
{
volScalarField C
autoPtr<volScalarField> CPtr;
IOobject CIO
(
IOobject
(
"C",
runTime.timeName(0),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
"C",
runTime.timeName(0),
mesh,
dimensionedScalar("0", dimensionSet(0, 2, -2 , -1, 0), 0.0)
IOobject::NO_READ,
IOobject::NO_WRITE
);
const dictionary& CDict(thermalProperties.subDict("C"));
word CType(CDict.lookup("C"));
if (CType == "CInf")
word CType(CDict.lookup("type"));
if (CType == "uniform")
{
C = CDict.lookup("CInf");
scalar CValue(readScalar(CDict.lookup("value")));
CPtr.reset
(
new volScalarField
(
CIO,
mesh,
dimensionedScalar
(
"C",
dimensionSet(0, 2, -2 , -1, 0),
CValue
),
zeroGradientFvPatchField<scalar>::typeName
)
);
}
else if(CType == "field")
{
CIO.readOpt() = IOobject::MUST_READ;
CPtr.reset
(
new volScalarField
(
CIO,
mesh
)
);
}
else
{
FatalErrorIn
(
"readThermalProperties.H"
) << "Valid type entries are uniform or field for C"
<< abort(FatalError);
}
volScalarField& C = CPtr();
volScalarField rhoK
autoPtr<volScalarField> rhoKPtr;
IOobject rhoKIO
(
IOobject
(
"k",
runTime.timeName(0),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
"k",
runTime.timeName(0),
mesh,
dimensionedScalar("0", dimensionSet(1, 1, -3 , -1, 0), 0.0)
IOobject::NO_READ,
IOobject::NO_WRITE
);
const dictionary& kDict(thermalProperties.subDict("k"));
word kType(kDict.lookup("k"));
if (kType == "kInf")
word kType(kDict.lookup("type"));
if (kType == "uniform")
{
rhoK = kDict.lookup("kInf");
scalar rhoKValue(readScalar(kDict.lookup("value")));
rhoKPtr.reset
(
new volScalarField
(
rhoKIO,
mesh,
dimensionedScalar
(
"rhoK",
dimensionSet(1, 1, -3 , -1, 0),
rhoKValue
),
zeroGradientFvPatchField<scalar>::typeName
)
);
}
else if (kType == "field")
{
rhoKIO.readOpt() = IOobject::MUST_READ;
rhoKPtr.reset
(
new volScalarField
(
rhoKIO,
mesh
)
);
}
else
{
FatalErrorIn
(
"readThermalProperties.H"
) << "Valid type entries are uniform or field for K"
<< abort(FatalError);
}
volScalarField alpha
volScalarField& rhoK = rhoKPtr();
autoPtr<volScalarField> alphaPtr;
IOobject alphaIO
(
IOobject
(
"alpha",
runTime.timeName(0),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
"alpha",
runTime.timeName(0),
mesh,
dimensionedScalar("0", dimensionSet(0, 0, 0 , -1, 0), 0.0)
IOobject::NO_READ,
IOobject::NO_WRITE
);
const dictionary& alphaDict(thermalProperties.subDict("alpha"));
word alphaType(alphaDict.lookup("alpha"));
if (alphaType == "alphaInf")
const dictionary& alphaDict(thermalProperties.subDict("alpha"));
word alphaType(alphaDict.lookup("type"));
if (alphaType == "uniform")
{
alpha = alphaDict.lookup("alphaInf");
scalar alphaValue(readScalar(alphaDict.lookup("value")));
alphaPtr.reset
(
new volScalarField
(
alphaIO,
mesh,
dimensionedScalar
(
"alpha",
inv(dimTemperature),
alphaValue
),
zeroGradientFvPatchField<scalar>::typeName
)
);
}
else if (alphaType == "field")
{
alphaIO.readOpt() = IOobject::MUST_READ;
alphaPtr.reset
(
new volScalarField
(
alphaIO,
mesh
)
);
}
else
{
FatalErrorIn
(
"readThermalProperties.H"
) << "Valid type entries are uniform or field for alpha"
<< abort(FatalError);
}
volScalarField& alpha = alphaPtr();
Info<< "Normalising k : k/rho\n" << endl;
volScalarField k = rhoK/rho;

21
applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/tractionDisplacementCorrection/tractionDisplacementCorrectionFvPatchVectorField.C
View File

@ -150,13 +150,18 @@ void tractionDisplacementCorrectionFvPatchVectorField::updateCoeffs()
"mechanicalProperties"
);
dimensionedScalar rho(mechanicalProperties.lookup("rho"));
dimensionedScalar rhoE(mechanicalProperties.lookup("E"));
dimensionedScalar nu(mechanicalProperties.lookup("nu"));
const fvPatchField<scalar>& rho =
patch().lookupPatchField<volScalarField, scalar>("rho");
dimensionedScalar E = rhoE/rho;
dimensionedScalar mu = E/(2.0*(1.0 + nu));
dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
const fvPatchField<scalar>& rhoE =
patch().lookupPatchField<volScalarField, scalar>("E");
const fvPatchField<scalar>& nu =
patch().lookupPatchField<volScalarField, scalar>("nu");
scalarField E = rhoE/rho;
scalarField mu = E/(2.0*(1.0 + nu));
scalarField lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
Switch planeStress(mechanicalProperties.lookup("planeStress"));
@ -175,8 +180,8 @@ void tractionDisplacementCorrectionFvPatchVectorField::updateCoeffs()
gradient() =
(
(traction_ + pressure_*n)/rho.value() - (n & (sigmaD + sigmaExp))
)/(2.0*mu + lambda).value();
(traction_ + pressure_*n)/rho - (n & (sigmaD + sigmaExp))
)/(2.0*mu + lambda);
fixedGradientFvPatchVectorField::updateCoeffs();
}

2
applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C
View File

@ -715,7 +715,7 @@ int main(int argc, char *argv[])
triSurfaceSearch querySurf(surf);
// Search engine on mesh. No face decomposition since mesh unwarped.
meshSearch queryMesh(mesh, false);
meshSearch queryMesh(mesh, polyMesh::FACEPLANES);
// Check all 'outside' points
forAll(outsidePts, outsideI)

2
applications/utilities/mesh/advanced/selectCells/selectCells.C
View File

@ -381,7 +381,7 @@ int main(int argc, char *argv[])
(void)edgeCalc.minLen(Info);
// Search engine on mesh. Face decomposition since faces might be warped.
meshSearch queryMesh(mesh, true);
meshSearch queryMesh(mesh, polyMesh::FACEDIAGTETS);
// Check all 'outside' points
forAll(outsidePts, outsideI)

2
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
View File

@ -971,7 +971,7 @@ Foam::backgroundMeshDecomposition::distribute
{
// Map cellVertices, cellVertexIndices and cellVertexTypes
meshSearch cellSearch(mesh_);
meshSearch cellSearch(mesh_, polyMesh::FACEPLANES);
const labelList& reverseCellMap = map().reverseCellMap();

2
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

@ -1405,7 +1405,7 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
zeroGradientFvPatchScalarField::typeName
);
meshSearch cellSearch(bMesh);
meshSearch cellSearch(bMesh, polyMesh::FACEPLANES);
List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
List<DynamicList<label> > cellVertexIndices(bMesh.nCells());

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

@ -2194,7 +2194,7 @@ int main(int argc, char *argv[])
label regionI = -1;
label cellI = mesh.findCell(insidePoint);
label cellI = mesh.findCell(insidePoint, polyMesh::FACEDIAGTETS);
Info<< nl << "Found point " << insidePoint << " in cell " << cellI
<< endl;

6
applications/utilities/parallelProcessing/decomposePar/decomposePar.C
View File

@ -528,7 +528,11 @@ int main(int argc, char *argv[])
<< "Cell number should be between 0 and "
<< mesh.nCells()-1 << nl
<< "On this mesh the particle should be in cell "
<< mesh.findCell(iter().position())
<< mesh.findCell
(
iter().position(),
polyMesh::FACEDIAGTETS
)
<< exit(FatalError);
}

6
applications/utilities/preProcessing/mapFields/mapLagrangian.C
View File

@ -56,7 +56,11 @@ static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
// See if particle on face by finding nearest face and shifting
// particle.
meshSearch meshSearcher(mesh, false);
meshSearch meshSearcher
(
mesh,
polyMesh::FACEPLANES // no decomposition needed
);
label faceI = meshSearcher.findNearestBoundaryFace(pt);

26
src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C
View File

@ -25,7 +25,7 @@ License
#include "treeDataCell.H"
#include "indexedOctree.H"
#include "primitiveMesh.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -83,13 +83,15 @@ void Foam::treeDataCell::update()
Foam::treeDataCell::treeDataCell
(
const bool cacheBb,
const primitiveMesh& mesh,
const labelUList& cellLabels
const polyMesh& mesh,
const labelUList& cellLabels,
const polyMesh::cellRepresentation decompMode
)
:
mesh_(mesh),
cellLabels_(cellLabels),
cacheBb_(cacheBb)
cacheBb_(cacheBb),
decompMode_(decompMode)
{
update();
}
@ -98,13 +100,15 @@ Foam::treeDataCell::treeDataCell
Foam::treeDataCell::treeDataCell
(
const bool cacheBb,
const primitiveMesh& mesh,
const Xfer<labelList>& cellLabels
const polyMesh& mesh,
const Xfer<labelList>& cellLabels,
const polyMesh::cellRepresentation decompMode
)
:
mesh_(mesh),
cellLabels_(cellLabels),
cacheBb_(cacheBb)
cacheBb_(cacheBb),
decompMode_(decompMode)
{
update();
}
@ -113,12 +117,14 @@ Foam::treeDataCell::treeDataCell
Foam::treeDataCell::treeDataCell
(
const bool cacheBb,
const primitiveMesh& mesh
const polyMesh& mesh,
const polyMesh::cellRepresentation decompMode
)
:
mesh_(mesh),
cellLabels_(identity(mesh_.nCells())),
cacheBb_(cacheBb)
cacheBb_(cacheBb),
decompMode_(decompMode)
{
update();
}
@ -162,7 +168,7 @@ bool Foam::treeDataCell::contains
const point& sample
) const
{
return mesh_.pointInCell(sample, cellLabels_[index]);
return mesh_.pointInCell(sample, cellLabels_[index], decompMode_);
}

30
src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H
View File

@ -36,6 +36,7 @@ SourceFiles
#ifndef treeDataCell_H
#define treeDataCell_H
#include "polyMesh.H"
#include "treeBoundBoxList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,7 +45,6 @@ namespace Foam
{
// Forward declaration of classes
class primitiveMesh;
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
@ -55,7 +55,7 @@ class treeDataCell
{
// Private data
const primitiveMesh& mesh_;
const polyMesh& mesh_;
//- Subset of cells to work on
const labelList cellLabels_;
@ -63,6 +63,9 @@ class treeDataCell
//- Whether to precalculate and store cell bounding box
const bool cacheBb_;
//- How to decide if point is inside cell
const polyMesh::cellRepresentation decompMode_;
//- cell bounding boxes (valid only if cacheBb_)
treeBoundBoxList bbs_;
@ -87,20 +90,27 @@ public:
treeDataCell
(
const bool cacheBb,
const primitiveMesh&,
const labelUList&
const polyMesh&,
const labelUList&,
const polyMesh::cellRepresentation decompMode
);
//- Construct from mesh and subset of cells, transferring contents
treeDataCell
(
const bool cacheBb,
const primitiveMesh&,
const Xfer<labelList>&
const polyMesh&,
const Xfer<labelList>&,
const polyMesh::cellRepresentation decompMode
);
//- Construct from mesh. Uses all cells in mesh.
treeDataCell(const bool cacheBb, const primitiveMesh&);
treeDataCell
(
const bool cacheBb,
const polyMesh&,
const polyMesh::cellRepresentation decompMode
);
// Member Functions
@ -112,11 +122,15 @@ public:
return cellLabels_;
}
inline const primitiveMesh& mesh() const
inline const polyMesh& mesh() const
{
return mesh_;
}
inline polyMesh::cellRepresentation decompMode() const
{
return decompMode_;
}
inline label size() const
{

180
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C
View File

@ -25,7 +25,8 @@ License
#include "cyclicGAMGInterface.H"
#include "addToRunTimeSelectionTable.H"
#include "Map.H"
#include "labelPair.H"
#include "HashTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -62,180 +63,65 @@ Foam::cyclicGAMGInterface::cyclicGAMGInterface
),
fineCyclicInterface_(refCast<const cyclicLduInterface>(fineInterface))
{
// Make a lookup table of entries for owner/neighbour
Map<SLList<label> > neighboursTable
// From coarse face to coarse cell
DynamicList<label> dynFaceCells(localRestrictAddressing.size());
// From fine face to coarse face
DynamicList<label> dynFaceRestrictAddressing
(
localRestrictAddressing.size()
);
// Table of face-sets to be agglomerated
Map<SLList<SLList<label> > > faceFaceTable
// From coarse cell pair to coarse face
HashTable<label, labelPair, labelPair::Hash<> > cellsToCoarseFace
(
localRestrictAddressing.size()
2*localRestrictAddressing.size()
);
label nCoarseFaces = 0;
forAll(localRestrictAddressing, ffi)
{
label curMaster = -1;
label curSlave = -1;
labelPair cellPair;
// Do switching on master/slave indexes based on the owner/neighbour of
// the processor index such that both sides get the same answer.
if (owner())
{
// Master side
curMaster = localRestrictAddressing[ffi];
curSlave = neighbourRestrictAddressing[ffi];
cellPair = labelPair
(
localRestrictAddressing[ffi],
neighbourRestrictAddressing[ffi]
);
}
else
{
// Slave side
curMaster = neighbourRestrictAddressing[ffi];
curSlave = localRestrictAddressing[ffi];
cellPair = labelPair
(
neighbourRestrictAddressing[ffi],
localRestrictAddressing[ffi]
);
}
// Look for the master cell. If it has already got a face,
// add the coefficient to the face. If not, create a new face.
if (neighboursTable.found(curMaster))
HashTable<label, labelPair, labelPair::Hash<> >::const_iterator fnd =
cellsToCoarseFace.find(cellPair);
if (fnd == cellsToCoarseFace.end())
{
// Check all current neighbours to see if the current slave already
// exists and if so, add the fine face to the agglomeration.
SLList<label>& curNbrs = neighboursTable.find(curMaster)();
SLList<SLList<label> >& curFaceFaces =
faceFaceTable.find(curMaster)();
bool nbrFound = false;
SLList<label>::iterator nbrsIter = curNbrs.begin();
SLList<SLList<label> >::iterator faceFacesIter =
curFaceFaces.begin();
for
(
;
nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
++nbrsIter, ++faceFacesIter
)
{
if (nbrsIter() == curSlave)
{
nbrFound = true;
faceFacesIter().append(ffi);
break;
}
}
if (!nbrFound)
{
curNbrs.append(curSlave);
curFaceFaces.append(ffi);
// New coarse face created
nCoarseFaces++;
}
// New coarse face
label coarseI = dynFaceCells.size();
dynFaceRestrictAddressing.append(coarseI);
dynFaceCells.append(localRestrictAddressing[ffi]);
cellsToCoarseFace.insert(cellPair, coarseI);
}
else
{
// This master has got no neighbours yet. Add a neighbour
// and a coefficient, thus creating a new face
neighboursTable.insert(curMaster, SLList<label>(curSlave));
faceFaceTable.insert(curMaster, SLList<SLList<label> >(ffi));
// New coarse face created
nCoarseFaces++;
}
} // end for all fine faces
faceCells_.setSize(nCoarseFaces, -1);
faceRestrictAddressing_.setSize(localRestrictAddressing.size());
labelList contents = neighboursTable.toc();
// Reset face counter for re-use
nCoarseFaces = 0;
if (owner())
{
// On master side, the owner addressing is stored in table of contents
forAll(contents, masterI)
{
SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
SLList<SLList<label> >& curFaceFaces =
faceFaceTable.find(contents[masterI])();
SLList<label>::iterator nbrsIter = curNbrs.begin();
SLList<SLList<label> >::iterator faceFacesIter =
curFaceFaces.begin();
for
(
;
nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
++nbrsIter, ++faceFacesIter
)
{
faceCells_[nCoarseFaces] = contents[masterI];
for
(
SLList<label>::iterator facesIter = faceFacesIter().begin();
facesIter != faceFacesIter().end();
++facesIter
)
{
faceRestrictAddressing_[facesIter()] = nCoarseFaces;
}
nCoarseFaces++;
}
// Already have coarse face
dynFaceRestrictAddressing.append(fnd());
}
}
else
{
// On slave side, the owner addressing is stored in linked lists
forAll(contents, masterI)
{
SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
SLList<SLList<label> >& curFaceFaces =
faceFaceTable.find(contents[masterI])();
SLList<label>::iterator nbrsIter = curNbrs.begin();
SLList<SLList<label> >::iterator faceFacesIter =
curFaceFaces.begin();
for
(
;
nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
++nbrsIter, ++faceFacesIter
)
{
faceCells_[nCoarseFaces] = nbrsIter();
for
(
SLList<label>::iterator facesIter = faceFacesIter().begin();
facesIter != faceFacesIter().end();
++facesIter
)
{
faceRestrictAddressing_[facesIter()] = nCoarseFaces;
}
nCoarseFaces++;
}
}
}
faceCells_.transfer(dynFaceCells);
faceRestrictAddressing_.transfer(dynFaceRestrictAddressing);
}

171
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/processorGAMGInterface/processorGAMGInterface.C
View File

@ -25,7 +25,8 @@ License
#include "processorGAMGInterface.H"
#include "addToRunTimeSelectionTable.H"
#include "Map.H"
#include "HashTable.H"
#include "labelPair.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -62,171 +63,65 @@ Foam::processorGAMGInterface::processorGAMGInterface
),
fineProcInterface_(refCast<const processorLduInterface>(fineInterface))
{
// Make a lookup table of entries for owner/neighbour
Map<SLList<label> > neighboursTable
// From coarse face to coarse cell
DynamicList<label> dynFaceCells(localRestrictAddressing.size());
// From fine face to coarse face
DynamicList<label> dynFaceRestrictAddressing
(
localRestrictAddressing.size()
);
// Table of face-sets to be agglomerated
Map<SLList<SLList<label> > > faceFaceTable
// From coarse cell pair to coarse face
HashTable<label, labelPair, labelPair::Hash<> > cellsToCoarseFace
(
localRestrictAddressing.size()
2*localRestrictAddressing.size()
);
label nCoarseFaces = 0;
forAll(localRestrictAddressing, ffi)
{
label curMaster = -1;
label curSlave = -1;
labelPair cellPair;
// Do switching on master/slave indexes based on the owner/neighbour of
// the processor index such that both sides get the same answer.
if (myProcNo() < neighbProcNo())
{
// Master side
curMaster = localRestrictAddressing[ffi];
curSlave = neighbourRestrictAddressing[ffi];
cellPair = labelPair
(
localRestrictAddressing[ffi],
neighbourRestrictAddressing[ffi]
);
}
else
{
// Slave side
curMaster = neighbourRestrictAddressing[ffi];
curSlave = localRestrictAddressing[ffi];
cellPair = labelPair
(
neighbourRestrictAddressing[ffi],
localRestrictAddressing[ffi]
);
}
// Look for the master cell. If it has already got a face,
// add the coefficient to the face. If not, create a new face.
if (neighboursTable.found(curMaster))
HashTable<label, labelPair, labelPair::Hash<> >::const_iterator fnd =
cellsToCoarseFace.find(cellPair);
if (fnd == cellsToCoarseFace.end())
{
// Check all current neighbours to see if the current slave already
// exists and if so, add the fine face to the agglomeration.
SLList<label>& curNbrs = neighboursTable.find(curMaster)();
SLList<SLList<label> >& curFaceFaces =
faceFaceTable.find(curMaster)();
bool nbrFound = false;
SLList<label>::iterator nbrsIter = curNbrs.begin();
SLList<SLList<label> >::iterator faceFacesIter =
curFaceFaces.begin();
for
(
;
nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
++nbrsIter, ++faceFacesIter
)
{
if (nbrsIter() == curSlave)
{
nbrFound = true;
faceFacesIter().append(ffi);
break;
}
}
if (!nbrFound)
{
curNbrs.append(curSlave);
curFaceFaces.append(ffi);
// New coarse face created
nCoarseFaces++;
}
// New coarse face
label coarseI = dynFaceCells.size();
dynFaceRestrictAddressing.append(coarseI);
dynFaceCells.append(localRestrictAddressing[ffi]);
cellsToCoarseFace.insert(cellPair, coarseI);
}
else
{
// This master has got no neighbours yet. Add a neighbour
// and a coefficient, thus creating a new face
neighboursTable.insert(curMaster, SLList<label>(curSlave));
faceFaceTable.insert(curMaster, SLList<SLList<label> >(ffi));
// New coarse face created
nCoarseFaces++;
}
} // end for all fine faces
faceCells_.setSize(nCoarseFaces, -1);
faceRestrictAddressing_.setSize(localRestrictAddressing.size());
labelList contents = neighboursTable.toc();
// Reset face counter for re-use
nCoarseFaces = 0;
if (myProcNo() < neighbProcNo())
{
// On master side, the owner addressing is stored in table of contents
forAll(contents, masterI)
{
SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
SLList<SLList<label> >& curFaceFaces =
faceFaceTable.find(contents[masterI])();
SLList<label>::iterator nbrsIter = curNbrs.begin();
SLList<SLList<label> >::iterator faceFacesIter =
curFaceFaces.begin();
for
(
;
nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
++nbrsIter, ++faceFacesIter
)
{
faceCells_[nCoarseFaces] = contents[masterI];
forAllConstIter(SLList<label>, faceFacesIter(), facesIter)
{
faceRestrictAddressing_[facesIter()] = nCoarseFaces;
}
nCoarseFaces++;
}
// Already have coarse face
dynFaceRestrictAddressing.append(fnd());
}
}
else
{
// On slave side, the owner addressing is stored in linked lists
forAll(contents, masterI)
{
SLList<label>& curNbrs = neighboursTable.find(contents[masterI])();
SLList<SLList<label> >& curFaceFaces =
faceFaceTable.find(contents[masterI])();
SLList<label>::iterator nbrsIter = curNbrs.begin();
SLList<SLList<label> >::iterator faceFacesIter =
curFaceFaces.begin();
for
(
;
nbrsIter != curNbrs.end(), faceFacesIter != curFaceFaces.end();
++nbrsIter, ++faceFacesIter
)
{
faceCells_[nCoarseFaces] = nbrsIter();
forAllConstIter(SLList<label>, faceFacesIter(), facesIter)
{
faceRestrictAddressing_[facesIter()] = nCoarseFaces;
}
nCoarseFaces++;
}
}
}
faceCells_.transfer(dynFaceCells);
faceRestrictAddressing_.transfer(dynFaceRestrictAddressing);
}

181
src/OpenFOAM/meshes/polyMesh/polyMesh.C
View File

@ -32,7 +32,6 @@ License
#include "globalMeshData.H"
#include "processorPolyPatch.H"
#include "OSspecific.H"
#include "demandDrivenData.H"
#include "polyMeshTetDecomposition.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
@ -209,6 +208,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
cellTreePtr_(NULL),
pointZones_
(
IOobject
@ -401,6 +401,7 @@ Foam::polyMesh::polyMesh
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
cellTreePtr_(NULL),
pointZones_
(
IOobject
@ -558,6 +559,7 @@ Foam::polyMesh::polyMesh
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
cellTreePtr_(NULL),
pointZones_
(
IOobject
@ -859,7 +861,7 @@ Foam::label Foam::polyMesh::nSolutionD() const
const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
{
if (!tetBasePtIsPtr_)
if (tetBasePtIsPtr_.empty())
{
if (debug)
{
@ -869,13 +871,51 @@ const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
<< endl;
}
tetBasePtIsPtr_ = new labelList
tetBasePtIsPtr_.reset
(
polyMeshTetDecomposition::findFaceBasePts(*this)
new labelList
(
polyMeshTetDecomposition::findFaceBasePts(*this)
)
);
}
return *tetBasePtIsPtr_;
return tetBasePtIsPtr_();
}
const Foam::indexedOctree<Foam::treeDataCell>&
Foam::polyMesh::cellTree() const
{
if (cellTreePtr_.empty())
{
treeBoundBox overallBb(points());
Random rndGen(261782);
overallBb = overallBb.extend(rndGen, 1E-4);
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
cellTreePtr_.reset
(
new indexedOctree<treeDataCell>
(
treeDataCell
(
false, // not cache bb
*this,
FACEDIAGTETS // use tetDecomposition for any inside test
),
overallBb,
8, // maxLevel
10, // leafsize
5.0 // duplicity
)
);
}
return cellTreePtr_();
}
@ -911,7 +951,7 @@ void Foam::polyMesh::addPatches
// recalculation. Problem: should really be done in removeBoundary but
// there is some info in parallelData which might be interesting inbetween
// removeBoundary and addPatches.
deleteDemandDrivenData(globalMeshDataPtr_);
globalMeshDataPtr_.clear();
if (validBoundary)
{
@ -1033,7 +1073,7 @@ const Foam::labelList& Foam::polyMesh::faceNeighbour() const
// Return old mesh motion points
const Foam::pointField& Foam::polyMesh::oldPoints() const
{
if (!oldPointsPtr_)
if (oldPointsPtr_.empty())
{
if (debug)
{
@ -1042,11 +1082,11 @@ const Foam::pointField& Foam::polyMesh::oldPoints() const
<< endl;
}
oldPointsPtr_ = new pointField(points_);
oldPointsPtr_.reset(new pointField(points_));
curMotionTimeIndex_ = time().timeIndex();
}
return *oldPointsPtr_;
return oldPointsPtr_();
}
@ -1068,8 +1108,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
if (curMotionTimeIndex_ != time().timeIndex())
{
// Mesh motion in the new time step
deleteDemandDrivenData(oldPointsPtr_);
oldPointsPtr_ = new pointField(points_);
oldPointsPtr_.clear();
oldPointsPtr_.reset(new pointField(points_));
curMotionTimeIndex_ = time().timeIndex();
}
@ -1099,9 +1139,9 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
);
// Adjust parallel shared points
if (globalMeshDataPtr_)
if (globalMeshDataPtr_.valid())
{
globalMeshDataPtr_->movePoints(points_);
globalMeshDataPtr_().movePoints(points_);
}
// Force recalculation of all geometric data with new points
@ -1141,14 +1181,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
void Foam::polyMesh::resetMotion() const
{
curMotionTimeIndex_ = 0;
deleteDemandDrivenData(oldPointsPtr_);
oldPointsPtr_.clear();
}
// Return parallel info
const Foam::globalMeshData& Foam::polyMesh::globalData() const
{
if (!globalMeshDataPtr_)
if (globalMeshDataPtr_.empty())
{
if (debug)
{
@ -1158,10 +1198,10 @@ const Foam::globalMeshData& Foam::polyMesh::globalData() const
<< endl;
}
// Construct globalMeshData using processorPatch information only.
globalMeshDataPtr_ = new globalMeshData(*this);
globalMeshDataPtr_.reset(new globalMeshData(*this));
}
return *globalMeshDataPtr_;
return globalMeshDataPtr_();
}
@ -1308,4 +1348,111 @@ void Foam::polyMesh::findTetFacePt
}
bool Foam::polyMesh::pointInCell
(
const point& p,
label cellI,
const cellRepresentation decompMode
) const
{
switch (decompMode)
{
case FACEPLANES:
{
return primitiveMesh::pointInCell(p, cellI);
}
break;
case FACECENTRETETS:
{
const point& cc = cellCentres()[cellI];
const cell& cFaces = cells()[cellI];
forAll(cFaces, cFaceI)
{
label faceI = cFaces[cFaceI];
const face& f = faces_[faceI];
const point& fc = faceCentres()[faceI];
bool isOwn = (owner_[faceI] == cellI);
forAll(f, fp)
{
label pointI;
label nextPointI;
if (isOwn)
{
pointI = f[fp];
nextPointI = f.nextLabel(fp);
}
else
{
pointI = f.nextLabel(fp);
nextPointI = f[fp];
}
if
(
tetPointRef
(
points()[nextPointI],
points()[pointI],
fc,
cc
).inside(p)
)
{
return true;
}
}
}
return false;
}
break;
case FACEDIAGTETS:
{
label tetFaceI, tetPtI;
findTetFacePt(cellI, p, tetFaceI, tetPtI);
return tetFaceI != -1;
}
break;
}
return false;
}
Foam::label Foam::polyMesh::findCell
(
const point& location,
const cellRepresentation decompMode
) const
{
if (nCells() == 0)
{
return -1;
}
// Find the nearest cell centre to this location
label cellI = findNearestCell(location);
// If point is in the nearest cell return
if (pointInCell(location, cellI, decompMode))
{
return cellI;
}
else // point is not in the nearest cell so search all cells
{
for (label cellI = 0; cellI < nCells(); cellI++)
{
if (pointInCell(location, cellI, decompMode))
{
return cellI;
}
}
return -1;
}
}
// ************************************************************************* //

43
src/OpenFOAM/meshes/polyMesh/polyMesh.H
View File

@ -34,7 +34,6 @@ SourceFiles
polyMeshFromShapeMesh.C
polyMeshIO.C
polyMeshUpdate.C
polyMeshFindCell.C
\*---------------------------------------------------------------------------*/
@ -61,9 +60,12 @@ SourceFiles
namespace Foam
{
// Forward declaration of classes
class globalMeshData;
class mapPolyMesh;
class polyMeshTetDecomposition;
class treeDataCell;
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
Class polyMesh Declaration
@ -91,6 +93,15 @@ public:
TOPO_PATCH_CHANGE
};
//- Enumeration defining the representation of the cell for
// inside checking
enum cellRepresentation
{
FACEPLANES, // cell bound by planes of faces
FACECENTRETETS, // tet decomposition using facectr and cellctr
FACEDIAGTETS // tet decomposition using face diagonal and cellctr
};
private:
@ -130,7 +141,10 @@ private:
mutable Vector<label> solutionD_;
//- Base point for face decomposition into tets
mutable labelList* tetBasePtIsPtr_;
mutable autoPtr<labelList> tetBasePtIsPtr_;
//- Search tree to allow spatial cell searching
mutable autoPtr<indexedOctree<treeDataCell> > cellTreePtr_;
// Zoning information
@ -146,7 +160,7 @@ private:
//- Parallel info
mutable globalMeshData* globalMeshDataPtr_;
mutable autoPtr<globalMeshData> globalMeshDataPtr_;
// Mesh motion related data
@ -161,7 +175,7 @@ private:
mutable label curMotionTimeIndex_;
//- Old points (for the last mesh motion)
mutable pointField* oldPointsPtr_;
mutable autoPtr<pointField> oldPointsPtr_;
// Private Member Functions
@ -365,6 +379,9 @@ public:
//- Return the tetBasePtIs
const labelList& tetBasePtIs() const;
//- Return the cell search tree
const indexedOctree<treeDataCell>& cellTree() const;
//- Return point zone mesh
const pointZoneMesh& pointZones() const
{
@ -505,6 +522,9 @@ public:
//- Clear primitive data (points, faces and cells)
void clearPrimitives();
//- Clear cell tree data
void clearCellTree();
//- Remove all files from mesh instance
void removeFiles(const fileName& instanceDir) const;
@ -532,6 +552,21 @@ public:
label& tetFaceI,
label& tetPtI
) const;
//- Is the point in the cell
bool pointInCell
(
const point&,
label cellI,
const cellRepresentation
) const;
//- Find cell enclosing this location (-1 if not in mesh)
label findCell
(
const point&,
const cellRepresentation
) const;
};

15
src/OpenFOAM/meshes/polyMesh/polyMeshClear.C
View File

@ -26,9 +26,10 @@ License
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "globalMeshData.H"
#include "demandDrivenData.H"
#include "pointMesh.H"
#include "Time.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -72,7 +73,9 @@ void Foam::polyMesh::clearGeom()
solutionD_ = Vector<label>::zero;
// Remove the stored tet base points
deleteDemandDrivenData(tetBasePtIsPtr_);
tetBasePtIsPtr_.clear();
// Remove the cell tree
cellTreePtr_.clear();
pointMesh::Delete(*this);
}
@ -91,7 +94,7 @@ void Foam::polyMesh::clearAddressing()
// parallelData depends on the processorPatch ordering so force
// recalculation
deleteDemandDrivenData(globalMeshDataPtr_);
globalMeshDataPtr_.clear();
// Reset valid directions
geometricD_ = Vector<label>::zero;
@ -121,4 +124,10 @@ void Foam::polyMesh::clearOut()
}
void Foam::polyMesh::clearCellTree()
{
cellTreePtr_.clear();
}
// ************************************************************************* //

7
src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C
View File

@ -30,6 +30,9 @@ Description
#include "Time.H"
#include "primitiveMesh.H"
#include "DynamicList.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "globalMeshData.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -504,6 +507,7 @@ Foam::polyMesh::polyMesh
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
cellTreePtr_(NULL),
pointZones_
(
IOobject
@ -548,6 +552,7 @@ Foam::polyMesh::polyMesh
),
globalMeshDataPtr_(NULL),
moving_(false),
changing_(false),
curMotionTimeIndex_(time().timeIndex()),
oldPointsPtr_(NULL)
{
@ -778,6 +783,7 @@ Foam::polyMesh::polyMesh
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
tetBasePtIsPtr_(NULL),
cellTreePtr_(NULL),
pointZones_
(
IOobject
@ -822,6 +828,7 @@ Foam::polyMesh::polyMesh
),
globalMeshDataPtr_(NULL),
moving_(false),
changing_(false),
curMotionTimeIndex_(time().timeIndex()),
oldPointsPtr_(NULL)
{

8
src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
View File

@ -45,7 +45,7 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
cellZones_.clearAddressing();
// Update parallel data
if (globalMeshDataPtr_)
if (globalMeshDataPtr_.valid())
{
globalMeshDataPtr_->updateMesh();
}
@ -53,12 +53,12 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
setInstance(time().timeName());
// Map the old motion points if present
if (oldPointsPtr_)
if (oldPointsPtr_.valid())
{
// Make a copy of the original points
pointField oldMotionPoints = *oldPointsPtr_;
pointField oldMotionPoints = oldPointsPtr_();
pointField& newMotionPoints = *oldPointsPtr_;
pointField& newMotionPoints = oldPointsPtr_();
// Resize the list to new size
newMotionPoints.setSize(points_.size());

38
src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C
View File

@ -25,8 +25,6 @@ License
#include "primitiveMesh.H"
#include "demandDrivenData.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -63,8 +61,6 @@ Foam::primitiveMesh::primitiveMesh()
ppPtr_(NULL),
cpPtr_(NULL),
cellTreePtr_(NULL),
labels_(0),
cellCentresPtr_(NULL),
@ -107,8 +103,6 @@ Foam::primitiveMesh::primitiveMesh
ppPtr_(NULL),
cpPtr_(NULL),
cellTreePtr_(NULL),
labels_(0),
cellCentresPtr_(NULL),
@ -351,36 +345,4 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const
}
const Foam::indexedOctree<Foam::treeDataCell>&
Foam::primitiveMesh::cellTree() const
{
if (!cellTreePtr_)
{
treeBoundBox overallBb(points());
Random rndGen(261782);
overallBb = overallBb.extend(rndGen, 1E-4);
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
cellTreePtr_ =
new indexedOctree<treeDataCell>
(
treeDataCell
(
false, // not cache bb
*this
),
overallBb,
8, // maxLevel
10, // leafsize
5.0 // duplicity
);
}
return *cellTreePtr_;
}
// ************************************************************************* //

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

@ -69,10 +69,6 @@ SourceFiles
namespace Foam
{
// Forward declaration of classes
class treeDataCell;
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
Class primitiveMesh Declaration
\*---------------------------------------------------------------------------*/
@ -155,9 +151,6 @@ class primitiveMesh
//- Cell-points
mutable labelListList* cpPtr_;
//- Search tree to allow spatial tet searching
mutable indexedOctree<treeDataCell>* cellTreePtr_;
// On-the-fly edge addresing storage
@ -485,10 +478,6 @@ public:
const labelListList& cellPoints() const;
//- Build (if necessary) and return the cell search tree
const indexedOctree<treeDataCell>& cellTree() const;
// Geometric data (raw!)
const vectorField& cellCentres() const;
@ -821,9 +810,6 @@ public:
//- Clear topological data
void clearAddressing();
//- Clear cell tree data
void clearCellTree();
//- Clear all geometry and addressing unnecessary for CFD
void clearOut();
};

17
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C
View File

@ -24,11 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "primitiveMesh.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "demandDrivenData.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -102,11 +98,6 @@ void Foam::primitiveMesh::printAllocated() const
Pout<< " Cell-point" << endl;
}
if (cellTreePtr_)
{
Pout<< " Cell-tree" << endl;
}
// Geometry
if (cellCentresPtr_)
{
@ -173,14 +164,6 @@ void Foam::primitiveMesh::clearAddressing()
deleteDemandDrivenData(pePtr_);
deleteDemandDrivenData(ppPtr_);
deleteDemandDrivenData(cpPtr_);
deleteDemandDrivenData(cellTreePtr_);
}
void Foam::primitiveMesh::clearCellTree()
{
deleteDemandDrivenData(cellTreePtr_);
}

26
src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
View File

@ -680,13 +680,17 @@ void Foam::polyTopoChange::getFaceOrder
// First unassigned face
label newFaceI = 0;
labelList nbr;
labelList order;
forAll(cellMap_, cellI)
{
label startOfCell = cellFaceOffsets[cellI];
label nFaces = cellFaceOffsets[cellI+1] - startOfCell;
// Neighbouring cells
SortableList<label> nbr(nFaces);
//SortableList<label> nbr(nFaces);
nbr.setSize(nFaces);
for (label i = 0; i < nFaces; i++)
{
@ -725,14 +729,24 @@ void Foam::polyTopoChange::getFaceOrder
}
}
nbr.sort();
//nbr.sort();
order.setSize(nFaces);
sortedOrder(nbr, order);
forAll(nbr, i)
//forAll(nbr, i)
//{
// if (nbr[i] != -1)
// {
// oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
// newFaceI++;
// }
//}
forAll(order, i)
{
if (nbr[i] != -1)
label index = order[i];
if (nbr[index] != -1)
{
oldToNew[cellFaces[startOfCell + nbr.indices()[i]]] =
newFaceI++;
oldToNew[cellFaces[startOfCell + index]] = newFaceI++;
}
}
}

2
src/engine/ignition/ignitionSite.C
View File

@ -36,7 +36,7 @@ void Foam::ignitionSite::findIgnitionCells(const fvMesh& mesh)
const volVectorField& centres = mesh.C();
const scalarField& vols = mesh.V();
label ignCell = mesh.findCell(location_);
label ignCell = mesh.findCell(location_, polyMesh::FACEDIAGTETS);
if (ignCell == -1)
{
return;

6
src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C
View File

@ -107,7 +107,11 @@ void Foam::basicSource::setCellSet()
forAll(points_, i)
{
label cellI = mesh_.findCell(points_[i]);
label cellI = mesh_.findCell
(
points_[i],
polyMesh::FACEDIAGTETS
);
if (cellI >= 0)
{
selectedCells.insert(cellI);

6
src/finiteVolume/cfdTools/general/fieldSources/timeActivatedExplicitSource/TimeActivatedExplicitSource.C
View File

@ -234,7 +234,11 @@ void Foam::TimeActivatedExplicitSource<Type>::setCellSet()
forAll(points_, i)
{
label cellI = mesh_.findCell(points_[i]);
label cellI = mesh_.findCell
(
points_[i],
polyMesh::FACEDIAGTETS
);
if (cellI >= 0)
{
selectedCells.insert(cellI);

6
src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
View File

@ -76,7 +76,11 @@ void Foam::setRefCell
else if (dict.found(refPointName))
{
point refPointi(dict.lookup(refPointName));
refCelli = field.mesh().findCell(refPointi);
refCelli = field.mesh().findCell
(
refPointi,
polyMesh::FACEDIAGTETS
);
label hasRef = (refCelli >= 0 ? 1 : 0);
label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());
if (sumHasRef != 1)

10
src/lagrangian/basic/InteractionLists/InteractionLists.C
View File

@ -163,7 +163,13 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
indexedOctree<treeDataCell> coupledPatchRangeTree
(
treeDataCell(true, mesh_, coupledPatchRangeCells),
treeDataCell
(
true, // cache cell bb
mesh_,
coupledPatchRangeCells, // subset of mesh
polyMesh::FACEDIAGTETS // consistent with tracking
),
procBbRndExt,
8, // maxLevel,
10, // leafSize,
@ -382,7 +388,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
indexedOctree<treeDataCell> allCellsTree
(
treeDataCell(true, mesh_),
treeDataCell(true, mesh_, polyMesh::FACEDIAGTETS),
procBbRndExt,
8, // maxLevel,
10, // leafSize,

9
src/lagrangian/dieselSpray/spray/findInjectorCell.H
View File

@ -9,7 +9,8 @@ if (injectorCell >= 0)
foundCell = mesh_.pointInCell
(
injectionPosition,
injectorCell
injectorCell,
polyMesh::FACEDIAGTETS
);
}
@ -39,7 +40,8 @@ if (!foundCell)
foundCell = mesh_.pointInCell
(
injectionPosition,
injectorCell
injectorCell,
polyMesh::FACEDIAGTETS
);
}
reduce(foundCell, orOp<bool>());
@ -61,7 +63,8 @@ if (!foundCell)
foundCell = mesh_.pointInCell
(
injectionPosition,
injectorCell
injectorCell,
polyMesh::FACEDIAGTETS
);
}
reduce(foundCell, orOp<bool>());

10
src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
View File

@ -205,7 +205,15 @@ bool Foam::InjectionModel<CloudType>::findCellAtPosition
{
position += SMALL*(cellCentres[cellI] - position);
if (this->owner().mesh().pointInCell(position, cellI))
if
(
this->owner().mesh().pointInCell
(
position,
cellI,
polyMesh::FACEDIAGTETS
)
)
{
procI = Pstream::myProcNo();
}

2
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
View File

@ -86,7 +86,7 @@ Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
{
const point& keepPoint = keepPoints_[i];
label localCellI = mesh.findCell(keepPoint);
label localCellI = mesh.findCell(keepPoint, polyMesh::FACEDIAGTETS);
label globalCellI = -1;

2
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
View File

@ -1825,7 +1825,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
label regionI = -1;
label cellI = mesh_.findCell(keepPoint);
label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
if (cellI != -1)
{

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

@ -1248,7 +1248,7 @@ void Foam::meshRefinement::findCellZoneInsideWalk
// Find the region containing the insidePoint
label keepRegionI = -1;
label cellI = mesh_.findCell(insidePoint);
label cellI = mesh_.findCell(insidePoint, polyMesh::FACEDIAGTETS);
if (cellI != -1)
{
@ -1418,7 +1418,7 @@ void Foam::meshRefinement::findCellZoneTopo
// Find the region containing the keepPoint
label keepRegionI = -1;
label cellI = mesh_.findCell(keepPoint);
label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
if (cellI != -1)
{
@ -1959,7 +1959,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
// Find the region containing the keepPoint
label keepRegionI = -1;
label cellI = mesh_.findCell(keepPoint);
label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
if (cellI != -1)
{

4
src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
View File

@ -174,7 +174,7 @@ void Foam::mappedPatchBase::findSamples
}
// Octree based search engine
meshSearch meshSearchEngine(mesh, false);
meshSearch meshSearchEngine(mesh, polyMesh::FACEDIAGTETS);
forAll(samples, sampleI)
{
@ -291,7 +291,7 @@ void Foam::mappedPatchBase::findSamples
}
// Octree based search engine
meshSearch meshSearchEngine(mesh, false);
meshSearch meshSearchEngine(mesh, polyMesh::FACEDIAGTETS);
forAll(samples, sampleI)
{

191
src/meshTools/meshSearch/meshSearch.C
View File

@ -302,7 +302,7 @@ Foam::label Foam::meshSearch::findCellLinear(const point& location) const
while ((!cellFound) && (n < mesh_.nCells()))
{
if (pointInCell(location, n))
if (mesh_.pointInCell(location, n, cellDecompMode_))
{
cellFound = true;
cellI = n;
@ -338,7 +338,7 @@ Foam::label Foam::meshSearch::findCellWalk
) << "illegal seedCell:" << seedCellI << exit(FatalError);
}
if (pointInCell(location, seedCellI))
if (mesh_.pointInCell(location, seedCellI, cellDecompMode_))
{
return seedCellI;
}
@ -368,7 +368,7 @@ Foam::label Foam::meshSearch::findCellWalk
}
// Check if this is the correct cell
if (pointInCell(location, cellI))
if (mesh_.pointInCell(location, cellI, cellDecompMode_))
{
return cellI;
}
@ -500,10 +500,14 @@ Foam::vector Foam::meshSearch::offset
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::meshSearch::meshSearch(const polyMesh& mesh, const bool faceDecomp)
Foam::meshSearch::meshSearch
(
const polyMesh& mesh,
const polyMesh::cellRepresentation cellDecompMode
)
:
mesh_(mesh),
faceDecomp_(faceDecomp)
cellDecompMode_(cellDecompMode)
{}
@ -512,11 +516,11 @@ Foam::meshSearch::meshSearch
(
const polyMesh& mesh,
const treeBoundBox& bb,
const bool faceDecomp
const polyMesh::cellRepresentation cellDecompMode
)
:
mesh_(mesh),
faceDecomp_(faceDecomp)
cellDecompMode_(cellDecompMode)
{
overallBbPtr_.reset(new treeBoundBox(bb));
}
@ -615,8 +619,9 @@ const
(
treeDataCell
(
false, // not cache bb
mesh_
false, // not cache bb
mesh_,
cellDecompMode_ // cell decomposition mode for inside tests
),
overallBbPtr_(),
8, // maxLevel
@ -630,90 +635,90 @@ const
}
// Is the point in the cell
// Works by checking if there is a face inbetween the point and the cell
// centre.
// Check for internal uses proper face decomposition or just average normal.
bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
{
if (faceDecomp_)
{
const point& ctr = mesh_.cellCentres()[cellI];
vector dir(p - ctr);
scalar magDir = mag(dir);
// Check if any faces are hit by ray from cell centre to p.
// If none -> p is in cell.
const labelList& cFaces = mesh_.cells()[cellI];
// Make sure half_ray does not pick up any faces on the wrong
// side of the ray.
scalar oldTol = intersection::setPlanarTol(0.0);
forAll(cFaces, i)
{
label faceI = cFaces[i];
pointHit inter = mesh_.faces()[faceI].ray
(
ctr,
dir,
mesh_.points(),
intersection::HALF_RAY,
intersection::VECTOR
);
if (inter.hit())
{
scalar dist = inter.distance();
if (dist < magDir)
{
// Valid hit. Hit face so point is not in cell.
intersection::setPlanarTol(oldTol);
return false;
}
}
}
intersection::setPlanarTol(oldTol);
// No face inbetween point and cell centre so point is inside.
return true;
}
else
{
const labelList& f = mesh_.cells()[cellI];
const labelList& owner = mesh_.faceOwner();
const vectorField& cf = mesh_.faceCentres();
const vectorField& Sf = mesh_.faceAreas();
forAll(f, facei)
{
label nFace = f[facei];
vector proj = p - cf[nFace];
vector normal = Sf[nFace];
if (owner[nFace] == cellI)
{
if ((normal & proj) > 0)
{
return false;
}
}
else
{
if ((normal & proj) < 0)
{
return false;
}
}
}
return true;
}
}
//// Is the point in the cell
//// Works by checking if there is a face inbetween the point and the cell
//// centre.
//// Check for internal uses proper face decomposition or just average normal.
//bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
//{
// if (faceDecomp_)
// {
// const point& ctr = mesh_.cellCentres()[cellI];
//
// vector dir(p - ctr);
// scalar magDir = mag(dir);
//
// // Check if any faces are hit by ray from cell centre to p.
// // If none -> p is in cell.
// const labelList& cFaces = mesh_.cells()[cellI];
//
// // Make sure half_ray does not pick up any faces on the wrong
// // side of the ray.
// scalar oldTol = intersection::setPlanarTol(0.0);
//
// forAll(cFaces, i)
// {
// label faceI = cFaces[i];
//
// pointHit inter = mesh_.faces()[faceI].ray
// (
// ctr,
// dir,
// mesh_.points(),
// intersection::HALF_RAY,
// intersection::VECTOR
// );
//
// if (inter.hit())
// {
// scalar dist = inter.distance();
//
// if (dist < magDir)
// {
// // Valid hit. Hit face so point is not in cell.
// intersection::setPlanarTol(oldTol);
//
// return false;
// }
// }
// }
//
// intersection::setPlanarTol(oldTol);
//
// // No face inbetween point and cell centre so point is inside.
// return true;
// }
// else
// {
// const labelList& f = mesh_.cells()[cellI];
// const labelList& owner = mesh_.faceOwner();
// const vectorField& cf = mesh_.faceCentres();
// const vectorField& Sf = mesh_.faceAreas();
//
// forAll(f, facei)
// {
// label nFace = f[facei];
// vector proj = p - cf[nFace];
// vector normal = Sf[nFace];
// if (owner[nFace] == cellI)
// {
// if ((normal & proj) > 0)
// {
// return false;
// }
// }
// else
// {
// if ((normal & proj) < 0)
// {
// return false;
// }
// }
// }
//
// return true;
// }
//}
Foam::label Foam::meshSearch::findNearestCell

25
src/meshTools/meshSearch/meshSearch.H
View File

@ -38,6 +38,7 @@ SourceFiles
#include "pointIndexHit.H"
#include "pointField.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -45,7 +46,6 @@ namespace Foam
{
// Forward declaration of classes
class polyMesh;
class treeDataCell;
class treeDataFace;
template<class Type> class indexedOctree;
@ -62,8 +62,8 @@ class meshSearch
//- Reference to mesh
const polyMesh& mesh_;
//- Whether to use face decomposition for all geometric tests
const bool faceDecomp_;
//- Whether to use cell decomposition for all geometric tests
const polyMesh::cellRepresentation cellDecompMode_;
//- data bounding box
mutable autoPtr<treeBoundBox> overallBbPtr_;
@ -168,7 +168,11 @@ public:
//- Construct from components. Constructs bb slightly bigger than
// mesh points bb.
meshSearch(const polyMesh& mesh, const bool faceDecomp = true);
meshSearch
(
const polyMesh& mesh,
const polyMesh::cellRepresentation = polyMesh::FACEDIAGTETS
);
//- Construct with a custom bounding box. Any mesh element outside
// bb will not be found. Up to user to make sure bb
@ -177,7 +181,7 @@ public:
(
const polyMesh& mesh,
const treeBoundBox& bb,
const bool faceDecomp = true
const polyMesh::cellRepresentation = polyMesh::FACEDIAGTETS
);
//- Destructor
@ -193,6 +197,11 @@ public:
return mesh_;
}
polyMesh::cellRepresentation decompMode() const
{
return cellDecompMode_;
}
//- Get (demand driven) reference to octree holding all
// boundary faces
const indexedOctree<treeDataFace>& boundaryTree() const;
@ -203,10 +212,6 @@ public:
// Queries
//- test for point in cell. Does not handle cells with center
// outside cell.
bool pointInCell(const point& p, const label celli) const;
//- Find nearest cell in terms of cell centre.
// Options:
// - use octree
@ -227,7 +232,7 @@ public:
const bool useTreeSearch = true
) const;
//- Find cell containing (using pointInCell) location.
//- Find cell containing location.
// If seed provided walks and falls back to linear/tree search.
// (so handles holes correctly)s
// Returns -1 if not in domain.

7
src/meshTools/octree/octreeDataCell.C
View File

@ -116,7 +116,12 @@ bool Foam::octreeDataCell::contains
const point& sample
) const
{
return mesh_.pointInCell(sample, cellLabels_[index]);
return mesh_.pointInCell
(
sample,
cellLabels_[index],
polyMesh::FACEDIAGTETS
);
}

2
src/meshTools/sets/cellSources/nearestToCell/nearestToCell.C
View File

@ -25,8 +25,6 @@ License
#include "nearestToCell.H"
#include "polyMesh.H"
#include "meshSearch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

2
src/meshTools/sets/cellSources/regionToCell/regionToCell.C
View File

@ -60,7 +60,7 @@ Foam::topoSetSource::addToUsageTable Foam::regionToCell::usage_
void Foam::regionToCell::combine(topoSet& set, const bool add) const
{
label cellI = mesh_.findCell(insidePoint_);
label cellI = mesh_.findCell(insidePoint_, polyMesh::FACEDIAGTETS);
// Load the subset of cells
boolList blockedFace(mesh_.nFaces(), false);

2
src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
View File

@ -166,7 +166,7 @@ void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
// Construct search engine on mesh
meshSearch queryMesh(mesh_, true);
meshSearch queryMesh(mesh_, polyMesh::FACEDIAGTETS);
// Check all 'outside' points

2
src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C
View File

@ -25,8 +25,6 @@ License
#include "nearestToPoint.H"
#include "polyMesh.H"
#include "meshSearch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

1
src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
View File

@ -25,7 +25,6 @@ License
#include "surfaceToPoint.H"
#include "polyMesh.H"
#include "meshSearch.H"
#include "triSurfaceSearch.H"
#include "cpuTime.H"

2
src/meshTools/surfaceSets/surfaceSets.C
View File

@ -235,7 +235,7 @@ void Foam::surfaceSets::getSurfaceSets
)
{
// Construct search engine on mesh
meshSearch queryMesh(mesh, true);
meshSearch queryMesh(mesh, polyMesh::FACEDIAGTETS);
// Cut faces with surface and classify cells
cellClassification cellType

2
src/postProcessing/functionObjects/field/streamLine/streamLine.C
View File

@ -363,7 +363,7 @@ void Foam::streamLine::read(const dictionary& dict)
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
meshSearchPtr_.reset(new meshSearch(mesh, false));
meshSearchPtr_.reset(new meshSearch(mesh, polyMesh::FACEDIAGTETS));
const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
sampledSetPtr_ = sampledSet::New

3
src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C
View File

@ -129,7 +129,8 @@ Foam::codedFunctionObject::codedFunctionObject
:
functionObject(name),
codedBase(),
time_(time)
time_(time),
dict_(dict)
{
if (readNow)
{

24
src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
View File

@ -102,7 +102,7 @@ void Foam::meshToMesh::calcAddressing()
indexedOctree<treeDataCell> oc
(
treeDataCell(false, fromMesh_),
treeDataCell(false, fromMesh_, polyMesh::FACEDIAGTETS),
shiftedBb, // overall bounding box
8, // maxLevel
10, // leafsize
@ -267,7 +267,7 @@ void Foam::meshToMesh::cellAddresses
cellAddressing_[toI] = -1;
// Check point is actually in the nearest cell
if (fromMesh.pointInCell(p, curCell))
if (fromMesh.pointInCell(p, curCell, polyMesh::FACEDIAGTETS))
{
cellAddressing_[toI] = curCell;
}
@ -292,7 +292,15 @@ void Foam::meshToMesh::cellAddresses
{
// search through all the neighbours.
// If point is in neighbour reset current cell
if (fromMesh.pointInCell(p, neighbours[nI]))
if
(
fromMesh.pointInCell
(
p,
neighbours[nI],
polyMesh::FACEDIAGTETS
)
)
{
cellAddressing_[toI] = neighbours[nI];
found = true;
@ -316,7 +324,15 @@ void Foam::meshToMesh::cellAddresses
{
// search through all the neighbours.
// If point is in neighbour reset current cell
if (fromMesh.pointInCell(p, nn[nI]))
if
(
fromMesh.pointInCell
(
p,
nn[nI],
polyMesh::FACEDIAGTETS
)
)
{
cellAddressing_[toI] = nn[nI];
found = true;

1
src/sampling/probes/patchProbes.C
View File

@ -28,7 +28,6 @@ License
#include "IOmanip.H"
// For 'nearInfo' helper class only
#include "mappedPatchBase.H"
//#include "meshSearch.H"
#include "treeBoundBox.H"
#include "treeDataFace.H"

2
src/sampling/probes/probes.C
View File

@ -45,7 +45,7 @@ void Foam::probes::findElements(const fvMesh& mesh)
{
const vector& location = operator[](probeI);
elementList_[probeI] = mesh.findCell(location);
elementList_[probeI] = mesh.findCell(location, polyMesh::FACEDIAGTETS);
if (debug && elementList_[probeI] != -1)
{

17
src/sampling/sampledSet/sampledSet/sampledSet.C
View File

@ -68,7 +68,7 @@ Foam::label Foam::sampledSet::getCell
{
label cellI = getBoundaryCell(faceI);
if (!mesh().pointInCell(sample, cellI))
if (!mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
{
FatalErrorIn
(
@ -85,7 +85,7 @@ Foam::label Foam::sampledSet::getCell
label cellI = mesh().faceOwner()[faceI];
if (mesh().pointInCell(sample, cellI))
if (mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
{
return cellI;
}
@ -93,7 +93,7 @@ Foam::label Foam::sampledSet::getCell
{
cellI = mesh().faceNeighbour()[faceI];
if (mesh().pointInCell(sample, cellI))
if (mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
{
return cellI;
}
@ -261,12 +261,17 @@ bool Foam::sampledSet::getTrackingPoint
if (bFaceI == -1)
{
// No boundary intersection. Try and find cell samplePt is in
trackCellI = mesh().findCell(samplePt);
trackCellI = mesh().findCell(samplePt, searchEngine_.decompMode());
if
(
(trackCellI == -1)
|| !mesh().pointInCell(samplePt, trackCellI)
|| !mesh().pointInCell
(
samplePt,
trackCellI,
searchEngine_.decompMode()
)
)
{
// Line samplePt - end_ does not intersect domain at all.
@ -316,7 +321,7 @@ bool Foam::sampledSet::getTrackingPoint
// samplePt inside or marginally outside.
trackPt = samplePt;
trackFaceI = -1;
trackCellI = mesh().findCell(trackPt);
trackCellI = mesh().findCell(trackPt, searchEngine_.decompMode());
isGoodSample = true;
}

2
src/sampling/sampledSet/sampledSets/sampledSets.C
View File

@ -138,7 +138,7 @@ Foam::sampledSets::sampledSets
mesh_(refCast<const fvMesh>(obr)),
loadFromFiles_(loadFromFiles),
outputPath_(fileName::null),
searchEngine_(mesh_, true),
searchEngine_(mesh_, polyMesh::FACEDIAGTETS),
interpolationScheme_(word::null),
writeFormat_(word::null)
{

6
src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
View File

@ -53,7 +53,11 @@ void Foam::triSurfaceMeshPointSet::calcSamples
{
forAll(sampleCoords_, sampleI)
{
label cellI = searchEngine().findCell(sampleCoords_[sampleI]);
label cellI = searchEngine().findCell
(
sampleCoords_[sampleI],
polyMesh::FACEDIAGTETS
);
if (cellI != -1)
{

12
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
View File

@ -241,7 +241,7 @@ bool Foam::sampledTriSurfaceMesh::update()
const pointField& fc = surface_.faceCentres();
// Mesh search engine, no triangulation of faces.
meshSearch meshSearcher(mesh(), false);
meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
List<nearInfo> nearest(fc.size());
@ -435,7 +435,15 @@ bool Foam::sampledTriSurfaceMesh::update()
sampleElements_[pointI] = cellI;
// Check if point inside cell
if (mesh().pointInCell(pt, sampleElements_[pointI]))
if
(
mesh().pointInCell
(
pt,
sampleElements_[pointI],
meshSearcher.decompMode()
)
)
{
samplePoints_[pointI] = pt;
}

3
src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H
View File

@ -135,6 +135,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

10
src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H
View File

@ -89,6 +89,16 @@ Foam::eConstThermo<EquationOfState>::New(const dictionary& dict)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState>
inline Foam::scalar Foam::eConstThermo<EquationOfState>::limit
(
const scalar T
) const
{
return T;
}
template<class EquationOfState>
inline Foam::scalar Foam::eConstThermo<EquationOfState>::cp
(

3
src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H
View File

@ -133,6 +133,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

10
src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H
View File

@ -89,6 +89,16 @@ Foam::hConstThermo<equationOfState>::New(const dictionary& dict)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState>
inline Foam::scalar Foam::hConstThermo<EquationOfState>::limit
(
const scalar T
) const
{
return T;
}
template<class equationOfState>
inline Foam::scalar Foam::hConstThermo<equationOfState>::cp
(

3
src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H
View File

@ -151,6 +151,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

10
src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H
View File

@ -82,6 +82,16 @@ inline Foam::hPolynomialThermo<EquationOfState, PolySize>::hPolynomialThermo
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState, int PolySize>
inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::limit
(
const scalar T
) const
{
return T;
}
template<class EquationOfState, int PolySize>
inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::cp
(

6
src/thermophysicalModels/specie/thermo/janaf/janafThermo.H
View File

@ -119,9 +119,6 @@ private:
//- Check that input data is valid
void checkInputData() const;
//- Check given temperature is within the range of the fitted coeffs
inline void checkT(const scalar T) const;
//- Return the coefficients corresponding to the given temperature
inline const coeffArray& coeffs(const scalar T) const;
@ -153,6 +150,9 @@ public:
// Member Functions
//- Limit the temperature to be in the range Tlow_ to Thigh_
inline scalar limit(const scalar T) const;
// Fundamental properties
//- Heat capacity at constant pressure [J/(kmol K)]

103
src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H
View File

@ -52,22 +52,6 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
}
template<class EquationOfState>
inline void Foam::janafThermo<EquationOfState>::checkT(const scalar T) const
{
if (T < Tlow_ || T > Thigh_)
{
FatalErrorIn
(
"janafThermo<EquationOfState>::checkT(const scalar T) const"
) << "attempt to use janafThermo<EquationOfState>"
" out of temperature range "
<< Tlow_ << " -> " << Thigh_ << "; T = " << T
<< abort(FatalError);
}
}
template<class EquationOfState>
inline const typename Foam::janafThermo<EquationOfState>::coeffArray&
Foam::janafThermo<EquationOfState>::coeffs
@ -75,8 +59,6 @@ Foam::janafThermo<EquationOfState>::coeffs
const scalar T
) const
{
checkT(T);
if (T < Tcommon_)
{
return lowCpCoeffs_;
@ -112,6 +94,31 @@ inline Foam::janafThermo<EquationOfState>::janafThermo
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::limit
(
const scalar T
) const
{
if (T < Tlow_ || T > Thigh_)
{
WarningIn
(
"janafThermo<EquationOfState>::limit(const scalar T) const"
) << "attempt to use janafThermo<EquationOfState>"
" out of temperature range "
<< Tlow_ << " -> " << Thigh_ << "; T = " << T
<< endl;
return min(max(T, Tlow_), Thigh_);
}
else
{
return T;
}
}
template<class EquationOfState>
inline Foam::scalar Foam::janafThermo<EquationOfState>::cp
(
@ -196,7 +203,19 @@ inline void Foam::janafThermo<EquationOfState>::operator+=
Tlow_ = max(Tlow_, jt.Tlow_);
Thigh_ = min(Thigh_, jt.Thigh_);
Tcommon_ = molr1*Tcommon_ + molr2*jt.Tcommon_;
if (notEqual(Tcommon_, jt.Tcommon_))
{
FatalErrorIn
(
"janafThermo<EquationOfState>::operator+="
"(const janafThermo<EquationOfState>& jt) const"
) << "Tcommon " << Tcommon_ << " for "
<< (this->name().size() ? this->name() : "others")
<< " != " << jt.Tcommon_ << " for "
<< (jt.name().size() ? jt.name() : "others")
<< exit(FatalError);
}
for
(
@ -231,7 +250,19 @@ inline void Foam::janafThermo<EquationOfState>::operator-=
Tlow_ = max(Tlow_, jt.Tlow_);
Thigh_ = min(Thigh_, jt.Thigh_);
Tcommon_ = molr1*Tcommon_ - molr2*jt.Tcommon_;
if (notEqual(Tcommon_, jt.Tcommon_))
{
FatalErrorIn
(
"janafThermo<EquationOfState>::operator-="
"(const janafThermo<EquationOfState>& jt) const"
) << "Tcommon " << Tcommon_ << " for "
<< (this->name().size() ? this->name() : "others")
<< " != " << jt.Tcommon_ << " for "
<< (jt.name().size() ? jt.name() : "others")
<< exit(FatalError);
}
for
(
@ -285,12 +316,26 @@ inline Foam::janafThermo<EquationOfState> Foam::operator+
+ molr2*jt2.lowCpCoeffs_[coefLabel];
}
if (notEqual(jt1.Tcommon_, jt2.Tcommon_))
{
FatalErrorIn
(
"operator+"
"(const janafThermo<EquationOfState>& jt1,"
" const janafThermo<EquationOfState>& jt2)"
) << "Tcommon " << jt1.Tcommon_ << " for "
<< (jt1.name().size() ? jt1.name() : "others")
<< " != " << jt2.Tcommon_ << " for "
<< (jt2.name().size() ? jt2.name() : "others")
<< exit(FatalError);
}
return janafThermo<EquationOfState>
(
eofs,
max(jt1.Tlow_, jt2.Tlow_),
min(jt1.Thigh_, jt2.Thigh_),
molr1*jt1.Tcommon_ + molr2*jt2.Tcommon_,
jt1.Tcommon_,
highCpCoeffs,
lowCpCoeffs
);
@ -329,12 +374,26 @@ inline Foam::janafThermo<EquationOfState> Foam::operator-
- molr2*jt2.lowCpCoeffs_[coefLabel];
}
if (notEqual(jt1.Tcommon_, jt2.Tcommon_))
{
FatalErrorIn
(
"operator-"
"(const janafThermo<EquationOfState>& jt1,"
" const janafThermo<EquationOfState>& jt2)"
) << "Tcommon " << jt1.Tcommon_ << " for "
<< (jt1.name().size() ? jt1.name() : "others")
<< " != " << jt2.Tcommon_ << " for "
<< (jt2.name().size() ? jt2.name() : "others")
<< exit(FatalError);
}
return janafThermo<EquationOfState>
(
eofs,
max(jt1.Tlow_, jt2.Tlow_),
min(jt1.Thigh_, jt2.Thigh_),
molr1*jt1.Tcommon_ - molr2*jt2.Tcommon_,
jt1.Tcommon_,
highCpCoeffs,
lowCpCoeffs
);

5
src/thermophysicalModels/specie/thermo/specieThermo/specieThermo.H
View File

@ -103,14 +103,15 @@ class specieThermo
// Private Member Functions
//- return the temperature corresponding to the value of the
//- Return the temperature corresponding to the value of the
// thermodynamic property f, given the function f = F(T) and dF(T)/dT
inline scalar T
(
scalar f,
scalar T0,
scalar (specieThermo::*F)(const scalar) const,
scalar (specieThermo::*dFdT)(const scalar) const
scalar (specieThermo::*dFdT)(const scalar) const,
scalar (specieThermo::*limit)(const scalar) const
) const;

33
src/thermophysicalModels/specie/thermo/specieThermo/specieThermoI.H
View File

@ -43,7 +43,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
scalar f,
scalar T0,
scalar (specieThermo<Thermo>::*F)(const scalar) const,
scalar (specieThermo<Thermo>::*dFdT)(const scalar) const
scalar (specieThermo<Thermo>::*dFdT)(const scalar) const,
scalar (specieThermo<Thermo>::*limit)(const scalar) const
) const
{
scalar Test = T0;
@ -54,7 +55,8 @@ inline Foam::scalar Foam::specieThermo<Thermo>::T
do
{
Test = Tnew;
Tnew = Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test);
Tnew =
(this->*limit)(Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test));
if (iter++ > maxIter_)
{
@ -276,7 +278,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TH
const scalar T0
) const
{
return T(h, T0, &specieThermo<Thermo>::H, &specieThermo<Thermo>::Cp);
return T
(
h,
T0,
&specieThermo<Thermo>::H,
&specieThermo<Thermo>::Cp,
&specieThermo<Thermo>::limit
);
}
@ -287,7 +296,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::THs
const scalar T0
) const
{
return T(hs, T0, &specieThermo<Thermo>::Hs, &specieThermo<Thermo>::Cp);
return T
(
hs,
T0,
&specieThermo<Thermo>::Hs,
&specieThermo<Thermo>::Cp,
&specieThermo<Thermo>::limit
);
}
@ -298,7 +314,14 @@ inline Foam::scalar Foam::specieThermo<Thermo>::TE
const scalar T0
) const
{
return T(e, T0, &specieThermo<Thermo>::E, &specieThermo<Thermo>::Cv);
return T
(
e,
T0,
&specieThermo<Thermo>::E,
&specieThermo<Thermo>::Cv,
&specieThermo<Thermo>::limit
);
}

15
src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
View File

@ -70,17 +70,10 @@ Foam::porousBafflePressureFvPatchField<Type>::porousBafflePressureFvPatchField
I_(readScalar(dict.lookup("I"))),
length_(readScalar(dict.lookup("length")))
{
if (dict.found("value"))
{
fvPatchField<Type>::operator=
(
Field<Type>("value", dict, p.size())
);
}
else
{
this->evaluate(Pstream::blocking);
}
fvPatchField<Type>::operator=
(
Field<Type>("value", dict, p.size())
);
}

1
src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
View File

@ -426,6 +426,7 @@ void LienCubicKELowRe::correct()
volScalarField G
(
"RASModel::G",
Cmu_*fMu*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU_)
);

1
src/turbulenceModels/incompressible/RAS/Make/files
View File

@ -16,6 +16,7 @@ LienCubicKELowRe/LienCubicKELowRe.C
NonlinearKEShih/NonlinearKEShih.C
LienLeschzinerLowRe/LienLeschzinerLowRe.C
LamBremhorstKE/LamBremhorstKE.C
kkLOmega/kkLOmega.C
/* Wall functions */
wallFunctions = derivedFvPatchFields/wallFunctions

779
src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.C Normal file
View File

@ -0,0 +1,779 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "kkLOmega.H"
#include "addToRunTimeSelectionTable.H"
#include "backwardsCompatibilityWallFunctions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(kkLOmega, 0);
addToRunTimeSelectionTable(RASModel, kkLOmega, dictionary);
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<volScalarField> kkLOmega::fv(const volScalarField& Ret) const
{
return(1.0 - exp(-sqrt(Ret)/Av_));
}
tmp<volScalarField> kkLOmega::fINT() const
{
return
(
min
(
kl_/(Cint_*(kl_ + kt_ + kMin_)),
dimensionedScalar("1.0", dimless, 1.0)
)
);
}
tmp<volScalarField> kkLOmega::fSS(const volScalarField& omega) const
{
return(exp(-sqr(Css_*nu()*omega/(kt_ + kMin_))));
}
tmp<volScalarField> kkLOmega::Cmu(const volScalarField& S) const
{
return(1.0/(A0_ + As_*(S/(omega_ + omegaMin_))));
}
tmp<volScalarField> kkLOmega::BetaTS(const volScalarField& Rew) const
{
return(scalar(1.0) - exp(-sqr(max(Rew - CtsCrit_, 0.0))/Ats_));
}
tmp<volScalarField> kkLOmega::fTaul
(
const volScalarField& lambdaEff,
const volScalarField& ktL
) const
{
return
(
scalar(1.0)
- exp
(
-CtauL_*ktL
/
(
sqr
(
lambdaEff*omega_
+ dimensionedScalar
(
"ROTVSMALL",
dimLength*inv(dimTime),
ROOTVSMALL
)
)
)
)
);
}
tmp<volScalarField> kkLOmega::alphaT
(
const volScalarField& lambdaEff,
const volScalarField& fv,
const volScalarField& ktS
) const
{
return(fv*CmuStd_*sqrt(ktS)*lambdaEff);
}
tmp<volScalarField> kkLOmega::fOmega
(
const volScalarField& lambdaEff,
const volScalarField& lambdaT
) const
{
return
(
scalar(1.0)
- exp
(
-0.41
* pow4
(
lambdaEff
/ (
lambdaT
+ dimensionedScalar
(
"ROTVSMALL",
lambdaT.dimensions(),
ROOTVSMALL
)
)
)
)
);
}
tmp<volScalarField> kkLOmega::gammaBP(const volScalarField& omega) const
{
return
(
max
(
kt_/nu()
/
(
omega
+ dimensionedScalar("ROTVSMALL", omega.dimensions(), ROOTVSMALL)
)
-
CbpCrit_
,
0.0
)
);
}
tmp<volScalarField> kkLOmega::gammaNAT
(
const volScalarField& ReOmega,
const volScalarField& fNatCrit
) const
{
return
(
max
(
ReOmega
- CnatCrit_
/ (
fNatCrit + dimensionedScalar("ROTVSMALL", dimless, ROOTVSMALL)
)
,
0.0
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
kkLOmega::kkLOmega
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName,
const word& modelName
)
:
RASModel(modelName, U, phi, transport, turbulenceModelName),
A0_
(
dimensioned<scalar>::lookupOrAddToDict
(
"A0",
coeffDict_,
4.04
)
),
As_
(
dimensioned<scalar>::lookupOrAddToDict
(
"As",
coeffDict_,
2.12
)
),
Av_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Av",
coeffDict_,
6.75
)
),
Abp_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Abp",
coeffDict_,
0.6
)
),
Anat_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Anat",
coeffDict_,
200
)
),
Ats_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ats",
coeffDict_,
200
)
),
CbpCrit_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CbpCrit",
coeffDict_,
1.2
)
),
Cnc_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cnc",
coeffDict_,
0.1
)
),
CnatCrit_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CnatCrit",
coeffDict_,
1250
)
),
Cint_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cint",
coeffDict_,
0.75
)
),
CtsCrit_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CtsCrit",
coeffDict_,
1000
)
),
CrNat_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CrNat",
coeffDict_,
0.02
)
),
C11_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C11",
coeffDict_,
3.4e-6
)
),
C12_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C12",
coeffDict_,
1.0e-10
)
),
CR_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CR",
coeffDict_,
0.12
)
),
CalphaTheta_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CalphaTheta",
coeffDict_,
0.035
)
),
Css_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Css",
coeffDict_,
1.5
)
),
CtauL_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CtauL",
coeffDict_,
4360
)
),
Cw1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw1",
coeffDict_,
0.44
)
),
Cw2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw2",
coeffDict_,
0.92
)
),
Cw3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw3",
coeffDict_,
0.3
)
),
CwR_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CwR",
coeffDict_,
1.5
)
),
Clambda_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Clambda",
coeffDict_,
2.495
)
),
CmuStd_
(
dimensioned<scalar>::lookupOrAddToDict
(
"CmuStd",
coeffDict_,
0.09
)
),
Prtheta_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prtheta",
coeffDict_,
0.85
)
),
Sigmak_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Sigmak",
coeffDict_,
1
)
),
Sigmaw_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Sigmaw",
coeffDict_,
1.17
)
),
kt_
(
IOobject
(
"kt",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateK("kt", mesh_)
),
omega_
(
IOobject
(
"omega",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateOmega("omega", mesh_)
),
kl_
(
IOobject
(
"kl",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateK("kl", mesh_)
),
nut_
(
IOobject
(
"nut",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateNut("nut", mesh_)
),
y_(mesh_)
{
bound(kt_, kMin_);
bound(kl_, kMin_);
bound(omega_, omegaMin_);
nut_ = kt_/(omega_ + omegaMin_);
nut_.correctBoundaryConditions();
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<volSymmTensorField> kkLOmega::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*(kt_) - nut_*twoSymm(fvc::grad(U_)),
kt_.boundaryField().types()
)
);
}
tmp<volSymmTensorField> kkLOmega::devReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
-nuEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}
tmp<fvVectorMatrix> kkLOmega::divDevReff(volVectorField& U) const
{
return
(
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(T(fvc::grad(U))))
);
}
bool kkLOmega::read()
{
if (RASModel::read())
{
A0_.readIfPresent(coeffDict());
As_.readIfPresent(coeffDict());
Av_.readIfPresent(coeffDict());
Abp_.readIfPresent(coeffDict());
Anat_.readIfPresent(coeffDict());
Abp_.readIfPresent(coeffDict());
Ats_.readIfPresent(coeffDict());
CbpCrit_.readIfPresent(coeffDict());
Cnc_.readIfPresent(coeffDict());
CnatCrit_.readIfPresent(coeffDict());
Cint_.readIfPresent(coeffDict());
CtsCrit_.readIfPresent(coeffDict());
CrNat_.readIfPresent(coeffDict());
C11_.readIfPresent(coeffDict());
C12_.readIfPresent(coeffDict());
CR_.readIfPresent(coeffDict());
CalphaTheta_.readIfPresent(coeffDict());
Css_.readIfPresent(coeffDict());
CtauL_.readIfPresent(coeffDict());
Cw1_.readIfPresent(coeffDict());
Cw2_.readIfPresent(coeffDict());
Cw3_.readIfPresent(coeffDict());
CwR_.readIfPresent(coeffDict());
Clambda_.readIfPresent(coeffDict());
CmuStd_.readIfPresent(coeffDict());
Prtheta_.readIfPresent(coeffDict());
Sigmak_.readIfPresent(coeffDict());
Sigmaw_.readIfPresent(coeffDict());
return true;
}
else
{
return false;
}
}
void kkLOmega::correct()
{
RASModel::correct();
if (!turbulence_)
{
return;
}
if (mesh_.changing())
{
y_.correct();
y_.boundaryField() = max(y_.boundaryField(), VSMALL);
}
const volScalarField kT(kt_ + kl_);
const volScalarField lambdaT(sqrt(kT)/(omega_ + omegaMin_));
const volScalarField lambdaEff(min(Clambda_*y_, lambdaT));
const volScalarField fw
(
lambdaEff/(lambdaT + dimensionedScalar("SMALL", dimLength, ROOTVSMALL))
);
const volTensorField gradU(fvc::grad(U_));
const volScalarField omega(sqrt(2.0)*mag(skew(gradU)));
const volScalarField S2(2.0*magSqr(symm(gradU)));
const volScalarField ktS(fSS(omega)*fw*kt_);
const volScalarField nuts
(
fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_))
*fINT()
*Cmu(sqrt(S2))*sqrt(ktS)*lambdaEff
);
const volScalarField Pkt(nuts*S2);
const volScalarField ktL(kt_ - ktS);
const volScalarField ReOmega(sqr(y_)*omega/nu());
const volScalarField nutl
(
min
(
C11_*fTaul(lambdaEff, ktL)*omega*sqr(lambdaEff)
* sqrt(ktL)*lambdaEff/nu()
+ C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*omega
,
0.5*(kl_ + ktL)/sqrt(S2)
)
);
const volScalarField Pkl(nutl*S2);
const volScalarField alphaTEff
(
alphaT(lambdaEff, fv(sqr(fw)*kt_/nu()/(omega_ + omegaMin_)), ktS)
);
// By pass s0urce term divided by kl_
const dimensionedScalar fwMin("SMALL", dimless, ROOTVSMALL);
const volScalarField Rbp
(
CR_*(1.0 - exp(-gammaBP(omega)()/Abp_))*omega_
/ (fw + fwMin)
);
const volScalarField fNatCrit(1.0 - exp(-Cnc_*sqrt(kl_)*y_/nu()));
// Natural source term divided by kl_
const volScalarField Rnat
(
CrNat_*(1.0 - exp(-gammaNAT(ReOmega, fNatCrit)/Anat_))*omega
);
const volScalarField Dt(nu()*magSqr(fvc::grad(sqrt(kt_))));
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> ktEqn
(
fvm::ddt(kt_)
+ fvm::div(phi_, kt_)
- fvm::Sp(fvc::div(phi_), kt_)
- fvm::laplacian(DkEff(alphaTEff), kt_, "laplacian(alphaTEff,kt)")
==
Pkt
+ (Rbp + Rnat)*kl_
- Dt
- fvm::Sp(omega_, kt_)
);
ktEqn().relax();
ktEqn().boundaryManipulate(kt_.boundaryField());
solve(ktEqn);
bound(kt_, kMin_);
const volScalarField Dl(nu()*magSqr(fvc::grad(sqrt(kl_))));
// Laminar kinetic energy equation
tmp<fvScalarMatrix> klEqn
(
fvm::ddt(kl_)
+ fvm::div(phi_, kl_)
- fvm::Sp(fvc::div(phi_), kl_)
- fvm::laplacian(nu(), kl_, "laplacian(nu,kl)")
==
Pkl
- fvm::Sp(Rbp, kl_)
- fvm::Sp(Rnat, kl_)
- Dl
);
klEqn().relax();
klEqn().boundaryManipulate(kl_.boundaryField());
solve(klEqn);
bound(kl_, kMin_);
omega_.boundaryField().updateCoeffs();
// Turbulence specific dissipation rate equation
tmp<fvScalarMatrix> omegaEqn
(
fvm::ddt(omega_)
+ fvm::div(phi_, omega_)
- fvm::Sp(fvc::div(phi_), omega_)
- fvm::laplacian
(
DomegaEff(alphaTEff),
omega_,
"laplacian(alphaTEff,omega)"
)
==
Cw1_*Pkt*omega_/(kt_ + kMin_)
+ fvm::SuSp
(
(CwR_/(fw + fwMin) - 1.0)*kl_*(Rbp + Rnat)/(kt_ + kMin_)
, omega_
)
- fvm::Sp(Cw2_*omega_, omega_)
+ Cw3_*fOmega(lambdaEff, lambdaT)*alphaTEff*sqr(fw)*sqrt(kt_)/pow3(y_)
);
omegaEqn().relax();
omegaEqn().boundaryManipulate(omega_.boundaryField());
solve(omegaEqn);
bound(omega_, omegaMin_);
// Re-calculate viscosity
nut_ = nuts + nutl;
nut_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

296
src/turbulenceModels/incompressible/RAS/kkLOmega/kkLOmega.H Normal file
View File

@ -0,0 +1,296 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::incompressible::RASModels::kkLOmega
Description
Low Reynolds-number k-kl-omega turbulence model for
incompressible flows.
Turbulence model described in:
\verbatim
D. Keith Walters, Davor Cokljat
"A Three-Equation Eddy-Viscosity Model for Reynold-Averaged
Navier-Stokes Simulations of Transitional Flow"
\endverbatim
The default model coefficients correspond to the following:
\verbatim
kkLOmegaCoeffs
{
A0 4.04
As 2.12
Av 6.75
Abp 0.6
Anat 200
Ats 200
CbpCrit 1.2
Cnc 0.1
CnatCrit 1250
Cint 0.75
CtsCrit 1000
CrNat 0.02
C11 3.4e-6
C12 1.0e-10
CR 0.12
CalphaTheta 0.035
Css 1.5
CtauL 4360
Cw1 0.44
Cw2 0.92
Cw3 0.3
CwR 1.5
Clambda 2.495
CmuStd 0.09
Prtheta 0.85
Sigmak 1
Sigmaw 1.17
}
\endverbatim
SourceFiles
kkLOmega.C
\*---------------------------------------------------------------------------*/
#ifndef kkLOmega_H
#define kkLOmega_H
#include "RASModel.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
/*---------------------------------------------------------------------------*\
Class kkLOmega Declaration
\*---------------------------------------------------------------------------*/
class kkLOmega
:
public RASModel
{
// Private memmber functions
tmp<volScalarField> fv(const volScalarField& Ret) const;
tmp<volScalarField> fINT() const;
tmp<volScalarField> fSS(const volScalarField& omega) const;
tmp<volScalarField> Cmu(const volScalarField& S) const;
tmp<volScalarField> BetaTS(const volScalarField& Rew) const;
tmp<volScalarField> fTaul
(
const volScalarField& lambdaEff,
const volScalarField& ktL
) const;
tmp<volScalarField> alphaT
(
const volScalarField& lambdaEff,
const volScalarField& fv,
const volScalarField& ktS
) const;
tmp<volScalarField> fOmega
(
const volScalarField& lambdaEff,
const volScalarField& lambdaT
) const;
tmp<volScalarField> gammaBP(const volScalarField& omega) const;
tmp<volScalarField> gammaNAT
(
const volScalarField& ReOmega,
const volScalarField& fNatCrit
) const;
protected:
// Protected data
// Model coefficients
dimensionedScalar A0_;
dimensionedScalar As_;
dimensionedScalar Av_;
dimensionedScalar Abp_;
dimensionedScalar Anat_;
dimensionedScalar Ats_;
dimensionedScalar CbpCrit_;
dimensionedScalar Cnc_;
dimensionedScalar CnatCrit_;
dimensionedScalar Cint_;
dimensionedScalar CtsCrit_;
dimensionedScalar CrNat_;
dimensionedScalar C11_;
dimensionedScalar C12_;
dimensionedScalar CR_;
dimensionedScalar CalphaTheta_;
dimensionedScalar Css_;
dimensionedScalar CtauL_;
dimensionedScalar Cw1_;
dimensionedScalar Cw2_;
dimensionedScalar Cw3_;
dimensionedScalar CwR_;
dimensionedScalar Clambda_;
dimensionedScalar CmuStd_;
dimensionedScalar Prtheta_;
dimensionedScalar Sigmak_;
dimensionedScalar Sigmaw_;
// Fields
volScalarField kt_;
volScalarField omega_;
volScalarField kl_;
volScalarField nut_;
wallDist y_;
public:
//- Runtime type information
TypeName("kkLOmega");
// Constructors
//- Construct from components
kkLOmega
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& turbulenceModelName = turbulenceModel::typeName,
const word& modelName = typeName
);
//- Destructor
virtual ~kkLOmega()
{}
// Member Functions
//- Return the turbulence viscosity
virtual tmp<volScalarField> nut() const
{
return nut_;
}
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& alphaT) const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaT/Sigmak_ + nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff(const volScalarField& alphaT) const
{
return tmp<volScalarField>
(
new volScalarField("DomegaEff", alphaT/Sigmaw_ + nu())
);
}
//- Return the laminar kinetic energy
virtual tmp<volScalarField> kl() const
{
return kl_;
}
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const
{
return kt_;
}
//- Return the turbulence specific dissipation rate
virtual tmp<volScalarField> omega() const
{
return omega_;
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
mesh_.time().timeName(),
mesh_
),
kt_*omega_,
omega_.boundaryField().types()
)
);
}
//- Return the Reynolds stress tensor
virtual tmp<volSymmTensorField> R() const;
//- Return the effective stress tensor including the laminar stress
virtual tmp<volSymmTensorField> devReff() const;
//- Return the source term for the momentum equation
virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
//- Read RASProperties dictionary
virtual bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

1
tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/Allrun
View File

@ -1,7 +1,6 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

2
tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/changeDictionaryDict
View File

@ -27,7 +27,7 @@ dictionaryReplacement
D 700;
I 500;
length 1.05;
jump uniform 0
jump uniform 0;
value uniform 0;
}
porous_half1

12
tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/mechanicalProperties
View File

@ -17,20 +17,20 @@ FoamFile
rho
{
rho rhoInf;
rhoInf rhoInf [ 1 -3 0 0 0 0 0 ] 7854;
type uniform;
value 7854;
}
nu
{
nu nuInf;
nuInf nuInf [ 0 0 0 0 0 0 0 ] 0.3;
type uniform;
value 0.3;
}
E
{
E EInf;
EInf EInf [ 1 -1 -2 0 0 0 0 ] 2e+11;
type uniform;
value 2e+11;
}
planeStress yes;

12
tutorials/stressAnalysis/solidDisplacementFoam/plateHole/constant/thermalProperties
View File

@ -17,20 +17,20 @@ FoamFile
C
{
C CInf;
CInf CInf [ 0 2 -2 -1 0 0 0 ] 434;
type uniform;
value 434;
}
k
{
k kInf;
kInf kInf [ 1 1 -3 -1 0 0 0 ] 60.5;
type uniform;
value 60.5;
}
alpha
{
alpha alphaInf;
alphaInf alphaInf [ 0 0 0 -1 0 0 0 ] 1.1e-05;
type uniform;
value 1.1e-05;
}
thermalStress no;

19
tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/mechanicalProperties
View File

@ -15,13 +15,24 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
rho rho [ 1 -3 0 0 0 0 0 ] 7854;
rho
{
type uniform;
value 7854;
}
E E [ 1 -1 -2 0 0 0 0 ] 2e+11;
nu
{
type uniform;
value 0.0;
}
nu nu [ 0 0 0 0 0 0 0 ] 0;
E
{
type uniform;
value 2e+11;
}
planeStress yes;
// ************************************************************************* //

18
tutorials/stressAnalysis/solidEquilibriumDisplacementFoam/beamEndLoad/constant/thermalProperties
View File

@ -15,11 +15,23 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
C C [ 0 2 -2 -1 0 0 0 ] 434;
C
{
type uniform;
value 434;
}
k k [ 1 1 -3 -1 0 0 0 ] 60.5;
k
{
type uniform;
value 60.5;
}
alpha alpha [ 0 0 0 -1 0 0 0 ] 1.1e-05;
alpha
{
type uniform;
value 1.1e-05;
}
thermalStress no;

4
wmake/rules/linux64Gcc/c++
View File

@ -17,5 +17,5 @@ cpptoo = $(Ctoo)
LINK_LIBS = $(c++DBUG)
LINKLIBSO = $(CC) $(c++FLAGS) -shared
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed
LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed

4
wmake/rules/linux64Gcc45/c++
View File

@ -17,5 +17,5 @@ cpptoo = $(Ctoo)
LINK_LIBS = $(c++DBUG)
LINKLIBSO = $(CC) $(c++FLAGS) -shared
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed
LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed

4
wmake/rules/linux64Gcc46/c++
View File

@ -17,5 +17,5 @@ cpptoo = $(Ctoo)
LINK_LIBS = $(c++DBUG)
LINKLIBSO = $(CC) $(c++FLAGS) -shared
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed
LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed

4
wmake/rules/linuxGcc/c++
View File

@ -17,5 +17,5 @@ cpptoo = $(Ctoo)
LINK_LIBS = $(c++DBUG)
LINKLIBSO = $(CC) $(c++FLAGS) -shared
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed
LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed

4
wmake/rules/linuxGcc45/c++
View File

@ -17,5 +17,5 @@ cpptoo = $(Ctoo)
LINK_LIBS = $(c++DBUG)
LINKLIBSO = $(CC) $(c++FLAGS) -shared
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed
LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed

4
wmake/rules/linuxGcc46/c++
View File

@ -17,5 +17,5 @@ cpptoo = $(Ctoo)
LINK_LIBS = $(c++DBUG)
LINKLIBSO = $(CC) $(c++FLAGS) -shared
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed
LINKLIBSO = $(CC) $(c++FLAGS) -shared -Xlinker --add-needed -Xlinker --no-as-needed
LINKEXE = $(CC) $(c++FLAGS) -Xlinker --add-needed -Xlinker --no-as-needed