Merge remote branch 'OpenCFD/master' into olesenm

- use GPL version 3
This commit is contained in:
Mark Olesen
2010-05-18 15:00:40 +02:00
9 changed files with 623 additions and 91 deletions

View File

@ -158,7 +158,6 @@ void Foam::createShellMesh::calcPointRegions
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::createShellMesh::createShellMesh
@ -184,7 +183,6 @@ Foam::createShellMesh::createShellMesh
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::createShellMesh::setRefinement
(
const pointField& thickness,

View File

@ -130,6 +130,7 @@ Usage
#include "volFields.H"
#include "surfaceFields.H"
#include "cyclicPolyPatch.H"
#include "syncTools.H"
using namespace Foam;
@ -595,6 +596,237 @@ void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
}
//XXXXXXXXX
label findUncoveredPatchFace
(
const fvMesh& mesh,
const UIndirectList<label>& extrudeMeshFaces,// mesh faces that are extruded
const label meshEdgeI // mesh edge
)
{
// Make set of extruded faces.
labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
forAll(extrudeMeshFaces, i)
{
extrudeFaceSet.insert(extrudeMeshFaces[i]);
}
label patchI = -1;
const labelList& eFaces = mesh.edgeFaces()[meshEdgeI];
forAll(eFaces, i)
{
label faceI = eFaces[i];
if (!mesh.isInternalFace(faceI) && !extrudeFaceSet.found(faceI))
{
patchI = mesh.boundaryMesh().whichPatch(faceI);
break;
}
}
return patchI;
}
// Count the number of faces in patches that need to be created
void countExtrudePatches
(
const fvMesh& mesh,
const primitiveFacePatch& extrudePatch,
const label nZones,
const labelList& zoneID,
const labelList& extrudeMeshFaces,
const labelList& extrudeMeshEdges,
labelList& zoneSidePatch,
labelList& zoneZonePatch
)
{
const labelListList& edgeFaces = extrudePatch.edgeFaces();
forAll(edgeFaces, edgeI)
{
const labelList& eFaces = edgeFaces[edgeI];
if (eFaces.size() == 2)
{
label zone0 = zoneID[eFaces[0]];
label zone1 = zoneID[eFaces[1]];
if (zone0 != zone1)
{
label minZone = min(zone0,zone1);
label maxZone = max(zone0,zone1);
zoneZonePatch[minZone*nZones+maxZone]++;
}
}
else
{
// Check whether we are on a mesh edge with external patches. If
// so choose any uncovered one. If none found put face in
// undetermined zone 'side' patch
label patchI = findUncoveredPatchFace
(
mesh,
UIndirectList<label>(extrudeMeshFaces, eFaces),
extrudeMeshEdges[edgeI]
);
if (patchI == -1)
{
// Determine the min zone of all connected zones.
label minZone = zoneID[eFaces[0]];
for (label i = 1; i < eFaces.size(); i++)
{
minZone = min(minZone, zoneID[eFaces[i]]);
}
zoneSidePatch[minZone]++;
}
}
}
Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneSidePatch);
Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneZonePatch);
}
bool lessThan(const point& x, const point& y)
{
for (direction dir = 0; dir < point::nComponents; dir++)
{
if (x[dir] < y[dir]) return true;
if (x[dir] > y[dir]) return false;
}
return false;
}
class minEqVectorListOp
{
public:
void operator()(List<vector>& x, const List<vector>& y) const
{
if (y.size())
{
if (x.size())
{
forAll(x, i)
{
if (lessThan(y[i], x[i]))
{
x[i] = y[i];
}
}
}
else
{
x = y;
}
}
}
};
// Constrain&sync normals on points that are on coupled patches.
void constrainCoupledNormals
(
const fvMesh& mesh,
const primitiveFacePatch& extrudePatch,
const labelList& regionToPoint,
vectorField& regionNormals
)
{
// Invert regionToPoint to create pointToRegions.
labelListList pointToRegions
(
invertOneToMany
(
extrudePatch.nPoints(),
regionToPoint
)
);
// Sort acc. to region so (hopefully) coupled points will do the same.
forAll(pointToRegions, pointI)
{
sort(pointToRegions[pointI]);
}
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Constrain displacement on cyclic patches
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Note: bit contentious to always do this on cyclic - should user use
// different patch type, e.g. 'cyclicSlip' instead?
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<cyclicPolyPatch>(pp))
{
forAll(pp.meshPoints(), pointI)
{
Map<label>::const_iterator fnd =
extrudePatch.meshPointMap().find
(
pp.meshPoints()[pointI]
);
if (fnd != extrudePatch.meshPointMap().end())
{
// fnd() is a point on this cyclic.
const vector& cycNormal = pp.pointNormals()[pointI];
const labelList& pRegions = pointToRegions[fnd()];
forAll(pRegions, i)
{
// Remove cyclic normal component.
vector& regionNormal = regionNormals[pRegions[i]];
regionNormal -= (regionNormal&cycNormal)*cycNormal;
}
}
}
}
}
// Synchronise regionNormals
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Re-work regionNormals into multiple normals per point
List<List<vector> > pointNormals(mesh.nPoints());
forAll(pointToRegions, pointI)
{
const labelList& pRegions = pointToRegions[pointI];
label meshPointI = extrudePatch.meshPoints()[pointI];
List<vector>& pNormals = pointNormals[meshPointI];
pNormals.setSize(pRegions.size());
forAll(pRegions, i)
{
pNormals[i] = regionNormals[pRegions[i]];
}
}
// Synchronise
syncTools::syncPointList
(
mesh,
pointNormals,
minEqVectorListOp(),
List<vector>(), // nullValue
false // applySeparation
);
// Re-work back into regionNormals
forAll(pointToRegions, pointI)
{
const labelList& pRegions = pointToRegions[pointI];
label meshPointI = extrudePatch.meshPoints()[pointI];
const List<vector>& pNormals = pointNormals[meshPointI];
forAll(pRegions, i)
{
regionNormals[pRegions[i]] = pNormals[i];
}
}
}
//XXXXXXXXX
tmp<pointField> calcOffset
(
const primitiveFacePatch& extrudePatch,
@ -739,6 +971,16 @@ int main(int argc, char *argv[])
<< endl;
// Determine corresponding mesh edges
const labelList extrudeMeshEdges
(
extrudePatch.meshEdges
(
mesh.edges(),
mesh.pointEdges()
)
);
// Check whether the zone is internal or external faces to determine
@ -772,7 +1014,7 @@ int main(int argc, char *argv[])
Info<< "FaceZone " << fz.name() << " has boundary faces" << endl;
}
}
Info<< endl;
@ -866,36 +1108,18 @@ int main(int argc, char *argv[])
labelList zoneSidePatch(faceZones.size(), 0);
labelList zoneZonePatch(faceZones.size()*faceZones.size(), 0);
forAll(edgeFaces, edgeI)
{
const labelList& eFaces = edgeFaces[edgeI];
if (eFaces.size() == 2)
{
label zone0 = zoneID[eFaces[0]];
label zone1 = zoneID[eFaces[1]];
countExtrudePatches
(
mesh,
extrudePatch,
faceZones.size(),
zoneID,
extrudeMeshFaces,
extrudeMeshEdges,
if (zone0 != zone1)
{
label minZone = min(zone0,zone1);
label maxZone = max(zone0,zone1);
zoneZonePatch[minZone*faceZones.size()+maxZone]++;
}
}
else
{
// Determine the min zone of all connected zones.
label minZone = zoneID[eFaces[0]];
for (label i = 1; i < eFaces.size(); i++)
{
minZone = min(minZone, zoneID[eFaces[i]]);
}
zoneSidePatch[minZone]++;
}
}
Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneSidePatch);
Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneZonePatch);
zoneSidePatch,
zoneZonePatch
);
// Now check which patches to add.
Info<< "Adding patches for edges on zones:" << nl << nl
@ -980,6 +1204,7 @@ int main(int argc, char *argv[])
// Is edge an non-manifold edge
PackedBoolList nonManifoldEdge(extrudePatch.nEdges());
// Note: logic has to be same as in countExtrudePatches.
forAll(edgeFaces, edgeI)
{
const labelList& eFaces = edgeFaces[edgeI];
@ -1004,63 +1229,29 @@ int main(int argc, char *argv[])
}
else
{
ePatches.setSize(eFaces.size());
forAll(eFaces, i)
label patchI = findUncoveredPatchFace
(
mesh,
UIndirectList<label>(extrudeMeshFaces, eFaces),
extrudeMeshEdges[edgeI]
);
if (patchI != -1)
{
ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
ePatches.setSize(eFaces.size(), patchI);
}
else
{
ePatches.setSize(eFaces.size());
forAll(eFaces, i)
{
ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
}
}
nonManifoldEdge[edgeI] = 1;
}
}
// Override constraint types
{
const edgeList& extrudeEdges = extrudePatch.edges();
const labelList& extrudeMeshPts = extrudePatch.meshPoints();
// Map from mesh edge to local patch edge index
EdgeMap<label> extrudeMeshEdges(extrudePatch.nEdges());
forAll(extrudeEdges, edgeI)
{
if (edgeFaces[edgeI].size() == 1)
{
const edge& e = extrudeEdges[edgeI];
const edge meshE(extrudeMeshPts[e[0]], extrudeMeshPts[e[1]]);
extrudeMeshEdges.insert(meshE, edgeI);
}
}
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (polyPatch::constraintType(pp.type()))
{
const edgeList& edges = pp.edges();
forAll(edges, ppEdgeI)
{
const edge& e = edges[ppEdgeI];
edge meshE(pp.meshPoints()[e[0]], pp.meshPoints()[e[1]]);
EdgeMap<label>::const_iterator iter = extrudeMeshEdges.find
(
meshE
);
if (iter != extrudeMeshEdges.end())
{
label extrudeEdgeI = iter();
extrudeEdgePatches[extrudeEdgeI] = labelList
(
edgeFaces[extrudeEdgeI].size(),
patchI
);
}
}
}
}
}
// Assign point regions
@ -1075,6 +1266,7 @@ int main(int argc, char *argv[])
(
extrudePatch,
nonManifoldEdge,
pointRegions,
regionPoints
);
@ -1096,6 +1288,18 @@ int main(int argc, char *argv[])
}
regionNormals /= mag(regionNormals);
// Constrain&sync normals on points that are on coupled patches.
constrainCoupledNormals
(
mesh,
extrudePatch,
regionPoints,
regionNormals
);
// For debugging: dump hedgehog plot of normals
{
OFstream str(runTime.path()/"regionNormals.obj");
@ -1112,7 +1316,7 @@ int main(int argc, char *argv[])
meshTools::writeOBJ(str, pt);
vertI++;
meshTools::writeOBJ(str, pt+0.01*regionNormals[region]);
meshTools::writeOBJ(str, pt+thickness*regionNormals[region]);
vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl;
}

View File

@ -346,7 +346,7 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
{
ParcelType& p = static_cast<ParcelType&>(*this);
// Invoke poost-processing mdoel
// Invoke post-processing model
td.cloud().postProcessing().postPatch(p, patchI);
// Invoke surface film model

View File

@ -212,7 +212,12 @@ Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
}
case pbNumber:
{
nP = massTotal_/(rho*volumeTotal_*parcels);
nP = massTotal_/(rho*volumeTotal_);
break;
}
case pbFixed:
{
nP = nParticlesFixed_;
break;
}
default:
@ -285,6 +290,7 @@ Foam::InjectionModel<CloudType>::InjectionModel(CloudType& owner)
nInjections_(0),
parcelsAddedTotal_(0),
parcelBasis_(pbNumber),
nParticlesFixed_(0.0),
time0_(0.0),
timeStep0_(0.0)
{
@ -310,6 +316,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
nInjections_(0),
parcelsAddedTotal_(0),
parcelBasis_(pbNumber),
nParticlesFixed_(0.0),
time0_(owner.db().time().value()),
timeStep0_(0.0)
{
@ -320,6 +327,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
<< endl;
const word parcelBasisType = coeffDict_.lookup("parcelBasisType");
if (parcelBasisType == "mass")
{
parcelBasis_ = pbMass;
@ -328,6 +336,16 @@ Foam::InjectionModel<CloudType>::InjectionModel
{
parcelBasis_ = pbNumber;
}
else if (parcelBasisType == "fixed")
{
parcelBasis_ = pbFixed;
Info<< " Choosing nParticles to be a fixed value, massTotal "
<< "variable now does not determine anything."
<< endl;
nParticlesFixed_ = readScalar(coeffDict_.lookup("nParticles"));
}
else
{
FatalErrorIn
@ -338,7 +356,7 @@ Foam::InjectionModel<CloudType>::InjectionModel
"CloudType&, "
"const word&"
")"
)<< "parcelBasisType must be either 'number' or 'mass'" << nl
)<< "parcelBasisType must be either 'number', 'mass' or 'fixed'" << nl
<< exit(FatalError);
}

View File

@ -73,7 +73,8 @@ public:
enum parcelBasis
{
pbNumber,
pbMass
pbMass,
pbFixed
};
@ -137,6 +138,10 @@ protected:
//- Parcel basis enumeration
parcelBasis parcelBasis_;
//- nParticles to assign to parcels when the 'fixed' basis
// is selected
scalar nParticlesFixed_;
//- Continuous phase time at start of injection time step [s]
scalar time0_;

View File

@ -76,9 +76,6 @@ class ManualInjection
//- Parcel size PDF model
const autoPtr<pdfs::pdf> parcelPDF_;
//- Number of particles represented by each parcel
scalar nParticlesPerParcel_;
protected:

View File

@ -43,6 +43,7 @@ $(kqRWallFunctions)/kqRWallFunction/kqRWallFunctionFvPatchFields.C
derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C
derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C

View File

@ -0,0 +1,163 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "fixedShearStressFvPatchVectorField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
fixedShearStressFvPatchVectorField::
fixedShearStressFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchVectorField(p, iF),
tau0_(vector::zero)
{}
fixedShearStressFvPatchVectorField::
fixedShearStressFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchVectorField(p, iF),
tau0_(dict.lookupOrDefault<vector>("tau", vector::zero))
{
fvPatchField<vector>::operator=(patchInternalField());
}
fixedShearStressFvPatchVectorField::
fixedShearStressFvPatchVectorField
(
const fixedShearStressFvPatchVectorField& ptf,
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchVectorField(ptf, p, iF, mapper),
tau0_(ptf.tau0_)
{}
fixedShearStressFvPatchVectorField::
fixedShearStressFvPatchVectorField
(
const fixedShearStressFvPatchVectorField& ptf
)
:
fixedValueFvPatchVectorField(ptf),
tau0_(ptf.tau0_)
{}
fixedShearStressFvPatchVectorField::
fixedShearStressFvPatchVectorField
(
const fixedShearStressFvPatchVectorField& ptf,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchVectorField(ptf, iF),
tau0_(ptf.tau0_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void fixedShearStressFvPatchVectorField::updateCoeffs()
{
if (updated())
{
return;
}
const label patchI = patch().index();
const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
const vectorField Ui = Uw.patchInternalField();
vector tauHat = tau0_/mag(tau0_);
const scalarField& ry = patch().deltaCoeffs();
scalarField nuEffw = rasModel.nuEff()().boundaryField()[patchI];
vectorField UwUpdated =
tauHat*(tauHat & (tau0_*(1.0/(ry*nuEffw)) + Ui));
operator==(UwUpdated);
if (debug)
{
vectorField nHat = this->patch().nf();
volSymmTensorField Reff = rasModel.devReff();
Info<< "tau : " << (nHat & Reff.boundaryField()[patchI])() << endl;
}
fixedValueFvPatchVectorField::updateCoeffs();
}
void fixedShearStressFvPatchVectorField::write(Ostream& os) const
{
fixedValueFvPatchVectorField::write(os);
os.writeKeyword("tau") << tau0_ << token::END_STATEMENT << nl;
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchVectorField,
fixedShearStressFvPatchVectorField
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

View File

@ -0,0 +1,146 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::fixedShearStressFvPatchVectorField
Description
Set a constant shear stress as tau0 = -nuEff dU/dn.
SourceFiles
fixedShearStressFvPatchVectorField.C
\*---------------------------------------------------------------------------*/
#ifndef fixedShearStressFvPatchVectorField_H
#define fixedShearStressFvPatchVectorField_H
#include "fvPatchFields.H"
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
/*---------------------------------------------------------------------------*\
Class fixedShearStressFvPatchVectorField Declaration
\*---------------------------------------------------------------------------*/
class fixedShearStressFvPatchVectorField
:
public fixedValueFvPatchVectorField
{
// Private data
//- Constant shear stress
vector tau0_;
public:
//- Runtime type information
TypeName("fixedShearStress");
// Constructors
//- Construct from patch and internal field
fixedShearStressFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&
);
//- Construct from patch, internal field and dictionary
fixedShearStressFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const dictionary&
);
//- Construct by mapping given
fixedShearStressFvPatchVectorField
(
const fixedShearStressFvPatchVectorField&,
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct as copy
fixedShearStressFvPatchVectorField
(
const fixedShearStressFvPatchVectorField&
);
//- Construct and return a clone
virtual tmp<fvPatchVectorField> clone() const
{
return tmp<fvPatchVectorField>
(
new fixedShearStressFvPatchVectorField(*this)
);
}
//- Construct as copy setting internal field reference
fixedShearStressFvPatchVectorField
(
const fixedShearStressFvPatchVectorField&,
const DimensionedField<vector, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchVectorField> clone
(
const DimensionedField<vector, volMesh>& iF
) const
{
return tmp<fvPatchVectorField>
(
new fixedShearStressFvPatchVectorField(*this, iF)
);
}
// Member functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //