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' into dsmc

Browse Source
This commit is contained in:
graham
2009-03-02 13:11:18 +00:00
parent 5e414f1146 b037b0e98d
commit 8b4c6f0662
56 changed files with 2862 additions and 790 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/Allwmake
View File

@ -5,8 +5,8 @@ set -x
# update OpenFOAM version strings if required
wmakePrintBuild -check || /bin/rm -f OpenFOAM/Make/$WM_OPTIONS/global.? 2>/dev/null
wmakeLnInclude -f OpenFOAM
wmakeLnInclude -f OSspecific/$WM_OS
wmakeLnInclude -f OpenFOAM -sf
wmakeLnInclude -f OSspecific/$WM_OS -sf
Pstream/Allwmake
wmake libo OSspecific/$WM_OS

21
src/OpenFOAM/db/Time/findInstance.C
View File

@ -58,8 +58,9 @@ Foam::word Foam::Time::findInstance
{
if (debug)
{
Info<< "Time::findInstance(const fileName&, const word&) : "
<< "found \"" << name
Info<< "Time::findInstance"
"(const fileName&, const word&, const IOobject::readOption)"
<< " : found \"" << name
<< "\" in " << timeName()/dir
<< endl;
}
@ -98,8 +99,8 @@ Foam::word Foam::Time::findInstance
if (debug)
{
Info<< "Time::findInstance"
"(const fileName&,const word&) : "
<< "found \"" << name
"(const fileName&, const word&, const IOobject::readOption)"
<< " : found \"" << name
<< "\" in " << ts[instanceI].name()/dir
<< endl;
}
@ -129,8 +130,8 @@ Foam::word Foam::Time::findInstance
if (debug)
{
Info<< "Time::findInstance"
"(const fileName&,const word&) : "
<< "found \"" << name
"(const fileName&, const word&, const IOobject::readOption)"
<< " : found \"" << name
<< "\" in " << constant()/dir
<< endl;
}
@ -141,10 +142,10 @@ Foam::word Foam::Time::findInstance
if (rOpt == IOobject::MUST_READ)
{
FatalErrorIn
(
"Time::findInstance(const fileName&,const word&)"
)
<< "Cannot find file \"" << name << "\" in directory "
(
"Time::findInstance"
"(const fileName&, const word&, const IOobject::readOption)"
) << "Cannot find file \"" << name << "\" in directory "
<< constant()/dir
<< exit(FatalError);
}

8
src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
View File

@ -154,7 +154,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
pow
(
gsf1.dimensions(),
dimensionedScalar("1", 1.0, gsf2.dimensions())
dimensionedScalar("1", gsf2.dimensions(), 1.0)
)
)
);
@ -183,7 +183,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
pow
(
gsf1.dimensions(),
dimensionedScalar("1", 1.0, gsf2.dimensions())
dimensionedScalar("1", gsf2.dimensions(), 1.0)
)
)
);
@ -214,7 +214,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
pow
(
gsf1.dimensions(),
dimensionedScalar("1", 1.0, gsf2.dimensions())
dimensionedScalar("1", gsf2.dimensions(), 1.0)
)
)
);
@ -247,7 +247,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
pow
(
gsf1.dimensions(),
dimensionedScalar("1", 1.0, gsf2.dimensions())
dimensionedScalar("1", gsf2.dimensions(), 1.0)
)
)
);

10
src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C
View File

@ -60,7 +60,7 @@ void Foam::lduMatrix::Amul
interfaceBouCoeffs,
interfaces,
psi,
Apsi,
Apsi,
cmpt
);
@ -79,7 +79,7 @@ void Foam::lduMatrix::Amul
register const label nFaces = upper().size();
#ifdef ICC_IA64_PREFETCH
#pragma swp
#pragma swp
#endif
for (register label face=0; face<nFaces; face++)
{
@ -227,7 +227,7 @@ void Foam::lduMatrix::sumA
}
#ifdef ICC_IA64_PREFETCH
#pragma swp
#pragma swp
#endif
for (register label face=0; face<nFaces; face++)
@ -316,7 +316,7 @@ void Foam::lduMatrix::residual
mBouCoeffs,
interfaces,
psi,
rA,
rA,
cmpt
);
@ -336,7 +336,7 @@ void Foam::lduMatrix::residual
register const label nFaces = upper().size();
#ifdef ICC_IA64_PREFETCH
#pragma swp
#pragma swp
#endif
for (register label face=0; face<nFaces; face++)
{

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

@ -54,40 +54,79 @@ void Foam::polyMesh::calcDirections() const
{
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
directions_[cmpt] = 1;
solutionD_[cmpt] = 1;
}
label nEmptyPatches = 0;
// Knock out empty and wedge directions. Note:they will be present on all
// domains.
vector dirVec = vector::zero;
label nEmptyPatches = 0;
label nWedgePatches = 0;
vector emptyDirVec = vector::zero;
vector wedgeDirVec = vector::zero;
forAll(boundaryMesh(), patchi)
{
if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
if (boundaryMesh()[patchi].size())
{
if (boundaryMesh()[patchi].size())
if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
{
nEmptyPatches++;
dirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
}
else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
{
const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
(
boundaryMesh()[patchi]
);
nWedgePatches++;
wedgeDirVec += cmptMag(wpp.centreNormal());
}
}
}
if (nEmptyPatches)
{
reduce(dirVec, sumOp<vector>());
reduce(emptyDirVec, sumOp<vector>());
dirVec /= mag(dirVec);
emptyDirVec /= mag(emptyDirVec);
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (dirVec[cmpt] > 1e-6)
if (emptyDirVec[cmpt] > 1e-6)
{
directions_[cmpt] = -1;
solutionD_[cmpt] = -1;
}
else
{
directions_[cmpt] = 1;
solutionD_[cmpt] = 1;
}
}
}
// Knock out wedge directions
geometricD_ = solutionD_;
if (nWedgePatches)
{
reduce(wedgeDirVec, sumOp<vector>());
wedgeDirVec /= mag(wedgeDirVec);
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (wedgeDirVec[cmpt] > 1e-6)
{
geometricD_[cmpt] = -1;
}
else
{
geometricD_[cmpt] = 1;
}
}
}
@ -163,7 +202,8 @@ Foam::polyMesh::polyMesh(const IOobject& io)
*this
),
bounds_(points_),
directions_(Vector<label>::zero),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
pointZones_
(
IOobject
@ -350,7 +390,8 @@ Foam::polyMesh::polyMesh
0
),
bounds_(points_, syncPar),
directions_(Vector<label>::zero),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
pointZones_
(
IOobject
@ -505,7 +546,8 @@ Foam::polyMesh::polyMesh
0
),
bounds_(points_, syncPar),
directions_(Vector<label>::zero),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
pointZones_
(
IOobject
@ -766,44 +808,37 @@ const Foam::fileName& Foam::polyMesh::facesInstance() const
}
const Foam::Vector<Foam::label>& Foam::polyMesh::directions() const
const Foam::Vector<Foam::label>& Foam::polyMesh::geometricD() const
{
if (directions_.x() == 0)
if (geometricD_.x() == 0)
{
calcDirections();
}
return directions_;
return geometricD_;
}
Foam::label Foam::polyMesh::nGeometricD() const
{
label nWedges = 0;
return cmptSum(geometricD() + Vector<label>::one)/2;
}
forAll(boundary_, patchi)
const Foam::Vector<Foam::label>& Foam::polyMesh::solutionD() const
{
if (solutionD_.x() == 0)
{
if (isA<wedgePolyPatch>(boundary_[patchi]))
{
nWedges++;
}
calcDirections();
}
if (nWedges != 0 && nWedges != 2 && nWedges != 4)
{
FatalErrorIn("label polyMesh::nGeometricD() const")
<< "Number of wedge patches " << nWedges << " is incorrect, "
"should be 0, 2 or 4"
<< exit(FatalError);
}
return nSolutionD() - nWedges/2;
return solutionD_;
}
Foam::label Foam::polyMesh::nSolutionD() const
{
return cmptSum(directions() + Vector<label>::one)/2;
return cmptSum(solutionD() + Vector<label>::one)/2;
}
@ -823,6 +858,10 @@ void Foam::polyMesh::addPatches
<< abort(FatalError);
}
// Reset valid directions
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
boundary_.setSize(p.size());
// Copy the patch pointers
@ -1037,6 +1076,10 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
faceZones_.movePoints(points_);
cellZones_.movePoints(points_);
// Reset valid directions (could change with rotation)
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
// Hack until proper callbacks. Below are all the polyMeh MeshObjects with a
// movePoints function.

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

@ -120,9 +120,13 @@ private:
// Created from points on construction, updated when the mesh moves
boundBox bounds_;
//- vector of non-constrained directions in mesh
// defined according to the presence of empty and wedge patches
mutable Vector<label> geometricD_;
//- vector of valid directions in mesh
// defined according to the presence of empty patches
mutable Vector<label> directions_;
mutable Vector<label> solutionD_;
// Zoning information
@ -309,17 +313,22 @@ public:
return bounds_;
}
//- Return the vector of valid directions in mesh.
// Defined according to the presence of empty patches.
// 1 indicates valid direction and -1 an invalid direction.
const Vector<label>& directions() const;
//- Return the vector of geometric directions in mesh.
// Defined according to the presence of empty and wedge patches.
// 1 indicates unconstrained direction and -1 a constrained
// direction.
const Vector<label>& geometricD() const;
//- Return the number of valid geometric dimensions in the mesh
label nGeometricD() const;
//- Return the number of valid solution dimensions in the mesh.
// For wedge cases this includes the circumferential direction
// in case of swirl.
//- Return the vector of solved-for directions in mesh.
// Differs from geometricD in that it includes for wedge cases
// the circumferential direction in case of swirl.
// 1 indicates valid direction and -1 an invalid direction.
const Vector<label>& solutionD() const;
//- Return the number of valid solved-for dimensions in the mesh
label nSolutionD() const;
//- Return point zone mesh

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

@ -68,6 +68,10 @@ void Foam::polyMesh::clearGeom()
boundary_[patchI].clearGeom();
}
// Reset valid directions (could change with rotation)
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
pointMesh::Delete(*this);
}
@ -87,6 +91,10 @@ void Foam::polyMesh::clearAddressing()
// recalculation
deleteDemandDrivenData(globalMeshDataPtr_);
// Reset valid directions
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
pointMesh::Delete(*this);
}

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

@ -217,7 +217,8 @@ Foam::polyMesh::polyMesh
boundaryFaces.size() + 1 // add room for a default patch
),
bounds_(points_, syncPar),
directions_(Vector<label>::zero),
geometricD_(Vector<label>::zero),
solutionD_(Vector<label>::zero),
pointZones_
(
IOobject

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

@ -68,6 +68,11 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
newMotionPoints.map(oldMotionPoints, mpm.pointMap());
}
// Reset valid directions (could change by faces put into empty patches)
geometricD_ = Vector<label>::zero;
solutionD_ = Vector<label>::zero;
// Hack until proper callbacks. Below are all the polyMesh-MeshObjects.
// pointMesh

41
src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
View File

@ -29,6 +29,7 @@ License
#include "mathematicalConstants.H"
#include "refinementSurfaces.H"
#include "searchableSurfaces.H"
#include "regExp.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -300,10 +301,44 @@ Foam::layerParameters::layerParameters
// readScalar(layerDict.lookup("minThickness"));
}
}
// Check whether layer specification matches any patches
const List<keyType> wildCards = layersDict.keys(true);
forAll(wildCards, i)
{
regExp re(wildCards[i]);
bool hasMatch = false;
forAll(boundaryMesh, patchI)
{
if (re.match(boundaryMesh[patchI].name()))
{
hasMatch = true;
break;
}
}
if (!hasMatch)
{
IOWarningIn("layerParameters::layerParameters(..)", layersDict)
<< "Wildcard layer specification for " << wildCards[i]
<< " does not match any patch." << endl;
}
}
const List<keyType> nonWildCards = layersDict.keys(false);
forAll(nonWildCards, i)
{
if (boundaryMesh.findPatchID(nonWildCards[i]) == -1)
{
IOWarningIn("layerParameters::layerParameters(..)", layersDict)
<< "Layer specification for " << nonWildCards[i]
<< " does not match any patch." << endl;
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

147
src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
View File

@ -265,101 +265,112 @@ Foam::refinementSurfaces::refinementSurfaces
zoneInside_(surfacesDict.size()),
regionOffset_(surfacesDict.size())
{
labelList globalMinLevel(surfacesDict.size(), 0);
labelList globalMaxLevel(surfacesDict.size(), 0);
scalarField globalAngle(surfacesDict.size(), -GREAT);
List<Map<label> > regionMinLevel(surfacesDict.size());
List<Map<label> > regionMaxLevel(surfacesDict.size());
List<Map<scalar> > regionAngle(surfacesDict.size());
// Wilcard specification : loop over all surface, all regions
// and try to find a match.
// Count number of surfaces.
label surfI = 0;
forAllConstIter(dictionary, surfacesDict, iter)
forAll(allGeometry.names(), geomI)
{
names_[surfI] = iter().keyword();
const word& geomName = allGeometry_.names()[geomI];
surfaces_[surfI] = allGeometry_.findSurfaceID(names_[surfI]);
if (surfaces_[surfI] == -1)
if (surfacesDict.found(geomName))
{
FatalErrorIn
(
"refinementSurfaces::refinementSurfaces"
"(const searchableSurfaces&, const dictionary>&"
) << "No surface called " << iter().keyword() << endl
<< "Valid surfaces are " << allGeometry_.names()
<< exit(FatalError);
surfI++;
}
const dictionary& dict = surfacesDict.subDict(iter().keyword());
}
const labelPair refLevel(dict.lookup("level"));
globalMinLevel[surfI] = refLevel[0];
globalMaxLevel[surfI] = refLevel[1];
// Size lists
surfaces_.setSize(surfI);
names_.setSize(surfI);
faceZoneNames_.setSize(surfI);
cellZoneNames_.setSize(surfI);
zoneInside_.setSize(surfI);
regionOffset_.setSize(surfI);
// Global zone names per surface
if (dict.found("faceZone"))
labelList globalMinLevel(surfI, 0);
labelList globalMaxLevel(surfI, 0);
scalarField globalAngle(surfI, -GREAT);
List<Map<label> > regionMinLevel(surfI);
List<Map<label> > regionMaxLevel(surfI);
List<Map<scalar> > regionAngle(surfI);
surfI = 0;
forAll(allGeometry.names(), geomI)
{
const word& geomName = allGeometry_.names()[geomI];
if (surfacesDict.found(geomName))
{
dict.lookup("faceZone") >> faceZoneNames_[surfI];
dict.lookup("cellZone") >> cellZoneNames_[surfI];
dict.lookup("zoneInside") >> zoneInside_[surfI];
}
const dictionary& dict = surfacesDict.subDict(geomName);
// Global perpendicular angle
if (dict.found("perpendicularAngle"))
{
globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
}
names_[surfI] = geomName;
surfaces_[surfI] = geomI;
if (dict.found("regions"))
{
const dictionary& regionsDict = dict.subDict("regions");
const wordList& regionNames =
allGeometry_[surfaces_[surfI]].regions();
const labelPair refLevel(dict.lookup("level"));
globalMinLevel[surfI] = refLevel[0];
globalMaxLevel[surfI] = refLevel[1];
forAllConstIter(dictionary, regionsDict, iter)
// Global zone names per surface
if (dict.found("faceZone"))
{
const word& key = iter().keyword();
dict.lookup("faceZone") >> faceZoneNames_[surfI];
dict.lookup("cellZone") >> cellZoneNames_[surfI];
dict.lookup("zoneInside") >> zoneInside_[surfI];
}
if (regionsDict.isDict(key))
// Global perpendicular angle
if (dict.found("perpendicularAngle"))
{
globalAngle[surfI] = readScalar
(
dict.lookup("perpendicularAngle")
);
}
if (dict.found("regions"))
{
const dictionary& regionsDict = dict.subDict("regions");
const wordList& regionNames =
allGeometry_[surfaces_[surfI]].regions();
forAll(regionNames, regionI)
{
// Get the dictionary for region iter.keyword()
const dictionary& regionDict = regionsDict.subDict(key);
label regionI = findIndex(regionNames, key);
if (regionI == -1)
if (regionsDict.found(regionNames[regionI]))
{
FatalErrorIn
// Get the dictionary for region
const dictionary& regionDict = regionsDict.subDict
(
"refinementSurfaces::refinementSurfaces"
"(const searchableSurfaces&, const dictionary>&"
) << "No region called " << key << " on surface "
<< allGeometry_[surfaces_[surfI]].name() << endl
<< "Valid regions are " << regionNames
<< exit(FatalError);
}
const labelPair refLevel(regionDict.lookup("level"));
regionMinLevel[surfI].insert(regionI, refLevel[0]);
regionMaxLevel[surfI].insert(regionI, refLevel[1]);
if (regionDict.found("perpendicularAngle"))
{
regionAngle[surfI].insert
(
regionI,
readScalar(regionDict.lookup("perpendicularAngle"))
regionNames[regionI]
);
const labelPair refLevel(regionDict.lookup("level"));
regionMinLevel[surfI].insert(regionI, refLevel[0]);
regionMaxLevel[surfI].insert(regionI, refLevel[1]);
if (regionDict.found("perpendicularAngle"))
{
regionAngle[surfI].insert
(
regionI,
readScalar
(
regionDict.lookup("perpendicularAngle")
)
);
}
}
}
}
surfI++;
}
surfI++;
}
// Calculate local to global region offset
label nRegions = 0;
forAll(surfacesDict, surfI)
forAll(surfaces_, surfI)
{
regionOffset_[surfI] = nRegions;
nRegions += allGeometry_[surfaces_[surfI]].regions().size();

80
src/finiteVolume/cfdTools/general/MRF/MRFZone.C
View File

@ -29,6 +29,10 @@ License
#include "volFields.H"
#include "surfaceFields.H"
#include "fvMatrices.H"
#include "PackedList.H"
#include "syncTools.H"
#include "faceSet.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -39,12 +43,15 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
dict_(is),
cellZoneID_(mesh_.cellZones().findZoneID(name_)),
faceZoneID_(mesh_.faceZones().findZoneID(name_)),
outsideFaces_(0),
patchNames_(dict_.lookup("patches")),
origin_(dict_.lookup("origin")),
axis_(dict_.lookup("axis")),
omega_(dict_.lookup("omega")),
Omega_("Omega", omega_*axis_)
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
axis_ = axis_/mag(axis_);
Omega_ = omega_*axis_;
@ -65,18 +72,71 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
if (!faceZoneFound)
{
FatalErrorIn
(
"Foam::MRFZone::MRFZone(const fvMesh& , const dictionary&)"
) << "cannot find MRF faceZone " << name_
<< exit(FatalError);
// Determine faces in cell zone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// (does not construct cells)
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
// Cells in zone
PackedBoolList zoneCell(mesh_.nCells());
if (cellZoneID_ != -1)
{
const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
forAll(cellLabels, i)
{
zoneCell[cellLabels[i]] = 1u;
}
}
// Faces in zone
PackedBoolList zoneFacesSet(mesh_.nFaces());
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
if
(
zoneCell.get(own[faceI]) == 1u
|| zoneCell.get(nei[faceI]) == 1u
)
{
zoneFacesSet[faceI] = 1u;
}
}
syncTools::syncFaceList(mesh_, zoneFacesSet, orEqOp<unsigned int>());
// Transfer to labelList
label n = zoneFacesSet.count();
outsideFaces_.setSize(n);
n = 0;
forAll(zoneFacesSet, faceI)
{
if (zoneFacesSet.get(faceI) == 1u)
{
outsideFaces_[n++] = faceI;
}
}
Info<< nl
<< "MRFZone " << name_ << " : did not find a faceZone; using "
<< returnReduce(outsideFaces_.size(), sumOp<label>())
<< " faces internal to cellZone instead." << endl;
// Flag use of outsideFaces
faceZoneID_ = -2;
}
patchLabels_.setSize(patchNames_.size());
forAll(patchNames_, i)
{
patchLabels_[i] = mesh_.boundaryMesh().findPatchID(patchNames_[i]);
patchLabels_[i] = patches.findPatchID(patchNames_[i]);
if (patchLabels_[i] == -1)
{
@ -125,7 +185,13 @@ void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
const vector& origin = origin_.value();
const vector& Omega = Omega_.value();
const labelList& faces = mesh_.faceZones()[faceZoneID_];
// Use either zone faces or outsideFaces_
const labelList& faces =
(
faceZoneID_ == -2
? outsideFaces_
: mesh_.faceZones()[faceZoneID_]
);
forAll(faces, i)
{

20
src/finiteVolume/cfdTools/general/MRF/MRFZone.H
View File

@ -26,7 +26,7 @@ Class
Foam::MRFZone
Description
MRF zone definition based on both cell and face zones and parameters
MRF zone definition based on cell zone and optional face zone and parameters
obtained from a control dictionary constructed from the given stream.
The rotation of the MRF region is defined by an origin and axis of
@ -68,18 +68,26 @@ class MRFZone
const fvMesh& mesh_;
word name_;
const word name_;
dictionary dict_;
const dictionary dict_;
label cellZoneID_;
//- label of face zone with faces on outside of cell zone.
// If -2 determines outside faces itself
label faceZoneID_;
wordList patchNames_;
//- outside faces (only if faceZoneID = -2)
labelList outsideFaces_;
const wordList patchNames_;
labelList patchLabels_;
dimensionedVector origin_;
const dimensionedVector origin_;
dimensionedVector axis_;
dimensionedScalar omega_;
const dimensionedScalar omega_;
dimensionedVector Omega_;

6
src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C
View File

@ -140,11 +140,11 @@ void Foam::quadraticFitSnGradData::findFaceDirs
#ifndef SPHERICAL_GEOMETRY
if (mesh.nGeometricD() <= 2) // find the normal direcion
{
if (mesh.directions()[0] == -1)
if (mesh.geometricD()[0] == -1)
{
kdir = vector(1, 0, 0);
}
else if (mesh.directions()[1] == -1)
else if (mesh.geometricD()[1] == -1)
{
kdir = vector(0, 1, 0);
}
@ -153,7 +153,7 @@ void Foam::quadraticFitSnGradData::findFaceDirs
kdir = vector(0, 0, 1);
}
}
else // 3D so find a direction in the place of the face
else // 3D so find a direction in the plane of the face
{
const face& f = mesh.faces()[faci];
kdir = mesh.points()[f[0]] - mesh.points()[f[1]];

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

@ -713,7 +713,7 @@ Foam::fvMatrix<Type>::H() const
(
pow
(
psi_.mesh().directions(),
psi_.mesh().solutionD(),
pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
)
);

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

@ -82,7 +82,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
(
pow
(
psi_.mesh().directions(),
psi_.mesh().solutionD(),
pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
)
);

6
src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
View File

@ -83,11 +83,11 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
# ifndef SPHERICAL_GEOMETRY
if (mesh.nGeometricD() <= 2) // find the normal direction
{
if (mesh.directions()[0] == -1)
if (mesh.geometricD()[0] == -1)
{
kdir = vector(1, 0, 0);
}
else if (mesh.directions()[1] == -1)
else if (mesh.geometricD()[1] == -1)
{
kdir = vector(0, 1, 0);
}
@ -115,7 +115,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
if (magk < SMALL)
{
FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
FatalErrorIn("findFaceDirs(..)") << " calculated kdir = zero"
<< exit(FatalError);
}
else

1
src/meshTools/Make/files
View File

@ -60,6 +60,7 @@ $(searchableSurface)/searchablePlane.C
$(searchableSurface)/searchablePlate.C
$(searchableSurface)/searchableSphere.C
$(searchableSurface)/searchableSurface.C
$(searchableSurface)/searchableSurfaceCollection.C
$(searchableSurface)/searchableSurfaces.C
$(searchableSurface)/searchableSurfacesQueries.C
$(searchableSurface)/searchableSurfaceWithGaps.C

12
src/meshTools/cellDist/wallPoint/wallPoint.H
View File

@ -51,6 +51,12 @@ namespace Foam
// Forward declaration of classes
class polyPatch;
class polyMesh;
class wallPoint;
// Forward declaration of friend functions and operators
Ostream& operator<<(Ostream&, const wallPoint&);
Istream& operator>>(Istream&, wallPoint&);
/*---------------------------------------------------------------------------*\
Class wallPoint Declaration
@ -78,12 +84,15 @@ class wallPoint
const scalar tol
);
public:
// Static data members
//- initial point far away.
static point greatPoint;
// Constructors
//- Construct null
@ -102,6 +111,7 @@ public:
const wallPoint&
);
// Member Functions
// Access
@ -184,11 +194,11 @@ public:
const scalar tol
);
// Member Operators
// Needed for List IO
inline bool operator==(const wallPoint&) const;
inline bool operator!=(const wallPoint&) const;

99
src/meshTools/meshTools/meshTools.C
View File

@ -25,7 +25,7 @@ License
\*---------------------------------------------------------------------------*/
#include "meshTools.H"
#include "primitiveMesh.H"
#include "polyMesh.H"
#include "hexMatcher.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -635,6 +635,103 @@ Foam::label Foam::meshTools::walkFace
}
void Foam::meshTools::constrainToMeshCentre(const polyMesh& mesh, point& pt)
{
const Vector<label>& dirs = mesh.geometricD();
const point& min = mesh.bounds().min();
const point& max = mesh.bounds().max();
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (dirs[cmpt] == -1)
{
pt[cmpt] = 0.5*(min[cmpt]+max[cmpt]);
}
}
}
void Foam::meshTools::constrainToMeshCentre
(
const polyMesh& mesh,
pointField& pts
)
{
const Vector<label>& dirs = mesh.geometricD();
const point& min = mesh.bounds().min();
const point& max = mesh.bounds().max();
bool isConstrained = false;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (dirs[cmpt] == -1)
{
isConstrained = true;
break;
}
}
if (isConstrained)
{
forAll(pts, i)
{
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (dirs[cmpt] == -1)
{
pts[i][cmpt] = 0.5*(min[cmpt]+max[cmpt]);
}
}
}
}
}
//- Set the constrained components of directions/velocity to zero
void Foam::meshTools::constrainDirection(const polyMesh& mesh, vector& d)
{
const Vector<label>& dirs = mesh.geometricD();
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (dirs[cmpt] == -1)
{
d[cmpt] = 0.0;
}
}
}
void Foam::meshTools::constrainDirection(const polyMesh& mesh, vectorField& d)
{
const Vector<label>& dirs = mesh.geometricD();
bool isConstrained = false;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (dirs[cmpt] == -1)
{
isConstrained = true;
break;
}
}
if (isConstrained)
{
forAll(d, i)
{
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
if (dirs[cmpt] == -1)
{
d[i][cmpt] = 0.0;
}
}
}
}
}
void Foam::meshTools::getParallelEdges
(
const primitiveMesh& mesh,

356
src/meshTools/meshTools/meshTools.H
View File

@ -50,7 +50,7 @@ namespace Foam
{
class primitiveMesh;
class polyMesh;
/*---------------------------------------------------------------------------*\
Namespace meshTools Declaration
@ -81,196 +81,212 @@ namespace meshTools
static const label pXpYpZMask = 1 << pXpYpZ;
//- Check if n is in same direction as normals of all faceLabels
bool visNormal
(
const vector& n,
const vectorField& faceNormals,
const labelList& faceLabels
);
// Normal handling
//- Calculate point normals on a 'box' mesh (all edges aligned with
// coordinate axes)
vectorField calcBoxPointNormals(const primitivePatch& pp);
//- Check if n is in same direction as normals of all faceLabels
bool visNormal
(
const vector& n,
const vectorField& faceNormals,
const labelList& faceLabels
);
//- Normalized edge vector
vector normEdgeVec(const primitiveMesh&, const label edgeI);
//- Calculate point normals on a 'box' mesh (all edges aligned with
// coordinate axes)
vectorField calcBoxPointNormals(const primitivePatch& pp);
//- Write obj representation of point
void writeOBJ
(
Ostream& os,
const point& pt
);
//- Normalized edge vector
vector normEdgeVec(const primitiveMesh&, const label edgeI);
//- Write obj representation of faces subset
void writeOBJ
(
Ostream& os,
const faceList&,
const pointField&,
const labelList& faceLabels
);
// OBJ writing
//- Write obj representation of faces
void writeOBJ
(
Ostream& os,
const faceList&,
const pointField&
);
//- Write obj representation of point
void writeOBJ
(
Ostream& os,
const point& pt
);
//- Write obj representation of cell subset
void writeOBJ
(
Ostream& os,
const cellList&,
const faceList&,
const pointField&,
const labelList& cellLabels
);
//- Write obj representation of faces subset
void writeOBJ
(
Ostream& os,
const faceList&,
const pointField&,
const labelList& faceLabels
);
//- Is edge used by cell
bool edgeOnCell
(
const primitiveMesh&,
const label cellI,
const label edgeI
);
//- Write obj representation of faces
void writeOBJ
(
Ostream& os,
const faceList&,
const pointField&
);
//- Is edge used by face
bool edgeOnFace
(
const primitiveMesh&,
const label faceI,
const label edgeI
);
//- Is face used by cell
bool faceOnCell
(
const primitiveMesh&,
const label cellI,
const label faceI
);
//- Write obj representation of cell subset
void writeOBJ
(
Ostream& os,
const cellList&,
const faceList&,
const pointField&,
const labelList& cellLabels
);
//- Return edge among candidates that uses the two vertices.
label findEdge
(
const edgeList& edges,
const labelList& candidates,
const label v0,
const label v1
);
//- Return edge between two vertices. Returns -1 if no edge.
label findEdge
(
const primitiveMesh&,
const label v0,
const label v1
);
// Cell/face/edge walking
//- Return edge shared by two faces. Throws error if no edge found.
label getSharedEdge
(
const primitiveMesh&,
const label f0,
const label f1
);
//- Return face shared by two cells. Throws error if none found.
label getSharedFace
(
const primitiveMesh&,
const label cell0,
const label cell1
);
//- Is edge used by cell
bool edgeOnCell
(
const primitiveMesh&,
const label cellI,
const label edgeI
);
//- Get faces on cell using edgeI. Throws error if no two found.
void getEdgeFaces
(
const primitiveMesh&,
const label cellI,
const label edgeI,
label& face0,
label& face1
);
//- Is edge used by face
bool edgeOnFace
(
const primitiveMesh&,
const label faceI,
const label edgeI
);
//- Return label of other edge (out of candidates edgeLabels)
// connected to vertex but not edgeI. Throws error if none found.
label otherEdge
(
const primitiveMesh&,
const labelList& edgeLabels,
const label edgeI,
const label vertI
);
//- Return face on cell using edgeI but not faceI. Throws error
// if none found.
label otherFace
(
const primitiveMesh&,
const label cellI,
const label faceI,
const label edgeI
);
//- Is face used by cell
bool faceOnCell
(
const primitiveMesh&,
const label cellI,
const label faceI
);
//- Return cell on other side of face. Throws error
// if face not internal.
label otherCell
(
const primitiveMesh&,
const label cellI,
const label faceI
);
//- Returns label of edge nEdges away from startEdge (in the direction
// of startVertI)
label walkFace
(
const primitiveMesh&,
const label faceI,
const label startEdgeI,
const label startVertI,
const label nEdges
);
//- Return edge among candidates that uses the two vertices.
label findEdge
(
const edgeList& edges,
const labelList& candidates,
const label v0,
const label v1
);
//- Return edge between two vertices. Returns -1 if no edge.
label findEdge
(
const primitiveMesh&,
const label v0,
const label v1
);
//- Return edge shared by two faces. Throws error if no edge found.
label getSharedEdge
(
const primitiveMesh&,
const label f0,
const label f1
);
//- Return face shared by two cells. Throws error if none found.
label getSharedFace
(
const primitiveMesh&,
const label cell0,
const label cell1
);
//- Get faces on cell using edgeI. Throws error if no two found.
void getEdgeFaces
(
const primitiveMesh&,
const label cellI,
const label edgeI,
label& face0,
label& face1
);
//- Return label of other edge (out of candidates edgeLabels)
// connected to vertex but not edgeI. Throws error if none found.
label otherEdge
(
const primitiveMesh&,
const labelList& edgeLabels,
const label edgeI,
const label vertI
);
//- Return face on cell using edgeI but not faceI. Throws error
// if none found.
label otherFace
(
const primitiveMesh&,
const label cellI,
const label faceI,
const label edgeI
);
//- Return cell on other side of face. Throws error
// if face not internal.
label otherCell
(
const primitiveMesh&,
const label cellI,
const label faceI
);
//- Returns label of edge nEdges away from startEdge (in the direction
// of startVertI)
label walkFace
(
const primitiveMesh&,
const label faceI,
const label startEdgeI,
const label startVertI,
const label nEdges
);
// Constraints on position
//- Set the constrained components of position to mesh centre
void constrainToMeshCentre(const polyMesh&, point&);
void constrainToMeshCentre(const polyMesh&, pointField&);
//- Set the constrained components of directions/velocity to zero
void constrainDirection(const polyMesh&, vector&);
void constrainDirection(const polyMesh&, vectorField&);
//
// Hex only functionality.
//
//- Given edge on hex find other 'parallel', non-connected edges.
void getParallelEdges
(
const primitiveMesh&,
const label cellI,
const label e0,
label&,
label&,
label&
);
//- Given edge on hex find other 'parallel', non-connected edges.
void getParallelEdges
(
const primitiveMesh&,
const label cellI,
const label e0,
label&,
label&,
label&
);
//- Given edge on hex find all 'parallel' (i.e. non-connected)
// edges and average direction of them
vector edgeToCutDir
(
const primitiveMesh&,
const label cellI,
const label edgeI
);
//- Given edge on hex find all 'parallel' (i.e. non-connected)
// edges and average direction of them
vector edgeToCutDir
(
const primitiveMesh&,
const label cellI,
const label edgeI
);
//- Reverse of edgeToCutDir: given direction find edge bundle and
// return one of them.
label cutDirToEdge
(
const primitiveMesh&,
const label cellI,
const vector& cutDir
);
//- Reverse of edgeToCutDir: given direction find edge bundle and
// return one of them.
label cutDirToEdge
(
const primitiveMesh&,
const label cellI,
const vector& cutDir
);
} // End namespace meshTools

125
src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
View File

@ -76,10 +76,8 @@ bool Foam::distributedTriSurfaceMesh::read()
// Distribution type
distType_ = distributionTypeNames_.read(dict_.lookup("distributionType"));
if (dict_.found("mergeDistance"))
{
dict_.lookup("mergeDistance") >> mergeDist_;
}
// Merge distance
mergeDist_ = readScalar(dict_.lookup("mergeDistance"));
return true;
}
@ -1341,25 +1339,44 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
dict_(io, dict)
{
read();
if (debug)
{
Info<< "Constructed from triSurface:" << endl;
writeStats(Info);
labelList nTris(Pstream::nProcs());
nTris[Pstream::myProcNo()] = triSurface::size();
Pstream::gatherList(nTris);
Pstream::scatterList(nTris);
Info<< endl<< "\tproc\ttris\tbb" << endl;
forAll(nTris, procI)
{
Info<< '\t' << procI << '\t' << nTris[procI]
<< '\t' << procBb_[procI] << endl;
}
Info<< endl;
}
}
Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
:
triSurfaceMesh(io),
// triSurfaceMesh
// (
// IOobject
// (
// io.name(),
// io.db().time().findInstanceDir(io.local()),
// io.local(),
// io.db(),
// io.readOpt(),
// io.writeOpt(),
// io.registerObject()
// )
// ),
//triSurfaceMesh(io),
triSurfaceMesh
(
IOobject
(
io.name(),
io.time().findInstance(io.local(), word::null),
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
io.registerObject()
)
),
dict_
(
IOobject
@ -1375,6 +1392,26 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
)
{
read();
if (debug)
{
Info<< "Read distributedTriSurface from " << io.objectPath()
<< ':' << endl;
writeStats(Info);
labelList nTris(Pstream::nProcs());
nTris[Pstream::myProcNo()] = triSurface::size();
Pstream::gatherList(nTris);
Pstream::scatterList(nTris);
Info<< endl<< "\tproc\ttris\tbb" << endl;
forAll(nTris, procI)
{
Info<< '\t' << procI << '\t' << nTris[procI]
<< '\t' << procBb_[procI] << endl;
}
Info<< endl;
}
}
@ -1384,21 +1421,21 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
const dictionary& dict
)
:
triSurfaceMesh(io, dict),
// triSurfaceMesh
// (
// IOobject
// (
// io.name(),
// io.db().time().findInstanceDir(io.local()),
// io.local(),
// io.db(),
// io.readOpt(),
// io.writeOpt(),
// io.registerObject()
// ),
// dict
// ),
//triSurfaceMesh(io, dict),
triSurfaceMesh
(
IOobject
(
io.name(),
io.time().findInstance(io.local(), word::null),
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
io.registerObject()
),
dict
),
dict_
(
IOobject
@ -1414,6 +1451,26 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
)
{
read();
if (debug)
{
Info<< "Read distributedTriSurface from " << io.objectPath()
<< " and dictionary:" << endl;
writeStats(Info);
labelList nTris(Pstream::nProcs());
nTris[Pstream::myProcNo()] = triSurface::size();
Pstream::gatherList(nTris);
Pstream::scatterList(nTris);
Info<< endl<< "\tproc\ttris\tbb" << endl;
forAll(nTris, procI)
{
Info<< '\t' << procI << '\t' << nTris[procI]
<< '\t' << procBb_[procI] << endl;
}
Info<< endl;
}
}
@ -2416,6 +2473,8 @@ void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
boundBox bb;
label nPoints;
calcBounds(bb, nPoints);
reduce(bb.min(), minOp<point>());
reduce(bb.max(), maxOp<point>());
os << "Triangles : " << returnReduce(triSurface::size(), sumOp<label>())
<< endl

10
src/meshTools/searchableSurface/distributedTriSurfaceMesh.H
View File

@ -91,15 +91,14 @@ public:
private:
// Private member data
//- Merging distance
scalar mergeDist_;
//- Decomposition used when independently decomposing surface.
autoPtr<decompositionMethod> decomposer_;
// Private member data
//- Bounding box settings
IOdictionary dict_;
@ -317,10 +316,11 @@ public:
const dictionary& dict
);
//- Construct read
//- Construct read. Does findInstance to find io.local().
distributedTriSurfaceMesh(const IOobject& io);
//- Construct from dictionary (used by searchableSurface)
//- Construct from dictionary (used by searchableSurface).
// Does read. Does findInstance to find io.local().
distributedTriSurfaceMesh
(
const IOobject& io,

9
src/meshTools/searchableSurface/searchableSphere.C
View File

@ -53,7 +53,14 @@ Foam::pointIndexHit Foam::searchableSphere::findNearest
if (nearestDistSqr > sqr(magN-radius_))
{
info.rawPoint() = centre_ + n/magN*radius_;
if (magN < ROOTVSMALL)
{
info.rawPoint() = centre_ + vector(1,0,0)/magN*radius_;
}
else
{
info.rawPoint() = centre_ + n/magN*radius_;
}
info.setHit();
info.setIndex(0);
}

525
src/meshTools/searchableSurface/searchableSurfaceCollection.C Normal file
View File

@ -0,0 +1,525 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "searchableSurfaceCollection.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
#include "Time.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableSurfaceCollection, 0);
addToRunTimeSelectionTable(searchableSurface, searchableSurfaceCollection, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::searchableSurfaceCollection::findNearest
(
const pointField& samples,
scalarField& minDistSqr,
List<pointIndexHit>& nearestInfo,
labelList& nearestSurf
) const
{
// Initialise
nearestInfo.setSize(samples.size());
nearestInfo = pointIndexHit();
nearestSurf.setSize(samples.size());
nearestSurf = -1;
List<pointIndexHit> hitInfo(samples.size());
const scalarField localMinDistSqr(samples.size(), GREAT);
forAll(subGeom_, surfI)
{
// Transform then divide
tmp<pointField> localSamples = cmptDivide
(
transform_[surfI].localPosition(samples),
scale_[surfI]
);
subGeom_[surfI].findNearest(localSamples, localMinDistSqr, hitInfo);
forAll(hitInfo, pointI)
{
if (hitInfo[pointI].hit())
{
// Rework back into global coordinate sys. Multiply then
// transform
point globalPt = transform_[surfI].globalPosition
(
cmptMultiply
(
hitInfo[pointI].rawPoint(),
scale_[surfI]
)
);
scalar distSqr = magSqr(globalPt - samples[pointI]);
if (distSqr < minDistSqr[pointI])
{
minDistSqr[pointI] = distSqr;
nearestInfo[pointI].setPoint(globalPt);
nearestInfo[pointI].setHit();
nearestInfo[pointI].setIndex(hitInfo[pointI].index());
nearestSurf[pointI] = surfI;
}
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableSurfaceCollection::searchableSurfaceCollection
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
instance_(dict.size()),
scale_(dict.size()),
transform_(dict.size()),
subGeom_(dict.size())
{
Info<< "SearchableCollection : " << name() << endl;
label surfI = 0;
forAllConstIter(dictionary, dict, iter)
{
if (dict.isDict(iter().keyword()))
{
instance_[surfI] = iter().keyword();
const dictionary& subDict = dict.subDict(instance_[surfI]);
scale_[surfI] = subDict.lookup("scale");
transform_.set
(
surfI,
coordinateSystem::New
(
"",
subDict.subDict("transform")
)
);
const word subGeomName(subDict.lookup("surface"));
//Pout<< "Trying to find " << subGeomName << endl;
const searchableSurface& s =
io.db().lookupObject<searchableSurface>(subGeomName);
subGeom_.set(surfI, &const_cast<searchableSurface&>(s));
Info<< " instance : " << instance_[surfI] << endl;
Info<< " surface : " << s.name() << endl;
Info<< " scale : " << scale_[surfI] << endl;
Info<< " coordsys : " << transform_[surfI] << endl;
surfI++;
}
}
instance_.setSize(surfI);
scale_.setSize(surfI);
transform_.setSize(surfI);
subGeom_.setSize(surfI);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchableSurfaceCollection::~searchableSurfaceCollection()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::wordList& Foam::searchableSurfaceCollection::regions() const
{
if (regions_.size() == 0)
{
regionOffset_.setSize(subGeom_.size());
DynamicList<word> allRegions;
forAll(subGeom_, surfI)
{
regionOffset_[surfI] = allRegions.size();
const wordList& subRegions = subGeom_[surfI].regions();
forAll(subRegions, i)
{
//allRegions.append(subRegions[i] + "_" + Foam::name(surfI));
allRegions.append(instance_[surfI] + "_" + subRegions[i]);
}
}
regions_.transfer(allRegions.shrink());
}
return regions_;
}
Foam::label Foam::searchableSurfaceCollection::size() const
{
label n = 0;
forAll(subGeom_, surfI)
{
n += subGeom_[surfI].size();
}
return n;
}
void Foam::searchableSurfaceCollection::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& nearestInfo
) const
{
// How to scale distance?
scalarField minDistSqr(nearestDistSqr);
labelList nearestSurf;
findNearest
(
samples,
minDistSqr,
nearestInfo,
nearestSurf
);
}
void Foam::searchableSurfaceCollection::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
info = pointIndexHit();
// Current nearest (to start) intersection
pointField nearest(end);
List<pointIndexHit> hitInfo(start.size());
forAll(subGeom_, surfI)
{
// Starting point
tmp<pointField> e0 = cmptDivide
(
transform_[surfI].localPosition
(
start
),
scale_[surfI]
);
// Current best end point
tmp<pointField> e1 = cmptDivide
(
transform_[surfI].localPosition
(
nearest
),
scale_[surfI]
);
subGeom_[surfI].findLine(e0, e1, hitInfo);
forAll(hitInfo, pointI)
{
if (hitInfo[pointI].hit())
{
// Transform back to global coordinate sys.
nearest[pointI] = transform_[surfI].globalPosition
(
cmptMultiply
(
hitInfo[pointI].rawPoint(),
scale_[surfI]
)
);
info[pointI] = hitInfo[pointI];
info[pointI].rawPoint() = nearest[pointI];
}
}
}
// Debug check
if (false)
{
forAll(info, pointI)
{
if (info[pointI].hit())
{
vector n(end[pointI] - start[pointI]);
scalar magN = mag(n);
if (magN > SMALL)
{
n /= mag(n);
scalar s = ((info[pointI].rawPoint()-start[pointI])&n);
if (s < 0 || s > 1)
{
FatalErrorIn
(
"searchableSurfaceCollection::findLine(..)"
) << "point:" << info[pointI]
<< " s:" << s
<< " outside vector "
<< " start:" << start[pointI]
<< " end:" << end[pointI]
<< abort(FatalError);
}
}
}
}
}
}
void Foam::searchableSurfaceCollection::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
// To be done ...
findLine(start, end, info);
}
void Foam::searchableSurfaceCollection::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
// To be done. Assume for now only one intersection.
List<pointIndexHit> nearestInfo;
findLine(start, end, nearestInfo);
info.setSize(start.size());
forAll(info, pointI)
{
if (nearestInfo[pointI].hit())
{
info[pointI].setSize(1);
info[pointI][0] = nearestInfo[pointI];
}
else
{
info[pointI].clear();
}
}
}
void Foam::searchableSurfaceCollection::getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
if (subGeom_.size() == 0)
{}
else if (subGeom_.size() == 1)
{
subGeom_[0].getRegion(info, region);
}
else
{
region.setSize(info.size());
region = -1;
// Which region did point come from. Retest for now to see which
// surface it originates from - crap solution! Should use global indices
// in index inside pointIndexHit to do this better.
pointField samples(info.size());
forAll(info, pointI)
{
if (info[pointI].hit())
{
samples[pointI] = info[pointI].hitPoint();
}
else
{
samples[pointI] = vector::zero;
}
}
//scalarField minDistSqr(info.size(), SMALL);
scalarField minDistSqr(info.size(), GREAT);
labelList nearestSurf;
List<pointIndexHit> nearestInfo;
findNearest
(
samples,
minDistSqr,
nearestInfo,
nearestSurf
);
// Check
{
forAll(info, pointI)
{
if (info[pointI].hit() && nearestSurf[pointI] == -1)
{
FatalErrorIn
(
"searchableSurfaceCollection::getRegion(..)"
) << "pointI:" << pointI
<< " sample:" << samples[pointI]
<< " nearest:" << nearestInfo[pointI]
<< " nearestsurf:" << nearestSurf[pointI]
<< abort(FatalError);
}
}
}
forAll(subGeom_, surfI)
{
// Collect points from my surface
labelList indices(findIndices(nearestSurf, surfI));
labelList surfRegion;
subGeom_[surfI].getRegion
(
IndirectList<pointIndexHit>(info, indices),
surfRegion
);
forAll(indices, i)
{
region[indices[i]] = regionOffset_[surfI] + surfRegion[i];
}
}
}
}
void Foam::searchableSurfaceCollection::getNormal
(
const List<pointIndexHit>& info,
vectorField& normal
) const
{
if (subGeom_.size() == 0)
{}
else if (subGeom_.size() == 1)
{
subGeom_[0].getNormal(info, normal);
}
else
{
normal.setSize(info.size());
// See above - crap retest to find surface point originates from.
pointField samples(info.size());
forAll(info, pointI)
{
if (info[pointI].hit())
{
samples[pointI] = info[pointI].hitPoint();
}
else
{
samples[pointI] = vector::zero;
}
}
//scalarField minDistSqr(info.size(), SMALL);
scalarField minDistSqr(info.size(), GREAT);
labelList nearestSurf;
List<pointIndexHit> nearestInfo;
findNearest
(
samples,
minDistSqr,
nearestInfo,
nearestSurf
);
forAll(subGeom_, surfI)
{
// Collect points from my surface
labelList indices(findIndices(nearestSurf, surfI));
vectorField surfNormal;
subGeom_[surfI].getNormal
(
IndirectList<pointIndexHit>(info, indices),
surfNormal
);
forAll(indices, i)
{
normal[indices[i]] = surfNormal[i];
}
}
}
}
void Foam::searchableSurfaceCollection::getVolumeType
(
const pointField& points,
List<volumeType>& volType
) const
{
FatalErrorIn
(
"searchableSurfaceCollection::getVolumeType(const pointField&"
", List<volumeType>&) const"
) << "Volume type not supported for collection."
<< exit(FatalError);
}
// ************************************************************************* //

212
src/meshTools/searchableSurface/searchableSurfaceCollection.H Normal file
View File

@ -0,0 +1,212 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::searchableSurfaceCollection
Description
Union of transformed searchableSurfaces
SourceFiles
searchableSurfaceCollection.C
\*---------------------------------------------------------------------------*/
#ifndef searchableSurfaceCollection_H
#define searchableSurfaceCollection_H
#include "searchableSurface.H"
#include "treeBoundBox.H"
#include "coordinateSystem.H"
#include "UPtrList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchableSurfaceCollection Declaration
\*---------------------------------------------------------------------------*/
class searchableSurfaceCollection
:
public searchableSurface
{
private:
// Private Member Data
// Per instance data
//- instance name
wordList instance_;
//- scaling vector
vectorField scale_;
//- transformation
PtrList<coordinateSystem> transform_;
UPtrList<searchableSurface> subGeom_;
//- Region names
mutable wordList regions_;
//- From individual regions to collection regions
mutable labelList regionOffset_;
// Private Member Functions
//- Find point nearest to sample. Updates minDistSqr. Sets nearestInfo
// and surface index
void findNearest
(
const pointField& samples,
scalarField& minDistSqr,
List<pointIndexHit>& nearestInfo,
labelList& nearestSurf
) const;
//- Disallow default bitwise copy construct
searchableSurfaceCollection(const searchableSurfaceCollection&);
//- Disallow default bitwise assignment
void operator=(const searchableSurfaceCollection&);
public:
//- Runtime type information
TypeName("searchableSurfaceCollection");
// Constructors
//- Construct from dictionary (used by searchableSurface)
searchableSurfaceCollection
(
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchableSurfaceCollection();
// Member Functions
virtual const wordList& regions() const;
//- Whether supports volume type below
virtual bool hasVolumeType() const
{
return false;
}
//- Range of local indices that can be returned.
virtual label size() const;
// Multiple point queries.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>&,
labelList& region
) const;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const;
// regIOobject implementation
bool writeData(Ostream&) const
{
notImplemented
(
"searchableSurfaceCollection::writeData(Ostream&) const"
);
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

10
src/meshTools/searchableSurface/searchableSurfaces.C
View File

@ -39,8 +39,6 @@ defineTypeNameAndDebug(searchableSurfaces, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct with length.
@ -187,14 +185,6 @@ Foam::searchableSurfaces::searchableSurfaces
// their object name. Maybe have stlTriSurfaceMesh which appends .stl
// when reading/writing?
namedIO().rename(key); // names_[surfI]
if (namedIO().local() != word::null)
{
namedIO().instance() = namedIO().time().findInstance
(
namedIO().local(),
namedIO().name()
);
}
// Create and hook surface
set

2
src/meshTools/searchableSurface/searchableSurfaces.H
View File

@ -69,6 +69,8 @@ class searchableSurfaces
labelList allSurfaces_;
// Private Member Functions
//- Disallow default bitwise copy construct
searchableSurfaces(const searchableSurfaces&);

127
src/meshTools/searchableSurface/triSurfaceMesh.C
View File

@ -43,6 +43,64 @@ addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//// Special version of Time::findInstance that does not check headerOk
//// to search for instances of raw files
//Foam::word Foam::triSurfaceMesh::findRawInstance
//(
// const Time& runTime,
// const fileName& dir,
// const word& name
//)
//{
// // Check current time first
// if (isFile(runTime.path()/runTime.timeName()/dir/name))
// {
// return runTime.timeName();
// }
// instantList ts = runTime.times();
// label instanceI;
//
// for (instanceI = ts.size()-1; instanceI >= 0; --instanceI)
// {
// if (ts[instanceI].value() <= runTime.timeOutputValue())
// {
// break;
// }
// }
//
// // continue searching from here
// for (; instanceI >= 0; --instanceI)
// {
// if (isFile(runTime.path()/ts[instanceI].name()/dir/name))
// {
// return ts[instanceI].name();
// }
// }
//
//
// // not in any of the time directories, try constant
//
// // Note. This needs to be a hard-coded constant, rather than the
// // constant function of the time, because the latter points to
// // the case constant directory in parallel cases
//
// if (isFile(runTime.path()/runTime.constant()/dir/name))
// {
// return runTime.constant();
// }
//
// FatalErrorIn
// (
// "searchableSurfaces::findRawInstance"
// "(const Time&, const fileName&, const word&)"
// ) << "Cannot find file \"" << name << "\" in directory "
// << runTime.constant()/dir
// << exit(FatalError);
//
// return runTime.constant();
//}
//- Check file existence
const Foam::fileName& Foam::triSurfaceMesh::checkFile
(
@ -149,7 +207,19 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
:
searchableSurface(io),
objectRegistry(io),
objectRegistry
(
IOobject
(
io.name(),
io.instance(),
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
false // searchableSurface already registered under name
)
),
triSurface(s),
surfaceClosed_(-1)
{}
@ -157,8 +227,34 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
:
// Find instance for triSurfaceMesh
searchableSurface(io),
objectRegistry(io),
//(
// IOobject
// (
// io.name(),
// io.time().findInstance(io.local(), word::null),
// io.local(),
// io.db(),
// io.readOpt(),
// io.writeOpt(),
// io.registerObject()
// )
//),
// Reused found instance in objectRegistry
objectRegistry
(
IOobject
(
io.name(),
static_cast<const searchableSurface&>(*this).instance(),
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
false // searchableSurface already registered under name
)
),
triSurface
(
checkFile
@ -178,7 +274,32 @@ Foam::triSurfaceMesh::triSurfaceMesh
)
:
searchableSurface(io),
objectRegistry(io),
//(
// IOobject
// (
// io.name(),
// io.time().findInstance(io.local(), word::null),
// io.local(),
// io.db(),
// io.readOpt(),
// io.writeOpt(),
// io.registerObject()
// )
//),
// Reused found instance in objectRegistry
objectRegistry
(
IOobject
(
io.name(),
static_cast<const searchableSurface&>(*this).instance(),
io.local(),
io.db(),
io.readOpt(),
io.writeOpt(),
false // searchableSurface already registered under name
)
),
triSurface
(
checkFile

12
src/meshTools/searchableSurface/triSurfaceMesh.H
View File

@ -76,6 +76,14 @@ private:
// Private Member Functions
////- Helper: find instance of files without header
//static word findRawInstance
//(
// const Time&,
// const fileName&,
// const word&
//);
//- Check file existence
static const fileName& checkFile
(
@ -105,10 +113,10 @@ public:
//- Construct from triSurface
triSurfaceMesh(const IOobject&, const triSurface&);
//- Construct read
//- Construct read.
triSurfaceMesh(const IOobject& io);
//- Construct from dictionary (used by searchableSurface)
//- Construct from IO and dictionary (used by searchableSurface).
// Dictionary may contain a 'scale' entry (eg, 0.001: mm -> m)
triSurfaceMesh
(

403
src/sampling/sampledSurface/isoSurface/isoSurface.C
View File

@ -30,6 +30,7 @@ License
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
#include "slicedVolFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -84,9 +85,74 @@ bool Foam::isoSurface::isEdgeOfFaceCut
}
// Get neighbour value and position.
void Foam::isoSurface::getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const surfaceScalarField& weights = mesh_.weights();
if (mesh_.isInternalFace(faceI))
{
label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
nbrValue = cVals[nbr];
nbrPoint = mesh_.cellCentres()[nbr];
}
else
{
label bFaceI = faceI-mesh_.nInternalFaces();
label patchI = boundaryRegion[bFaceI];
const polyPatch& pp = mesh_.boundaryMesh()[patchI];
label patchFaceI = faceI-pp.start();
if (isA<emptyPolyPatch>(pp))
{
// Assume zero gradient
nbrValue = cVals[own[faceI]];
nbrPoint = mesh_.faceCentres()[faceI];
}
else if (pp.coupled())
{
if (!refCast<const coupledPolyPatch>(pp).separated())
{
// Behave as internal face:
// other side value
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
// other side cell centre
nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
}
else
{
// Do some interpolation for now
const scalarField& w = weights.boundaryField()[patchI];
nbrPoint = mesh_.faceCentres()[faceI];
nbrValue =
(1-w[patchFaceI])*cVals[own[faceI]]
+ w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
}
}
else
{
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
nbrPoint = mesh_.faceCentres()[faceI];
}
}
}
// Determine for every face/cell whether it (possibly) generates triangles.
void Foam::isoSurface::calcCutTypes
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const scalarField& pVals
)
@ -102,7 +168,20 @@ void Foam::isoSurface::calcCutTypes
{
// CC edge.
bool ownLower = (cVals[own[faceI]] < iso_);
bool neiLower = (cVals[nei[faceI]] < iso_);
scalar nbrValue;
point nbrPoint;
getNeighbour
(
boundaryRegion,
cVals,
own[faceI],
faceI,
nbrValue,
nbrPoint
);
bool neiLower = (nbrValue < iso_);
if (ownLower != neiLower)
{
@ -120,53 +199,47 @@ void Foam::isoSurface::calcCutTypes
}
}
}
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
label faceI = pp.start();
if (isA<emptyPolyPatch>(pp))
forAll(pp, i)
{
// Assume zero gradient so owner and neighbour/boundary value equal
bool ownLower = (cVals[own[faceI]] < iso_);
forAll(pp, i)
scalar nbrValue;
point nbrPoint;
getNeighbour
(
boundaryRegion,
cVals,
own[faceI],
faceI,
nbrValue,
nbrPoint
);
bool neiLower = (nbrValue < iso_);
if (ownLower != neiLower)
{
bool ownLower = (cVals[own[faceI]] < iso_);
faceCutType_[faceI] = CUT;
}
else
{
// Mesh edge.
const face f = mesh_.faces()[faceI];
if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower))
if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
{
faceCutType_[faceI] = CUT;
}
faceI++;
}
}
else
{
forAll(pp, i)
{
bool ownLower = (cVals[own[faceI]] < iso_);
bool neiLower = (cVals.boundaryField()[patchI][i] < iso_);
if (ownLower != neiLower)
{
faceCutType_[faceI] = CUT;
}
else
{
// Mesh edge.
const face f = mesh_.faces()[faceI];
if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
{
faceCutType_[faceI] = CUT;
}
}
faceI++;
}
faceI++;
}
}
@ -330,59 +403,6 @@ Foam::pointIndexHit Foam::isoSurface::collapseSurface
}
// Get neighbour value and position.
void Foam::isoSurface::getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
if (mesh_.isInternalFace(faceI))
{
label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
nbrValue = cVals[nbr];
nbrPoint = mesh_.cellCentres()[nbr];
}
else
{
label bFaceI = faceI-mesh_.nInternalFaces();
label patchI = boundaryRegion[bFaceI];
const polyPatch& pp = mesh_.boundaryMesh()[patchI];
label patchFaceI = faceI-pp.start();
if (isA<emptyPolyPatch>(pp))
{
// Assume zero gradient
nbrValue = cVals[own[faceI]];
nbrPoint = mesh_.faceCentres()[faceI];
}
else if
(
pp.coupled()
&& !refCast<const coupledPolyPatch>(pp).separated()
)
{
// other side value
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
// other side cell centre
nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
}
else
{
nbrValue = cVals.boundaryField()[patchI][patchFaceI];
nbrPoint = mesh_.faceCentres()[faceI];
}
}
}
// Determine per cell centre whether all the intersections get collapsed
// to a single point
void Foam::isoSurface::calcSnappedCc
@ -606,13 +626,14 @@ void Foam::isoSurface::calcSnappedPoint
forAll(pFaces, pFaceI)
{
// Create points for all intersections close to point
// (i.e. from pyramid edges)
label faceI = pFaces[pFaceI];
const face& f = mesh_.faces()[faceI];
label own = mesh_.faceOwner()[faceI];
// Create points for all intersections close to point
// (i.e. from pyramid edges)
// Get neighbour value
scalar nbrValue;
point nbrPoint;
getNeighbour
@ -831,6 +852,7 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
}
// Determine 'flat hole' situation (see RMT paper).
// Two unconnected triangles get connected because (some of) the edges
// separating them get collapsed. Below only checks for duplicate triangles,
@ -846,22 +868,29 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
forAll(tris, triI)
{
const labelledTri& tri = tris[triI];
const labelList& pFaces = pointFaces[tri[0]];
// Find the minimum of any duplicates
// Find the maximum of any duplicates. Maximum since the tris
// below triI
// get overwritten so we cannot use them in a comparison.
label dupTriI = -1;
forAll(pFaces, i)
{
if (pFaces[i] < triI && tris[pFaces[i]] == tri)
label nbrTriI = pFaces[i];
if (nbrTriI > triI && (tris[nbrTriI] == tri))
{
dupTriI = pFaces[i];
//Pout<< "Duplicate : " << triI << " verts:" << tri
// << " to " << nbrTriI << " verts:" << tris[nbrTriI]
// << endl;
dupTriI = nbrTriI;
break;
}
}
if (dupTriI == -1)
{
// There is no lower triangle
// There is no (higher numbered) duplicate triangle
label newTriI = newToOldTri.size();
newToOldTri.append(triI);
tris[newTriI] = tris[triI];
@ -876,6 +905,43 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
Pout<< "isoSurface : merged from " << nTris
<< " down to " << tris.size() << " unique triangles." << endl;
}
if (debug)
{
triSurface surf(tris, geometricSurfacePatchList(0), newPoints);
forAll(surf, faceI)
{
const labelledTri& f = surf[faceI];
const labelList& fFaces = surf.faceFaces()[faceI];
forAll(fFaces, i)
{
label nbrFaceI = fFaces[i];
if (nbrFaceI <= faceI)
{
// lower numbered faces already checked
continue;
}
const labelledTri& nbrF = surf[nbrFaceI];
if (f == nbrF)
{
FatalErrorIn("validTri(const triSurface&, const label)")
<< "Check : "
<< " triangle " << faceI << " vertices " << f
<< " fc:" << f.centre(surf.points())
<< " has the same vertices as triangle " << nbrFaceI
<< " vertices " << nbrF
<< " fc:" << nbrF.centre(surf.points())
<< abort(FatalError);
}
}
}
}
}
return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
@ -1473,45 +1539,23 @@ Foam::isoSurface::isoSurface
{
if (debug)
{
Pout<< "isoSurface :" << nl
<< " isoField : " << cVals.name() << nl
<< " isoValue : " << iso << nl
<< " regularise : " << regularise << nl
<< " mergeTol : " << mergeTol << nl
Pout<< "isoSurface:" << nl
<< " isoField : " << cVals.name() << nl
<< " cell min/max : "
<< min(cVals.internalField()) << " / "
<< max(cVals.internalField()) << nl
<< " point min/max : "
<< min(pVals) << " / "
<< max(pVals) << nl
<< " isoValue : " << iso << nl
<< " regularise : " << regularise << nl
<< " mergeTol : " << mergeTol << nl
<< endl;
}
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Check
forAll(patches, patchI)
{
if (isA<emptyPolyPatch>(patches[patchI]))
{
FatalErrorIn
(
"isoSurface::isoSurface\n"
"(\n"
" const volScalarField& cVals,\n"
" const scalarField& pVals,\n"
" const scalar iso,\n"
" const bool regularise,\n"
" const scalar mergeTol\n"
")\n"
) << "Iso surfaces not supported on case with empty patches."
<< exit(FatalError);
}
}
// Determine if any cut through face/cell
calcCutTypes(cVals, pVals);
// Determine if point is on boundary. Points on boundaries are never
// snapped. Coupled boundaries are handled explicitly so not marked here.
PackedBoolList isBoundaryPoint(mesh_.nPoints());
// Pre-calculate patch-per-face to avoid whichPatch call.
labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
forAll(patches, patchI)
@ -1525,33 +1569,13 @@ Foam::isoSurface::isoSurface
boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
faceI++;
}
// Mark all points that are not physically coupled (so anything
// but collocated coupled patches)
if
(
!pp.coupled()
|| refCast<const coupledPolyPatch>(pp).separated()
)
{
label faceI = pp.start();
forAll(pp, i)
{
boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
const face& f = mesh_.faces()[faceI];
forAll(f, fp)
{
isBoundaryPoint.set(f[fp], 1);
}
faceI++;
}
}
}
// Determine if any cut through face/cell
calcCutTypes(boundaryRegion, cVals, pVals);
DynamicList<point> snappedPoints(nCutCells_);
// Per cc -1 or a point inside snappedPoints.
@ -1584,6 +1608,39 @@ Foam::isoSurface::isoSurface
label nCellSnaps = snappedPoints.size();
// Determine if point is on boundary. Points on boundaries are never
// snapped. Coupled boundaries are handled explicitly so not marked here.
PackedBoolList isBoundaryPoint(mesh_.nPoints());
forAll(patches, patchI)
{
// Mark all boundary points that are not physically coupled (so anything
// but collocated coupled patches)
const polyPatch& pp = patches[patchI];
if
(
!pp.coupled()
|| refCast<const coupledPolyPatch>(pp).separated()
)
{
label faceI = pp.start();
forAll(pp, i)
{
const face& f = mesh_.faces()[faceI];
forAll(f, fp)
{
isBoundaryPoint.set(f[fp], 1);
}
faceI++;
}
}
}
// Per point -1 or a point inside snappedPoints.
labelList snappedPoint;
if (regularise)
@ -1613,7 +1670,7 @@ Foam::isoSurface::isoSurface
// Generate field to interpolate. This is identical to the mesh.C()
// except on separated coupled patches.
// except on separated coupled patches and on empty patches.
slicedVolVectorField meshC
(
IOobject
@ -1635,10 +1692,12 @@ Foam::isoSurface::isoSurface
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
patches[patchI].coupled()
&& refCast<const coupledPolyPatch>(patches[patchI]).separated()
pp.coupled()
&& refCast<const coupledPolyPatch>(pp).separated()
)
{
fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
@ -1647,9 +1706,32 @@ Foam::isoSurface::isoSurface
);
pfld.operator==
(
patches[patchI].patchSlice(mesh_.faceCentres())
pp.patchSlice(mesh_.faceCentres())
);
}
else if (isA<emptyPolyPatch>(pp))
{
typedef slicedVolVectorField::GeometricBoundaryField bType;
bType& bfld = const_cast<bType&>(meshC.boundaryField());
// Clear old value. Cannot resize it since slice.
bfld.set(patchI, NULL);
// Set new value we can change
bfld.set
(
patchI,
new calculatedFvPatchField<vector>
(
mesh_.boundary()[patchI],
meshC
)
);
// Change to face centres
bfld[patchI] = pp.patchSlice(mesh_.faceCentres());
}
}
}
@ -1676,7 +1758,8 @@ Foam::isoSurface::isoSurface
if (debug)
{
Pout<< "isoSurface : generated " << triMeshCells.size()
<< " unmerged triangles." << endl;
<< " unmerged triangles from " << triPoints.size()
<< " unmerged points." << endl;
}
@ -1773,6 +1856,14 @@ Foam::isoSurface::isoSurface
orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
//}
if (debug)
{
fileName stlFile = mesh_.time().path() + ".stl";
Pout<< "Dumping surface to " << stlFile << endl;
triSurface::write(stlFile);
}
}

50
src/sampling/sampledSurface/isoSurface/isoSurface.H
View File

@ -31,9 +31,6 @@ Description
G.M. Treece, R.W. Prager and A.H. Gee.
Note:
- not possible on patches of type 'empty'. There are no values on
'empty' patch fields so even the api would have to change
(no volScalarField as argument). Too messy.
- in parallel the regularisation (coarsening) always takes place
and slightly different surfaces will be created compared to non-parallel.
The surface will still be continuous though!
@ -43,6 +40,7 @@ Description
- uses geometric merge with fraction of bounding box as distance.
- triangles can be between two cell centres so constant sampling
does not make sense.
- on empty patches behaves like zero gradient.
- does not do 2D correctly, creates non-flat iso surface.
SourceFiles
@ -135,9 +133,20 @@ class isoSurface
const bool neiLower
) const;
void getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const;
//- Set faceCutType,cellCutType.
void calcCutTypes
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const scalarField& pVals
);
@ -156,16 +165,6 @@ class isoSurface
DynamicList<labelledTri, 64>& localTris
);
void getNeighbour
(
const labelList& boundaryRegion,
const volScalarField& cVals,
const label cellI,
const label faceI,
scalar& nbrValue,
point& nbrPoint
) const;
//- Determine per cc whether all near cuts can be snapped to single
// point.
void calcSnappedCc
@ -193,37 +192,39 @@ class isoSurface
template<class Type>
Type generatePoint
(
const DynamicList<Type>& snappedPoints,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index
const bool hasSnap1,
const Type& snapP1
) const;
template<class Type>
void generateTriPoints
(
const DynamicList<Type>& snapped,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index,
const bool hasSnap1,
const Type& snapP1,
const scalar s2,
const Type& p2,
const label p2Index,
const bool hasSnap2,
const Type& snapP2,
const scalar s3,
const Type& p3,
const label p3Index,
const bool hasSnap3,
const Type& snapP3,
DynamicList<Type>& points
) const;
@ -244,7 +245,8 @@ class isoSurface
const scalar neiVal,
const Type& neiPt,
const label neiSnap,
const bool hasNeiSnap,
const Type& neiSnapPt,
DynamicList<Type>& triPoints,
DynamicList<label>& triMeshCells

336
src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
View File

@ -33,15 +33,15 @@ License
template<class Type>
Type Foam::isoSurface::generatePoint
(
const DynamicList<Type>& snappedPoints,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index
const bool hasSnap1,
const Type& snapP1
) const
{
scalar d = s1-s0;
@ -50,13 +50,13 @@ Type Foam::isoSurface::generatePoint
{
scalar s = (iso_-s0)/d;
if (s >= 0.5 && s <= 1 && p1Index != -1)
if (hasSnap1 && s >= 0.5 && s <= 1)
{
return snappedPoints[p1Index];
return snapP1;
}
else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
else if (hasSnap0 && s >= 0.0 && s <= 0.5)
{
return snappedPoints[p0Index];
return snapP0;
}
else
{
@ -75,23 +75,25 @@ Type Foam::isoSurface::generatePoint
template<class Type>
void Foam::isoSurface::generateTriPoints
(
const DynamicList<Type>& snapped,
const scalar s0,
const Type& p0,
const label p0Index,
const bool hasSnap0,
const Type& snapP0,
const scalar s1,
const Type& p1,
const label p1Index,
const bool hasSnap1,
const Type& snapP1,
const scalar s2,
const Type& p2,
const label p2Index,
const bool hasSnap2,
const Type& snapP2,
const scalar s3,
const Type& p3,
const label p3Index,
const bool hasSnap3,
const Type& snapP3,
DynamicList<Type>& points
) const
@ -123,29 +125,55 @@ void Foam::isoSurface::generateTriPoints
case 0x0E:
case 0x01:
points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
);
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
);
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
);
break;
case 0x0D:
case 0x02:
points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
);
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
);
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
);
break;
case 0x0C:
case 0x03:
{
Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
Type tp1 =
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
Type tp2 =
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
);
points.append(tp1);
points.append(tp2);
points.append(tp2);
points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
);
points.append(tp1);
}
break;
@ -153,23 +181,40 @@ void Foam::isoSurface::generateTriPoints
case 0x0B:
case 0x04:
{
points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
points.append
(
generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
);
points.append
(
generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
);
points.append
(
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
);
}
break;
case 0x0A:
case 0x05:
{
Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
Type tp0 =
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
Type tp1 =
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
points.append(tp0);
points.append(tp1);
points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
);
points.append(tp0);
points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
);
points.append(tp1);
}
break;
@ -177,23 +222,40 @@ void Foam::isoSurface::generateTriPoints
case 0x09:
case 0x06:
{
Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
Type tp0 =
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
Type tp1 =
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
points.append(tp0);
points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
);
points.append(tp1);
points.append(tp0);
points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
);
points.append(tp1);
}
break;
case 0x07:
case 0x08:
points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
points.append
(
generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
);
points.append
(
generatePoint(s3,p3,hasSnap3,snapP3,s2,p2,hasSnap2,snapP2)
);
points.append
(
generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
);
break;
}
}
@ -215,7 +277,8 @@ Foam::label Foam::isoSurface::generateFaceTriPoints
const scalar neiVal,
const Type& neiPt,
const label neiSnap,
const bool hasNeiSnap,
const Type& neiSnapPt,
DynamicList<Type>& triPoints,
DynamicList<label>& triMeshCells
@ -234,23 +297,37 @@ Foam::label Foam::isoSurface::generateFaceTriPoints
generateTriPoints
(
snappedPoints,
pVals[pointI],
pCoords[pointI],
snappedPoint[pointI],
snappedPoint[pointI] != -1,
(
snappedPoint[pointI] != -1
? snappedPoints[snappedPoint[pointI]]
: pTraits<Type>::zero
),
pVals[nextPointI],
pCoords[nextPointI],
snappedPoint[nextPointI],
snappedPoint[nextPointI] != -1,
(
snappedPoint[nextPointI] != -1
? snappedPoints[snappedPoint[nextPointI]]
: pTraits<Type>::zero
),
cVals[own],
cCoords[own],
snappedCc[own],
snappedCc[own] != -1,
(
snappedCc[own] != -1
? snappedPoints[snappedCc[own]]
: pTraits<Type>::zero
),
neiVal,
neiPt,
neiSnap,
hasNeiSnap,
neiSnapPt,
triPoints
);
@ -311,25 +388,6 @@ void Foam::isoSurface::generateTriPoints
<< abort(FatalError);
}
// Determine neighbouring snap status
labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
forAll(pp, i)
{
neiSnappedCc[faceI-mesh_.nInternalFaces()] =
snappedCc[own[faceI]];
faceI++;
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false);
// Generate triangle points
@ -356,7 +414,12 @@ void Foam::isoSurface::generateTriPoints
cVals[nei[faceI]],
cCoords[nei[faceI]],
snappedCc[nei[faceI]],
snappedCc[nei[faceI]] != -1,
(
snappedCc[nei[faceI]] != -1
? snappedPoints[snappedCc[nei[faceI]]]
: pTraits<Type>::zero
),
triPoints,
triMeshCells
@ -365,6 +428,34 @@ void Foam::isoSurface::generateTriPoints
}
// Determine neighbouring snap status
boolList neiSnapped(mesh_.nFaces()-mesh_.nInternalFaces(), false);
List<Type> neiSnappedPoint(neiSnapped.size(), pTraits<Type>::zero);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
forAll(pp, i)
{
label bFaceI = faceI-mesh_.nInternalFaces();
label snappedIndex = snappedCc[own[faceI]];
if (snappedIndex != -1)
{
neiSnapped[bFaceI] = true;
neiSnappedPoint[bFaceI] = snappedPoints[snappedIndex];
}
faceI++;
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiSnapped, false);
syncTools::swapBoundaryFaceList(mesh_, neiSnappedPoint, false);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
@ -372,45 +463,104 @@ void Foam::isoSurface::generateTriPoints
if
(
isA<processorPolyPatch>(pp)
&& refCast<const processorPolyPatch>(pp).owner()
&& !refCast<const processorPolyPatch>(pp).separated()
)
{
label faceI = pp.start();
forAll(pp, i)
//if (refCast<const processorPolyPatch>(pp).owner())
{
if (faceCutType_[faceI] != NOTCUT)
label faceI = pp.start();
forAll(pp, i)
{
generateFaceTriPoints
(
cVals,
pVals,
if (faceCutType_[faceI] != NOTCUT)
{
label bFaceI = faceI-mesh_.nInternalFaces();
if
(
neiSnapped[bFaceI]
&& (neiSnappedPoint[bFaceI]==pTraits<Type>::zero)
)
{
FatalErrorIn("isoSurface::generateTriPoints(..)")
<< "problem:" << abort(FatalError);
}
cCoords,
pCoords,
snappedPoints,
snappedCc,
snappedPoint,
faceI,
generateFaceTriPoints
(
cVals,
pVals,
cVals.boundaryField()[patchI][i],
cCoords.boundaryField()[patchI][i],
neiSnappedCc[faceI-mesh_.nInternalFaces()],
cCoords,
pCoords,
triPoints,
triMeshCells
);
snappedPoints,
snappedCc,
snappedPoint,
faceI,
cVals.boundaryField()[patchI][i],
cCoords.boundaryField()[patchI][i],
neiSnapped[faceI-mesh_.nInternalFaces()],
neiSnappedPoint[faceI-mesh_.nInternalFaces()],
triPoints,
triMeshCells
);
}
faceI++;
}
faceI++;
}
}
else if (isA<emptyPolyPatch>(pp))
{
// Assume zero-gradient. But what about coordinates?
// Check if field is there (when generating geometry the
// empty patches have been rewritten to be the face centres),
// otherwise use zero-gradient.
label faceI = pp.start();
const fvPatchScalarField& fvp = cVals.boundaryField()[patchI];
// Owner value of cVals
scalarField internalVals;
if (fvp.size() == 0)
{
internalVals.setSize(pp.size());
forAll(pp, i)
{
internalVals[i] = cVals[own[pp.start()+i]];
}
}
const scalarField& bVals =
(
fvp.size() > 0
? static_cast<const scalarField&>(fvp)
: internalVals
);
const fvPatchField<Type>& pc = cCoords.boundaryField()[patchI];
// Owner value of cCoords
Field<Type> internalCoords;
if (pc.size() == 0)
{
internalCoords.setSize(pp.size());
forAll(pp, i)
{
internalCoords[i] = cCoords[own[pp.start()+i]];
}
}
const Field<Type>& bCoords =
(
pc.size() > 0
? static_cast<const Field<Type>&>(pc)
: internalCoords
);
forAll(pp, i)
{
if (faceCutType_[faceI] != NOTCUT)
@ -428,9 +578,10 @@ void Foam::isoSurface::generateTriPoints
snappedPoint,
faceI,
cVals[own[faceI]],
cCoords.boundaryField()[patchI][i],
-1, // fc not snapped
bVals[i],
bCoords[i],
false, // fc not snapped
pTraits<Type>::zero,
triPoints,
triMeshCells
@ -462,7 +613,8 @@ void Foam::isoSurface::generateTriPoints
cVals.boundaryField()[patchI][i],
cCoords.boundaryField()[patchI][i],
-1, // fc not snapped
false, // fc not snapped
pTraits<Type>::zero,
triPoints,
triMeshCells

2
src/turbulenceModels/LES/Allwmake
View File

@ -2,7 +2,7 @@
cd ${0%/*} || exit 1 # run from this directory
set -x
wmakeLnInclude -f ../incompressible/LES
wmakeLnInclude -f ../incompressible/LES -sf
wmake libso LESfilters
wmake libso LESdeltas

10
src/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C
View File

@ -42,8 +42,7 @@ addToRunTimeSelectionTable(LESdelta, cubeRootVolDelta, dictionary);
void cubeRootVolDelta::calcDelta()
{
const Vector<label>& directions = mesh().directions();
label nD = (directions.nComponents + cmptSum(directions))/2;
label nD = mesh().nGeometricD();
if (nD == 3)
{
@ -55,14 +54,15 @@ void cubeRootVolDelta::calcDelta()
<< "Case is 2D, LES is not strictly applicable\n"
<< endl;
const Vector<label>& directions = mesh().geometricD();
scalar thickness = 0.0;
for (direction dir=0; dir<directions.nComponents; dir++)
{
if (directions[dir] == -1)
{
boundBox bb(mesh().points(), false);
thickness = bb.span()[dir];
thickness = mesh().bounds().span()[dir];
break;
}
}

3
src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C
View File

@ -41,8 +41,7 @@ namespace Foam
void Foam::IDDESDelta::calcDelta()
{
const Vector<label>& directions = mesh().directions();
label nD = (directions.nComponents + cmptSum(directions))/2;
label nD = mesh().nGeometricD();
// initialise hwn as wall distance
volScalarField hwn = wallDist(mesh()).y();