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 commit 'OpenCFD/master' into olesenm

Browse Source
This commit is contained in:
Mark Olesen
2008-11-05 10:42:05 +01:00
parent 5d0b3348c5 593db180ce
commit 345df1e970
33 changed files with 4079 additions and 333 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/OpenFOAM/containers/Lists/SortableList/SortableList.C
View File

@ -32,7 +32,7 @@ License
// Construct from List
template <class Type>
Foam::SortableList<Type>::SortableList(const List<Type>& values)
Foam::SortableList<Type>::SortableList(const UList<Type>& values)
:
List<Type>(values),
indices_(values.size())

2
src/OpenFOAM/containers/Lists/SortableList/SortableList.H
View File

@ -88,7 +88,7 @@ public:
//- Construct from List, sorting the elements.
// Starts with indices set to index in argument
explicit SortableList(const List<Type>&);
explicit SortableList(const UList<Type>&);
//- Construct from tranferred List, sorting the elements.
// Starts with indices set to index in argument

14
src/OpenFOAM/fields/Fields/Field/Field.C
View File

@ -575,20 +575,6 @@ void Field<Type>::replace
}
template<class Type>
void Field<Type>::transfer(Field<Type>& f)
{
List<Type>::transfer(f);
}
template<class Type>
void Field<Type>::transfer(List<Type>& lst)
{
List<Type>::transfer(lst);
}
template<class Type>
tmp<Field<Type> > Field<Type>::T() const
{

6
src/OpenFOAM/fields/Fields/Field/Field.H
View File

@ -297,12 +297,6 @@ public:
//- Replace a component field of the field
void replace(const direction, const cmptType&);
//- Transfer the contents of the argument Field into this Field
void transfer(Field<Type>&);
//- Transfer the contents of the argument List into this Field
void transfer(List<Type>&);
//- Return the field transpose (only defined for second rank tensors)
tmp<Field<Type> > T() const;

5
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointFaces.C
View File

@ -38,6 +38,11 @@ const labelListList& primitiveMesh::pointFaces() const
{
if (!pfPtr_)
{
if (debug)
{
Pout<< "primitiveMesh::pointFaces() : "
<< "calculating pointFaces" << endl;
}
// Invert faces()
pfPtr_ = new labelListList(nPoints());
invertManyToMany(nPoints(), faces(), *pfPtr_);

82
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
View File

@ -218,7 +218,8 @@ Foam::label Foam::autoSnapDriver::getCollocatedPoints
// Calculate displacement as average of patch points.
Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
(
const motionSmoother& meshMover
const motionSmoother& meshMover,
const List<labelPair>& baffles
) const
{
const indirectPrimitivePatch& pp = meshMover.patch();
@ -253,6 +254,34 @@ Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
const pointField& points = pp.points();
const polyMesh& mesh = meshMover.mesh();
// Get labels of faces to count (master of coupled faces and baffle pairs)
PackedList<1> isMasterFace(syncTools::getMasterFaces(mesh));
{
forAll(baffles, i)
{
label f0 = baffles[i].first();
label f1 = baffles[i].second();
if (isMasterFace.get(f0) == 1)
{
// Make f1 a slave
isMasterFace.set(f1, 0);
}
else if (isMasterFace.get(f1) == 1)
{
isMasterFace.set(f0, 0);
}
else
{
FatalErrorIn("autoSnapDriver::smoothPatchDisplacement(..)")
<< "Both sides of baffle consisting of faces " << f0
<< " and " << f1 << " are already slave faces."
<< abort(FatalError);
}
}
}
// Get average position of boundary face centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -266,9 +295,14 @@ Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
forAll(pFaces, pfI)
{
avgBoundary[patchPointI] += pp[pFaces[pfI]].centre(points);
label faceI = pFaces[pfI];
if (isMasterFace.get(pp.addressing()[faceI]) == 1)
{
avgBoundary[patchPointI] += pp[faceI].centre(points);
nBoundary[patchPointI]++;
}
}
nBoundary[patchPointI] = pFaces.size();
}
syncTools::syncPointList
@ -886,7 +920,7 @@ void Foam::autoSnapDriver::preSmoothPatch
checkFaces[faceI] = faceI;
}
pointField patchDisp(smoothPatchDisplacement(meshMover));
pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
// The current mesh is the starting mesh to smooth from.
meshMover.setDisplacement(patchDisp);
@ -1008,9 +1042,11 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
// Displacement per patch point
vectorField patchDisp(localPoints.size(), vector::zero);
if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
{
// Current surface snapped to
labelList snapSurf(localPoints.size(), -1);
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces;
labelList unzonedSurfaces;
@ -1039,16 +1075,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
patchDisp[pointI] =
hitInfo[pointI].hitPoint()
- localPoints[pointI];
snapSurf[pointI] = hitSurface[pointI];
}
//else
//{
// WarningIn("autoSnapDriver::calcNearestSurface(..)")
// << "For point:" << pointI
// << " coordinate:" << localPoints[pointI]
// << " did not find any surface within:"
// << 4*snapDist[pointI]
// << " meter." << endl;
//}
}
}
@ -1060,6 +1089,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
// Surfaces with zone information
const wordList& faceZoneNames = surfaces.faceZoneNames();
// Current best snap distance
scalarField minSnapDist(snapDist);
forAll(zonedSurfaces, i)
@ -1105,19 +1135,25 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
minSnapDist[pointI],
mag(patchDisp[pointI])
);
}
else
{
WarningIn("autoSnapDriver::calcNearestSurface(..)")
<< "For point:" << pointI
<< " coordinate:" << localPoints[pointI]
<< " did not find any surface within:"
<< 4*minSnapDist[pointI]
<< " meter." << endl;
snapSurf[pointI] = zoneSurfI;
}
}
}
// Check if all points are being snapped
forAll(snapSurf, pointI)
{
if (snapSurf[pointI] == -1)
{
WarningIn("autoSnapDriver::calcNearestSurface(..)")
<< "For point:" << pointI
<< " coordinate:" << localPoints[pointI]
<< " did not find any surface within:"
<< minSnapDist[pointI]
<< " meter." << endl;
}
}
{
scalarField magDisp(mag(patchDisp));

6
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
View File

@ -100,7 +100,11 @@ class autoSnapDriver
//- Calculate displacement per patch point to smooth out patch.
// Quite complicated in determining which points to move where.
pointField smoothPatchDisplacement(const motionSmoother&) const;
pointField smoothPatchDisplacement
(
const motionSmoother&,
const List<labelPair>&
) const;
//- Check that face zones are synced
void checkCoupledFaceZones() const;

4
src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -382,10 +382,12 @@ private:
//- Finds zone per cell for cells inside closed named surfaces.
// (uses geometric test for insideness)
// Adapts namedSurfaceIndex so all faces on boundary of cellZone
// have corresponding faceZone.
void findCellZoneGeometric
(
const labelList& closedNamedSurfaces,
const labelList& namedSurfaceIndex,
labelList& namedSurfaceIndex,
const labelList& surfaceToCellZone,
labelList& cellToZone
) const;

196
src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View File

@ -47,6 +47,7 @@ License
#include "motionSmoother.H"
#include "polyMeshGeometry.H"
#include "IOmanip.H"
#include "cellSet.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -1490,7 +1491,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
void Foam::meshRefinement::findCellZoneGeometric
(
const labelList& closedNamedSurfaces, // indices of closed surfaces
const labelList& namedSurfaceIndex, // per face index of named surface
labelList& namedSurfaceIndex, // per face index of named surface
const labelList& surfaceToCellZone, // cell zone index per surface
labelList& cellToZone
@ -1627,6 +1628,76 @@ void Foam::meshRefinement::findCellZoneGeometric
}
}
}
// Adapt the namedSurfaceIndex
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
// for if any cells were not completely covered.
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
{
// Give face the zone of the owner
namedSurfaceIndex[faceI] = findIndex
(
surfaceToCellZone,
max(ownZone, neiZone)
);
}
}
labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
for
(
label faceI = mesh_.nInternalFaces();
faceI < mesh_.nFaces();
faceI++
)
{
label own = mesh_.faceOwner()[faceI];
neiCellZone[faceI-mesh_.nInternalFaces()] = cellToZone[own];
}
syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
forAll(pp, i)
{
label faceI = pp.start()+i;
label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
{
// Give face the zone of the owner
namedSurfaceIndex[faceI] = findIndex
(
surfaceToCellZone,
max(ownZone, neiZone)
);
}
}
}
}
// Sync
syncTools::syncFaceList
(
mesh_,
namedSurfaceIndex,
maxEqOp<label>(),
false
);
}
@ -1656,7 +1727,7 @@ void Foam::meshRefinement::findCellZoneTopo
blockedFace[faceI] = true;
}
}
syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>(), false);
// No need to sync since namedSurfaceIndex already is synced
// Set region per cell based on walking
regionSplit cellRegion(mesh_, blockedFace);
@ -2149,8 +2220,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
{
FatalErrorIn
(
"meshRefinement::findCellZoneTopo"
"(const point&, const labelList&, const labelList&, labelList&)"
"meshRefinement::splitMesh"
"(const label, const labelList&, const point&)"
) << "Point " << keepPoint
<< " is not inside the mesh." << nl
<< "Bounding box of the mesh:" << mesh_.globalData().bb()
@ -2703,6 +2774,79 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
}
// Put the cells into the correct zone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Closed surfaces with cellZone specified.
labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
// Zone per cell:
// -2 : unset
// -1 : not in any zone
// >=0: zoneID
labelList cellToZone(mesh_.nCells(), -2);
// Set using geometric test
// ~~~~~~~~~~~~~~~~~~~~~~~~
if (closedNamedSurfaces.size() > 0)
{
findCellZoneGeometric
(
closedNamedSurfaces, // indices of closed surfaces
namedSurfaceIndex, // per face index of named surface
surfaceToCellZone, // cell zone index per surface
cellToZone
);
}
//{
// Pout<< "** finding out blocked faces." << endl;
//
// cellSet zonedCellsGeom(mesh_, "zonedCellsGeom", 100);
// forAll(cellToZone, cellI)
// {
// if (cellToZone[cellI] >= 0)
// {
// zonedCellsGeom.insert(cellI);
// }
// }
// Pout<< "Writing zoned cells to " << zonedCellsGeom.objectPath()
// << endl;
// zonedCellsGeom.write();
//
//
// faceSet zonedFaces(mesh_, "zonedFaces", 100);
// forAll(namedSurfaceIndex, faceI)
// {
// label surfI = namedSurfaceIndex[faceI];
//
// if (surfI != -1)
// {
// zonedFaces.insert(faceI);
// }
// }
// Pout<< "Writing zoned faces to " << zonedFaces.objectPath() << endl;
// zonedFaces.write();
//}
// Set using walking
// ~~~~~~~~~~~~~~~~~
//if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
{
// Topological walk
findCellZoneTopo
(
keepPoint,
namedSurfaceIndex,
surfaceToCellZone,
cellToZone
);
}
// Topochange container
polyTopoChange meshMod(mesh_);
@ -2770,50 +2914,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Put the cells into the correct zone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Closed surfaces with cellZone specified.
labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
// Zone per cell:
// -2 : unset
// -1 : not in any zone
// >=0: zoneID
labelList cellToZone(mesh_.nCells(), -2);
// Set using geometric test
// ~~~~~~~~~~~~~~~~~~~~~~~~
if (closedNamedSurfaces.size() > 0)
{
findCellZoneGeometric
(
closedNamedSurfaces, // indices of closed surfaces
namedSurfaceIndex, // per face index of named surface
surfaceToCellZone, // cell zone index per surface
cellToZone
);
}
// Set using walking
// ~~~~~~~~~~~~~~~~~
//if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
{
// Topological walk
findCellZoneTopo
(
keepPoint,
namedSurfaceIndex,
surfaceToCellZone,
cellToZone
);
}
// Actually move the cells to their zone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(cellToZone, cellI)
{
label zoneI = cellToZone[cellI];

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

@ -33,15 +33,68 @@ void Foam::setRefCell
const volScalarField& field,
const dictionary& dict,
label& refCelli,
scalar& refValue
scalar& refValue,
bool forceReference
)
{
if (field.needReference())
if (field.needReference() || forceReference)
{
word refCellName = field.name() + "RefCell";
word refPointName = field.name() + "RefPoint";
word refValueName = field.name() + "RefValue";
refCelli = readLabel(dict.lookup(refCellName));
if (dict.found(refCellName))
{
if (Pstream::master())
{
refCelli = readLabel(dict.lookup(refCellName));
}
else
{
refCelli = -1;
}
}
else if (dict.found(refPointName))
{
point refPointi(dict.lookup(refPointName));
refCelli = field.mesh().findCell(refPointi);
label hasRef = (refCelli >= 0 ? 1 : 0);
label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());
if (sumHasRef != 1)
{
FatalErrorIn
(
"void Foam::setRefCell"
"("
" const volScalarField&,"
" const dictionary&,"
" label& scalar&,"
" bool"
")"
)
<< "Unable to set reference cell for field " << field.name()
<< nl << " Reference point " << refPointName
<< " found on multiple domains" << nl << abort(FatalError);
}
}
else
{
FatalErrorIn
(
"void Foam::setRefCell"
"("
" const volScalarField&,"
" const dictionary&,"
" label& scalar&,"
" bool"
")"
)
<< "Unable to set reference cell for field" << field.name() << nl
<< " Please supply either " << refCellName
<< " or " << refPointName << nl << abort(FatalError);
}
refValue = readScalar(dict.lookup(refValueName));
}
}

3
src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
View File

@ -52,7 +52,8 @@ void setRefCell
const volScalarField& field,
const dictionary& dict,
label& refCelli,
scalar& refValue
scalar& refValue,
bool forceReference = false
);
}

10
src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
View File

@ -234,6 +234,16 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
(
mpp.samplePatch()
);
if (patchID < 0)
{
FatalErrorIn
(
"void directMappedFixedValueFvPatchField<Type>::"
"updateCoeffs()"
)<< "Unable to find sample patch " << mpp.samplePatch()
<< " for patch " << this->patch().name() << nl
<< abort(FatalError);
}
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
const word& fieldName = this->dimensionedInternalField().name();
const fieldType& sendField =

2
src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
View File

@ -483,7 +483,7 @@ void Foam::fvMatrix<Type>::setReference
{
if (psi_.needReference() || forceReference)
{
if (Pstream::master())
if (cell >= 0)
{
source()[cell] += diag()[cell]*value;
diag()[cell] += diag()[cell];

1
src/meshTools/Make/files
View File

@ -107,6 +107,7 @@ $(pointSources)/faceToPoint/faceToPoint.C
$(pointSources)/boxToPoint/boxToPoint.C
$(pointSources)/surfaceToPoint/surfaceToPoint.C
$(pointSources)/zoneToPoint/zoneToPoint.C
$(pointSources)/nearestToPoint/nearestToPoint.C
surfaceSets/surfaceSets.C

51
src/meshTools/searchableSurface/searchableSurfaceWithGaps.C
View File

@ -26,7 +26,6 @@ License
#include "searchableSurfaceWithGaps.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
#include "Time.H"
#include "ListOps.H"
@ -82,7 +81,7 @@ Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
// Do second offset vector perp to original edge and first offset vector
offsets[1] = n ^ offsets[0];
offsets[1] *= gap_/mag(offsets[1]);
offsets[1] *= gap_;
}
return offsets;
@ -207,6 +206,10 @@ void Foam::searchableSurfaceWithGaps::findLine
List<pointIndexHit>& info
) const
{
// Test with unperturbed vectors
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
surface().findLine(start, end, info);
// Count number of misses. Determine map
@ -215,6 +218,10 @@ void Foam::searchableSurfaceWithGaps::findLine
if (returnReduce(nMiss, sumOp<label>()) > 0)
{
//Pout<< "** retesting with offset0 " << nMiss << " misses out of "
// << start.size() << endl;
// extract segments according to map
pointField compactStart(start, compactMap);
pointField compactEnd(end, compactMap);
@ -228,20 +235,36 @@ void Foam::searchableSurfaceWithGaps::findLine
offset1
);
// Test with offset0 perturbed vectors
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// test in pairs: only if both perturbations hit something
// do we accept the hit.
const vectorField smallVec(SMALL*(compactEnd-compactStart));
List<pointIndexHit> plusInfo;
surface().findLine(compactStart+offset0, compactEnd+offset0, plusInfo);
surface().findLine
(
compactStart+offset0-smallVec,
compactEnd+offset0+smallVec,
plusInfo
);
List<pointIndexHit> minInfo;
surface().findLine(compactStart-offset0, compactEnd-offset0, minInfo);
surface().findLine
(
compactStart-offset0-smallVec,
compactEnd-offset0+smallVec,
minInfo
);
// Extract any hits
forAll(plusInfo, i)
{
if (plusInfo[i].hit() && minInfo[i].hit())
{
info[compactMap[i]] = plusInfo[i].hitPoint()-offset0[i];
info[compactMap[i]] = plusInfo[i];
info[compactMap[i]].rawPoint() -= offset0[i];
}
}
@ -250,6 +273,12 @@ void Foam::searchableSurfaceWithGaps::findLine
if (returnReduce(nMiss, sumOp<label>()) > 0)
{
//Pout<< "** retesting with offset1 " << nMiss << " misses out of "
// << start.size() << endl;
// Test with offset1 perturbed vectors
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Extract (inplace possible because of order)
forAll(plusMissMap, i)
{
@ -266,17 +295,18 @@ void Foam::searchableSurfaceWithGaps::findLine
offset0.setSize(plusMissMap.size());
offset1.setSize(plusMissMap.size());
const vectorField smallVec(SMALL*(compactEnd-compactStart));
surface().findLine
(
compactStart+offset1,
compactEnd+offset1,
compactStart+offset1-smallVec,
compactEnd+offset1+smallVec,
plusInfo
);
surface().findLine
(
compactStart-offset1,
compactEnd-offset1,
compactStart-offset1-smallVec,
compactEnd-offset1+smallVec,
minInfo
);
@ -285,7 +315,8 @@ void Foam::searchableSurfaceWithGaps::findLine
{
if (plusInfo[i].hit() && minInfo[i].hit())
{
info[compactMap[i]] = plusInfo[i].hitPoint()-offset1[i];
info[compactMap[i]] = plusInfo[i];
info[compactMap[i]].rawPoint() -= offset1[i];
}
}
}

8
src/meshTools/searchableSurface/searchableSurfaceWithGaps.H
View File

@ -27,7 +27,13 @@ Class
Description
searchableSurface using multiple slightly shifted underlying surfaces
to make sure pierces don't go through gaps.
to make sure pierces don't go through gaps:
- shift test vector with two small vectors (of size gap_) perpendicular
to the original.
Test with + and - this vector. Only if both register a hit is it seen
as one.
- extend the test vector slightly (with SMALL) to account for numerical
inaccuracies.
SourceFiles
searchableSurfaceWithGaps.C

154
src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.C Normal file
View File

@ -0,0 +1,154 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "nearestToPoint.H"
#include "polyMesh.H"
#include "meshSearch.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(nearestToPoint, 0);
addToRunTimeSelectionTable(topoSetSource, nearestToPoint, word);
addToRunTimeSelectionTable(topoSetSource, nearestToPoint, istream);
}
Foam::topoSetSource::addToUsageTable Foam::nearestToPoint::usage_
(
nearestToPoint::typeName,
"\n Usage: nearestToPoint (pt0 .. ptn)\n\n"
" Select the nearest point for each of the points pt0 ..ptn\n\n"
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::nearestToPoint::combine(topoSet& set, const bool add) const
{
// Do linear search since usually just a few points.
forAll(points_, pointI)
{
const pointField& pts = mesh_.points();
if (pts.size() > 0)
{
label minPointI = 0;
scalar minDistSqr = magSqr(pts[minPointI] - points_[pointI]);
for (label i = 1; i < pts.size(); i++)
{
scalar distSqr = magSqr(pts[i] - points_[pointI]);
if (distSqr < minDistSqr)
{
minDistSqr = distSqr;
minPointI = i;
}
}
addOrDelete(set, minPointI, add);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::nearestToPoint::nearestToPoint
(
const polyMesh& mesh,
const pointField& points
)
:
topoSetSource(mesh),
points_(points)
{}
// Construct from dictionary
Foam::nearestToPoint::nearestToPoint
(
const polyMesh& mesh,
const dictionary& dict
)
:
topoSetSource(mesh),
points_(dict.lookup("points"))
{}
// Construct from Istream
Foam::nearestToPoint::nearestToPoint
(
const polyMesh& mesh,
Istream& is
)
:
topoSetSource(mesh),
points_(checkIs(is))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::nearestToPoint::~nearestToPoint()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::nearestToPoint::applyToSet
(
const topoSetSource::setAction action,
topoSet& set
) const
{
if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
{
Info<< " Adding points nearest to " << points_ << endl;
combine(set, true);
}
else if (action == topoSetSource::DELETE)
{
Info<< " Removing points nearest to " << points_ << endl;
combine(set, false);
}
}
// ************************************************************************* //

122
src/meshTools/sets/pointSources/nearestToPoint/nearestToPoint.H Normal file
View File

@ -0,0 +1,122 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::nearestToPoint
Description
A topoSetSource to select points nearest to points.
SourceFiles
nearestToPoint.C
\*---------------------------------------------------------------------------*/
#ifndef nearestToPoint_H
#define nearestToPoint_H
#include "topoSetSource.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class nearestToPoint Declaration
\*---------------------------------------------------------------------------*/
class nearestToPoint
:
public topoSetSource
{
// Private data
//- Add usage string
static addToUsageTable usage_;
//- points to select nearest to
pointField points_;
// Private Member Functions
void combine(topoSet& set, const bool add) const;
public:
//- Runtime type information
TypeName("nearestToPoint");
// Constructors
//- Construct from components
nearestToPoint
(
const polyMesh& mesh,
const pointField& points
);
//- Construct from dictionary
nearestToPoint
(
const polyMesh& mesh,
const dictionary& dict
);
//- Construct from Istream
nearestToPoint
(
const polyMesh& mesh,
Istream&
);
// Destructor
virtual ~nearestToPoint();
// Member Functions
virtual void applyToSet
(
const topoSetSource::setAction action,
topoSet&
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

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

@ -164,7 +164,7 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::hc() const
(
(((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd
+ a[0]
)*Tstd
)*Tstd + a[5]
);
}

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

@ -27,6 +27,7 @@ dynMixedSmagorinsky/dynMixedSmagorinsky.C
/*Smagorinsky2/Smagorinsky2.C*/
kOmegaSSTSAS/kOmegaSSTSAS.C
/* Wall functions */
wallFunctions=derivedFvPatchFields/wallFunctions

447
src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.C Normal file
View File

@ -0,0 +1,447 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "kOmegaSSTSAS.H"
#include "addToRunTimeSelectionTable.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace LESModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(kOmegaSSTSAS, 0);
addToRunTimeSelectionTable(LESModel, kOmegaSSTSAS, dictionary);
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<volScalarField> kOmegaSSTSAS::F1(const volScalarField& CDkOmega) const
{
volScalarField CDkOmegaPlus = max
(
CDkOmega,
dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
);
volScalarField arg1 = min
(
min
(
max
(
(scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*nu()/(sqr(y_)*omega_)
),
(4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
),
scalar(10)
);
return tanh(pow4(arg1));
}
tmp<volScalarField> kOmegaSSTSAS::F2() const
{
volScalarField arg2 = min
(
max
(
(scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*nu()/(sqr(y_)*omega_)
),
scalar(100)
);
return tanh(sqr(arg2));
}
tmp<volScalarField> kOmegaSSTSAS::Lvk2
(
const volScalarField& S2
) const
{
return kappa_*sqrt(S2)
/(
mag(fvc::laplacian(U()))
+ dimensionedScalar
(
"ROOTVSMALL",
dimensionSet(0, -1 , -1, 0, 0, 0, 0),
ROOTVSMALL
)
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
kOmegaSSTSAS::kOmegaSSTSAS
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& modelName
)
:
LESModel(modelName, U, phi, transport),
alphaK1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK1",
coeffDict(),
0.85034
)
),
alphaK2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK2",
coeffDict(),
1.0
)
),
alphaOmega1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega1",
coeffDict(),
0.5
)
),
alphaOmega2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega2",
coeffDict(),
0.85616
)
),
gamma1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma1",
coeffDict(),
0.5532
)
),
gamma2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma2",
coeffDict(),
0.4403
)
),
beta1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta1",
coeffDict(),
0.075
)
),
beta2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta2",
coeffDict(),
0.0828
)
),
betaStar_
(
dimensioned<scalar>::lookupOrAddToDict
(
"betaStar",
coeffDict(),
0.09
)
),
a1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"a1",
coeffDict(),
0.31
)
),
c1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"c1",
coeffDict(),
10.0
)
),
alphaPhi_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaPhi",
coeffDict(),
0.666667
)
),
zetaTilda2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"zetaTilda2",
coeffDict(),
1.755
)
),
FSAS_
(
dimensioned<scalar>::lookupOrAddToDict
(
"FSAS",
coeffDict(),
1.25
)
),
omega0_("omega0", dimless/dimTime, SMALL),
omegaSmall_("omegaSmall", dimless/dimTime, SMALL),
y_(mesh_),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
coeffDict(),
0.09
)
),
kappa_
(
dimensioned<scalar>::lookupOrAddToDict
(
"kappa",
*this,
0.4187
)
),
k_
(
IOobject
(
"k",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
omega_
(
IOobject
(
"omega",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
nuSgs_
(
IOobject
(
"nuSgs",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
)
{
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void kOmegaSSTSAS::correct(const tmp<volTensorField>& gradU)
{
LESModel::correct(gradU);
if (mesh_.changing())
{
y_.correct();
}
volScalarField S2 = magSqr(symm(gradU()));
gradU.clear();
volVectorField gradK = fvc::grad(k_);
volVectorField gradOmega = fvc::grad(omega_);
volScalarField L = sqrt(k_)/(pow(Cmu_, 0.25)*(omega_ + omegaSmall_));
volScalarField CDkOmega = (2.0*alphaOmega2_)*(gradK & gradOmega)/(omega_ + omegaSmall_);
volScalarField F1 = this->F1(CDkOmega);
volScalarField G = nuSgs_*2.0*S2;
// Turbulent kinetic energy equation
solve
(
fvm::ddt(k_)
+ fvm::div(phi(), k_)
- fvm::Sp(fvc::div(phi()), k_)
- fvm::laplacian(DkEff(F1), k_)
==
min(G, c1_*betaStar_*k_*omega_)
- fvm::Sp(betaStar_*omega_, k_)
);
bound(k_, k0());
volScalarField grad_omega_k = max
(
magSqr(gradOmega)/
sqr(omega_ + omegaSmall_),
magSqr(gradK)/
sqr(k_ + k0())
);
// Turbulent frequency equation
solve
(
fvm::ddt(omega_)
+ fvm::div(phi(), omega_)
- fvm::Sp(fvc::div(phi()), omega_)
- fvm::laplacian(DomegaEff(F1), omega_)
==
gamma(F1)*2.0*S2
- fvm::Sp(beta(F1)*omega_, omega_)
- fvm::SuSp // cross diffusion term
(
(F1 - scalar(1))*CDkOmega/omega_,
omega_
)
+ FSAS_
*max
(
dimensionedScalar("zero",dimensionSet(0, 0 , -2, 0, 0),0. ),
zetaTilda2_*kappa_*S2*(L/Lvk2(S2))- 2.0/alphaPhi_*k_*grad_omega_k
)
);
bound(omega_, omega0_);
// Re-calculate viscosity
nuSgs_ == a1_*k_/max(a1_*omega_, F2()*sqrt(S2));
nuSgs_.correctBoundaryConditions();
}
tmp<volScalarField> kOmegaSSTSAS::epsilon() const
{
return 2.0*nuEff()*magSqr(symm(fvc::grad(U())));
}
tmp<volSymmTensorField> kOmegaSSTSAS::B() const
{
return ((2.0/3.0)*I)*k() - nuSgs()*twoSymm(fvc::grad(U()));
}
tmp<volSymmTensorField> kOmegaSSTSAS::devBeff() const
{
return -nuEff()*dev(twoSymm(fvc::grad(U())));
}
tmp<fvVectorMatrix> kOmegaSSTSAS::divDevBeff(volVectorField& U) const
{
return
(
- fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(fvc::grad(U)().T()))
);
}
bool kOmegaSSTSAS::read()
{
if (LESModel::read())
{
alphaK1_.readIfPresent(coeffDict());
alphaK2_.readIfPresent(coeffDict());
alphaOmega1_.readIfPresent(coeffDict());
alphaOmega2_.readIfPresent(coeffDict());
gamma1_.readIfPresent(coeffDict());
gamma2_.readIfPresent(coeffDict());
beta1_.readIfPresent(coeffDict());
beta2_.readIfPresent(coeffDict());
betaStar_.readIfPresent(coeffDict());
a1_.readIfPresent(coeffDict());
c1_.readIfPresent(coeffDict());
alphaPhi_.readIfPresent(coeffDict());
zetaTilda2_.readIfPresent(coeffDict());
FSAS_.readIfPresent(coeffDict());
return true;
}
else
{
return false;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace LESModels
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

258
src/turbulenceModels/LES/incompressible/kOmegaSSTSAS/kOmegaSSTSAS.H Normal file
View File

@ -0,0 +1,258 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::LESmodels::kOmegaSSTSAS
Description
kOmegaSSTSAS LES turbulence model for incompressible flows
Note: does not have an explicit dependency on spatial discretisation
i.e. LESdelta not used.
SourceFiles
kOmegaSSTSAS.C
\*---------------------------------------------------------------------------*/
#ifndef kOmegaSSTSAS_H
#define kOmegaSSTSAS_H
#include "LESModel.H"
#include "volFields.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace LESModels
{
/*---------------------------------------------------------------------------*\
Class kOmegaSSTSAS Declaration
\*---------------------------------------------------------------------------*/
class kOmegaSSTSAS
:
public LESModel
{
// Private member functions
// Disallow default bitwise copy construct and assignment
kOmegaSSTSAS(const kOmegaSSTSAS&);
kOmegaSSTSAS& operator=(const kOmegaSSTSAS&);
protected:
// Protected data
// Model constants
dimensionedScalar alphaK1_;
dimensionedScalar alphaK2_;
dimensionedScalar alphaOmega1_;
dimensionedScalar alphaOmega2_;
dimensionedScalar gamma1_;
dimensionedScalar gamma2_;
dimensionedScalar beta1_;
dimensionedScalar beta2_;
dimensionedScalar betaStar_;
dimensionedScalar a1_;
dimensionedScalar c1_;
dimensionedScalar alphaPhi_;
dimensionedScalar zetaTilda2_;
dimensionedScalar FSAS_;
dimensionedScalar omega0_;
dimensionedScalar omegaSmall_;
wallDist y_;
dimensionedScalar Cmu_;
dimensionedScalar kappa_;
// Fields
volScalarField k_;
volScalarField omega_;
volScalarField nuSgs_;
// Protected member functions
tmp<volScalarField> Lvk2
(
const volScalarField& S2
) const;
tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
tmp<volScalarField> F2() const;
tmp<volScalarField> blend
(
const volScalarField& F1,
const dimensionedScalar& psi1,
const dimensionedScalar& psi2
) const
{
return F1*(psi1 - psi2) + psi2;
}
tmp<volScalarField> alphaK
(
const volScalarField& F1
) const
{
return blend(F1, alphaK1_, alphaK2_);
}
tmp<volScalarField> alphaOmega
(
const volScalarField& F1
) const
{
return blend(F1, alphaOmega1_, alphaOmega2_);
}
tmp<volScalarField> beta
(
const volScalarField& F1
) const
{
return blend(F1, beta1_, beta2_);
}
tmp<volScalarField> gamma
(
const volScalarField& F1
) const
{
return blend(F1, gamma1_, gamma2_);
}
public:
//- Runtime type information
TypeName("kOmegaSSTSAS");
// Constructors
//- Constructor from components
kOmegaSSTSAS
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport,
const word& modelName = typeName
);
//- Destructor
virtual ~kOmegaSSTSAS()
{}
// Member Functions
//- Return SGS kinetic energy
tmp<volScalarField> k() const
{
return k_;
}
//- Return omega
tmp<volScalarField> omega() const
{
return omega_;
}
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaK(F1)*nuSgs_ + nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DomegaEff", alphaOmega(F1)*nuSgs_ + nu())
);
}
//- Return sub-grid disipation rate
tmp<volScalarField> epsilon() const;
//- Return SGS viscosity
tmp<volScalarField> nuSgs() const
{
return nuSgs_;
}
//- Return the sub-grid stress tensor.
tmp<volSymmTensorField> B() const;
//- Return the effective sub-grid turbulence stress tensor
// including the laminar stress
tmp<volSymmTensorField> devBeff() const;
//- Return the deviatoric part of the divergence of Beff
// i.e. the additional term in the filtered NSE.
tmp<fvVectorMatrix> divDevBeff(volVectorField& U) const;
//- Solve the turbulence equations (k-w) and correct the turbulence viscosity
virtual void correct(const tmp<volTensorField>& gradU);
//- Read turbulenceProperties dictionary
virtual bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace LESModels
} // End namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //