diff --git a/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.C b/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.C
index 098b4282a1b3c6659a90d9ea99856237eb2fe054..c7df352eb91e22343703902f58218e9fee335d40 100644
--- a/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.C
+++ b/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.C
@@ -43,14 +43,45 @@ namespace fv
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+Foam::volScalarField& Foam::fv::atmPlantCanopyTurbSource::getOrReadField
+(
+ const word& fieldName
+) const
+{
+ auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+ if (!ptr)
+ {
+ ptr = new volScalarField
+ (
+ IOobject
+ (
+ fieldName,
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_
+ );
+ mesh_.objectRegistry::store(ptr);
+ }
+
+ return *ptr;
+}
+
+
Foam::tmp<Foam::volScalarField::Internal>
Foam::fv::atmPlantCanopyTurbSource::calcPlantCanopyTerm
(
const volVectorField::Internal& U
) const
{
+ const volScalarField& Cd = getOrReadField(CdName_);
+ const volScalarField& LAD = getOrReadField(LADname_);
+
// (SP:Eq. 42)
- return 12.0*Foam::sqrt(Cmu_)*plantCd_()*leafAreaDensity_()*mag(U);
+ return 12.0*Foam::sqrt(Cmu_)*Cd()*LAD()*mag(U);
}
@@ -65,36 +96,15 @@ Foam::fv::atmPlantCanopyTurbSource::atmPlantCanopyTurbSource
)
:
fv::cellSetOption(sourceName, modelType, dict, mesh),
- isEpsilon_(false),
- rhoName_(coeffs_.getOrDefault<word>("rho", "rho")),
+ isEpsilon_(true),
Cmu_(Zero),
C1_(Zero),
C2_(Zero),
- plantCd_
- (
- IOobject
- (
- "plantCd",
- mesh.time().timeName(),
- mesh,
- IOobject::MUST_READ,
- IOobject::AUTO_WRITE
- ),
- mesh
- ),
- leafAreaDensity_
- (
- IOobject
- (
- "leafAreaDensity",
- mesh.time().timeName(),
- mesh,
- IOobject::MUST_READ,
- IOobject::AUTO_WRITE
- ),
- mesh
- )
+ CdName_(),
+ LADname_()
{
+ read(dict);
+
const auto* turbPtr =
mesh_.findObject<turbulenceModel>
(
@@ -113,9 +123,8 @@ Foam::fv::atmPlantCanopyTurbSource::atmPlantCanopyTurbSource
tmp<volScalarField> tepsilon = turbPtr->epsilon();
tmp<volScalarField> tomega = turbPtr->omega();
- if (tepsilon.is_reference())
+ if (!tepsilon.isTmp())
{
- isEpsilon_ = true;
fieldNames_[0] = tepsilon().name();
const dictionary& turbDict = turbPtr->coeffDict();
@@ -123,7 +132,7 @@ Foam::fv::atmPlantCanopyTurbSource::atmPlantCanopyTurbSource
C1_.read("C1", turbDict);
C2_.read("C2", turbDict);
}
- else if (tomega.is_reference())
+ else if (!tomega.isTmp())
{
isEpsilon_ = false;
fieldNames_[0] = tomega().name();
@@ -214,4 +223,21 @@ void Foam::fv::atmPlantCanopyTurbSource::addSup
}
+bool Foam::fv::atmPlantCanopyTurbSource::read(const dictionary& dict)
+{
+ if (!fv::cellSetOption::read(dict))
+ {
+ return false;
+ }
+
+ CdName_ = dict.getOrDefault<word>("Cd", "Cd");
+ LADname_ = dict.getOrDefault<word>("LAD", "LAD");
+
+ (void) getOrReadField(CdName_);
+ (void) getOrReadField(LADname_);
+
+ return true;
+}
+
+
// ************************************************************************* //
diff --git a/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.H b/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.H
index e6a053ee531cf5b9cea3cc17cef2728910847149..5d255efe326d63d02053e289033db4831071a216 100644
--- a/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.H
+++ b/src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 ENERCON GmbH
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -36,15 +36,16 @@ Description
Corrections applied to either of the below, if exist:
\verbatim
- epsilon | Turbulent kinetic energy dissipation rate [m2/s3]
+ epsilon | Turbulent kinetic energy dissipation rate [m^2/s^3]
omega | Specific dissipation rate [1/s]
\endverbatim
Required fields:
\verbatim
- epsilon/omega | Dissipation rate OR Spec. dissipation rate [m2/s3]/[1/s]
- plantCd | Plant canopy drag coefficient [-]
- leafAreaDensity | Leaf area density [1/m]
+ Cd | Canopy drag coefficient [-]
+ LAD | Leaf area density [m^2/m^3]
+ epsilon | Turbulent kinetic energy dissipation rate [m^2/s^3]
+ omega | Specific dissipation rate [1/s]
\endverbatim
References:
@@ -61,28 +62,24 @@ Usage
\verbatim
atmPlantCanopyTurbSource1
{
- // Mandatory entries (unmodifiable)
- type atmPlantCanopyTurbSource;
-
- atmPlantCanopyTurbSourceCoeffs
- {
- // Mandatory (inherited) entries (unmodifiable)
- selectionMode all;
+ // Mandatory entries
+ type atmPlantCanopyTurbSource;
- // Optional entries (unmodifiable)
- rho rho;
- }
+ // Optional entries
+ Cd <word>;
+ LAD <word>;
- // Optional (inherited) entries
+ // Inherited entries
...
}
\endverbatim
where the entries mean:
\table
- Property | Description | Type | Req'd | Dflt
- type | Type name: atmPlantCanopyTurbSource | word | yes | -
- rho | Name of density field | word | no | rho
+ Property | Description | Type | Reqd | Deflt
+ type | Type name: atmPlantCanopyTurbSource | word | yes | -
+ Cd | Name of operand canopy drag coefficient field | word | no | Cd
+ LAD | Name of operand leaf area density field | word | no | LAD
\endtable
The inherited entries are elaborated in:
@@ -109,7 +106,7 @@ namespace fv
{
/*---------------------------------------------------------------------------*\
- Class atmPlantCanopyTurbSource Declaration
+ Class atmPlantCanopyTurbSource Declaration
\*---------------------------------------------------------------------------*/
class atmPlantCanopyTurbSource
@@ -121,26 +118,24 @@ class atmPlantCanopyTurbSource
//- Internal flag to determine the working field is epsilon or omega
Switch isEpsilon_;
- //- Name of density field
- const word rhoName_;
-
//- Required turbulence model coefficients (copied from turb model)
dimensionedScalar Cmu_;
dimensionedScalar C1_;
dimensionedScalar C2_;
+ //- Name of operand canopy drag coefficient field
+ word CdName_;
- // Fields
-
- //- Plant canopy drag coefficient field [-]
- volScalarField plantCd_;
-
- //- Leaf area density field [1/m]
- volScalarField leafAreaDensity_;
+ //- Name of operand leaf area density field
+ word LADname_;
// Private Member Functions
+ //- Return requested field from the object registry
+ //- or read+register the field to the object registry
+ volScalarField& getOrReadField(const word& fieldName) const;
+
//- Return the modifier for plant canopy effects
tmp<volScalarField::Internal> calcPlantCanopyTerm
(
@@ -221,11 +216,8 @@ public:
const label fieldi
);
- //- Read source dictionary (effectively no-op)
- virtual bool read(const dictionary& dict)
- {
- return true;
- }
+ //- Read source dictionary
+ virtual bool read(const dictionary& dict);
};
diff --git a/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C b/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
index 15f448ff5735ac57a8377d8999c8b6ca2a40c5b0..158ac966478d005190c5c422b87997a49c7f6efb 100644
--- a/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
+++ b/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
@@ -42,6 +42,36 @@ namespace fv
}
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+Foam::volScalarField& Foam::fv::atmPlantCanopyUSource::getOrReadField
+(
+ const word& fieldName
+) const
+{
+ auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+ if (!ptr)
+ {
+ ptr = new volScalarField
+ (
+ IOobject
+ (
+ fieldName,
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_
+ );
+ mesh_.objectRegistry::store(ptr);
+ }
+
+ return *ptr;
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::fv::atmPlantCanopyUSource::atmPlantCanopyUSource
@@ -53,32 +83,11 @@ Foam::fv::atmPlantCanopyUSource::atmPlantCanopyUSource
)
:
fv::cellSetOption(sourceName, modelType, dict, mesh),
- rhoName_(coeffs_.getOrDefault<word>("rho", "rho")),
- plantCd_
- (
- IOobject
- (
- "plantCd",
- mesh.time().timeName(),
- mesh,
- IOobject::MUST_READ,
- IOobject::AUTO_WRITE
- ),
- mesh
- ),
- leafAreaDensity_
- (
- IOobject
- (
- "leafAreaDensity",
- mesh.time().timeName(),
- mesh,
- IOobject::MUST_READ,
- IOobject::AUTO_WRITE
- ),
- mesh
- )
+ CdName_(),
+ LADname_()
{
+ read(dict);
+
fieldNames_.resize(1, "U");
fv::option::resetApplied();
@@ -95,13 +104,17 @@ void Foam::fv::atmPlantCanopyUSource::addSup
const label fieldi
)
{
- const volVectorField& U = eqn.psi();
-
- if (V_ > VSMALL)
+ if (V_ < VSMALL)
{
- // (SP:Eq. 42)
- eqn -= fvm::Sp(plantCd_*leafAreaDensity_*mag(U), U);
+ return;
}
+
+ const volVectorField& U = eqn.psi();
+ const volScalarField& Cd = getOrReadField(CdName_);
+ const volScalarField& LAD = getOrReadField(LADname_);
+
+ // (SP:Eq. 42), (BSG:Eq. 7)
+ eqn -= fvm::Sp(Cd*LAD*mag(U), U);
}
@@ -112,12 +125,16 @@ void Foam::fv::atmPlantCanopyUSource::addSup
const label fieldi
)
{
- const volVectorField& U = eqn.psi();
-
- if (V_ > VSMALL)
+ if (V_ < VSMALL)
{
- eqn -= fvm::Sp(rho*plantCd_*leafAreaDensity_*mag(U), U);
+ return;
}
+
+ const volVectorField& U = eqn.psi();
+ const volScalarField& Cd = getOrReadField(CdName_);
+ const volScalarField& LAD = getOrReadField(LADname_);
+
+ eqn -= fvm::Sp(rho*Cd*LAD*mag(U), U);
}
@@ -129,12 +146,33 @@ void Foam::fv::atmPlantCanopyUSource::addSup
const label fieldi
)
{
+ if (V_ < VSMALL)
+ {
+ return;
+ }
+
const volVectorField& U = eqn.psi();
+ const volScalarField& Cd = getOrReadField(CdName_);
+ const volScalarField& LAD = getOrReadField(LADname_);
+
+ eqn -= fvm::Sp(alpha*rho*Cd*LAD*mag(U), U);
+}
+
- if (V_ > VSMALL)
+bool Foam::fv::atmPlantCanopyUSource::read(const dictionary& dict)
+{
+ if (!fv::cellSetOption::read(dict))
{
- eqn -= fvm::Sp(alpha*rho*plantCd_*leafAreaDensity_*mag(U), U);
+ return false;
}
+
+ CdName_ = dict.getOrDefault<word>("Cd", "Cd");
+ LADname_ = dict.getOrDefault<word>("LAD", "LAD");
+
+ (void) getOrReadField(CdName_);
+ (void) getOrReadField(LADname_);
+
+ return true;
}
diff --git a/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.H b/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.H
index 035204be7e14d70d099fa4e1aa10b28b02f3041a..842e80e6e720796c77f266ca3faf12f09ecda550 100644
--- a/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.H
+++ b/src/atmosphericModels/fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 ENERCON GmbH
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -31,7 +31,7 @@ Group
grpFvOptionsSources
Description
- Applies sources on velocity, i.e. \c U, to incorporate effects
+ Applies sources on velocity (i.e. \c U) to incorporate effects
of plant canopy for atmospheric boundary layer modelling.
Corrections applied to:
@@ -42,8 +42,8 @@ Description
Required fields:
\verbatim
U | Velocity [m/s]
- plantCd | Plant canopy drag coefficient [-]
- leafAreaDensity | Leaf area density [1/m]
+ Cd | Plant canopy drag coefficient [-]
+ LAD | Leaf area density [1/m]
\endverbatim
References:
@@ -53,6 +53,13 @@ Description
Modification of two-equation models to account for plant drag.
Boundary-Layer Meteorology, 121(2), 229-266.
DOI:10.1007/s10546-006-9073-5
+
+ Governing equations (tag:BSG):
+ Brozovsky, J., Simonsen, A., & Gaitani, N. (2021).
+ Validation of a CFD model for the evaluation of urban microclimate
+ at high latitudes: A case study in Trondheim, Norway.
+ Building and Environment, 205, 108175.
+ DOI:10.1016/j.buildenv.2021.108175
\endverbatim
Usage
@@ -60,28 +67,24 @@ Usage
\verbatim
atmPlantCanopyUSource1
{
- // Mandatory entries (unmodifiable)
- type atmPlantCanopyUSource;
-
- atmPlantCanopyUSourceCoeffs
- {
- // Mandatory (inherited) entries (unmodifiable)
- selectionMode all;
+ // Mandatory entries
+ type atmPlantCanopyUSource;
- // Optional entries (unmodifiable)
- rho rho;
- }
+ // Optional entries
+ Cd <word>;
+ LAD <word>;
- // Optional (inherited) entries
+ // Inherited entries
...
}
\endverbatim
where the entries mean:
\table
- Property | Description | Type | Req'd | Dflt
- type | Type name: atmPlantCanopyUSource | word | yes | -
- rho | Name of density field | word | no | rho
+ Property | Description | Type | Reqd | Deflt
+ type | Type name: atmPlantCanopyUSource | word | yes | -
+ Cd | Name of operand canopy drag coefficient field | word | no | Cd
+ LAD | Name of operand leaf area density field | word | no | LAD
\endtable
The inherited entries are elaborated in:
@@ -108,7 +111,7 @@ namespace fv
{
/*---------------------------------------------------------------------------*\
- Class atmPlantCanopyUSource Declaration
+ Class atmPlantCanopyUSource Declaration
\*---------------------------------------------------------------------------*/
class atmPlantCanopyUSource
@@ -117,17 +120,18 @@ class atmPlantCanopyUSource
{
// Private Data
- //- Name of density field
- const word rhoName_;
+ //- Name of operand canopy drag coefficient field
+ word CdName_;
+ //- Name of operand leaf area density field
+ word LADname_;
- // Fields
- //- Plant canopy drag coefficient field [-]
- volScalarField plantCd_;
+ // Private Member Functions
- //- Leaf area density field [1/m]
- volScalarField leafAreaDensity_;
+ //- Return requested field from the object registry
+ //- or read+register the field to the object registry
+ volScalarField& getOrReadField(const word& fieldName) const;
public:
@@ -180,11 +184,8 @@ public:
const label fieldi
);
- //- Read source dictionary (effectively no-op)
- virtual bool read(const dictionary& dict)
- {
- return true;
- }
+ //- Read source dictionary
+ virtual bool read(const dictionary& dict);
};
diff --git a/src/atmosphericModels/turbulenceModels/RAS/kL/kL.C b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.C
index 94e393c2eef709c06046a989648344a544811998..0b06f04c14df0a864f007bab479adb1b03485c29 100644
--- a/src/atmosphericModels/turbulenceModels/RAS/kL/kL.C
+++ b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.C
@@ -67,10 +67,8 @@ tmp<volScalarField> kL<BasicTurbulenceModel>::nutPrime() const
template<class BasicTurbulenceModel>
tmp<volScalarField> kL<BasicTurbulenceModel>::epsilonCanopy() const
{
- const auto* CdPtr =
- this->mesh_.template findObject<volScalarField>("plantCd");
- const auto* LADPtr =
- this->mesh_.template findObject<volScalarField>("leafAreaDensity");
+ const auto* CdPtr = this->mesh_.template findObject<volScalarField>("Cd");
+ const auto* LADPtr = this->mesh_.template findObject<volScalarField>("LAD");
const volVectorField& U = this->U_;
if (CdPtr && LADPtr)
diff --git a/src/atmosphericModels/turbulenceModels/RAS/kL/kL.H b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.H
index 870b38ac694171b3f9c96a22ebf39cc2199f6dd0..46caa1b5b431188196aebe41854b4eacc82cc297 100644
--- a/src/atmosphericModels/turbulenceModels/RAS/kL/kL.H
+++ b/src/atmosphericModels/turbulenceModels/RAS/kL/kL.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2021 OpenCFD Ltd.
+ Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -104,8 +104,8 @@ Usage
\heading Input fields (optional)
\plaintable
canopyHeight | Canopy height [m]
- plantCd | Plant canopy drag coefficient [-]
- leafAreaDensity | Leaf area density [1/m]
+ Cd | Plant canopy drag coefficient [-]
+ LAD | Leaf area density [1/m]
Rt | Turbulent Richardson number [-]
L | Characteristic length scale [m]
\endplaintable
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 5f3a5d7d252a32300437b840970fb996da99725d..4f974d5d2f2d93cfde5cf0513b6a96503d59f9cd 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -224,6 +224,7 @@ $(faceZoneSources)/setToFaceZone/setToFaceZone.C
$(faceZoneSources)/setAndNormalToFaceZone/setAndNormalToFaceZone.C
$(faceZoneSources)/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C
$(faceZoneSources)/planeToFaceZone/planeToFaceZone.C
+$(faceZoneSources)/cellToFaceZone/cellToFaceZone.C
cellZoneSources = topoSet/cellZoneSources
$(cellZoneSources)/topoSetCellZoneSource/topoSetCellZoneSource.C
diff --git a/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.C b/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.C
index b3401415a0c5d5310ea9ed40a113437b699313b3..9e88902e6c4ac7a68f9d4dbbc338703d6f546115 100644
--- a/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.C
+++ b/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.C
@@ -48,11 +48,12 @@ namespace Foam
Foam::topoSetSource::addToUsageTable Foam::cellToFace::usage_
(
cellToFace::typeName,
- "\n Usage: cellToFace <cellSet> all|both\n\n"
+ "\n Usage: cellToFace <cellSet> all|both|outside\n\n"
" Select -all : all faces of cells in the cellSet\n"
" -both: faces where both neighbours are in the cellSet\n\n"
);
+
const Foam::Enum
<
Foam::cellToFace::cellAction
@@ -61,6 +62,7 @@ Foam::cellToFace::cellActionNames_
({
{ cellAction::ALL, "all" },
{ cellAction::BOTH, "both" },
+ { cellAction::OUTSIDE, "outside" }
});
@@ -148,6 +150,65 @@ void Foam::cellToFace::combine
}
}
}
+ else if (option_ == OUTSIDE)
+ {
+ // Add all faces where only one neighbour is in set.
+
+ const label nInt = mesh_.nInternalFaces();
+ const labelList& own = mesh_.faceOwner();
+ const labelList& nei = mesh_.faceNeighbour();
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+
+ // Check all internal faces
+ for (label facei = 0; facei < nInt; ++facei)
+ {
+ if (cellLabels.found(own[facei]) != cellLabels.found(nei[facei]))
+ {
+ addOrDelete(set, facei, add);
+ }
+ }
+
+
+ // Get coupled cell status
+ boolList neiInSet(mesh_.nBoundaryFaces(), false);
+
+ for (const polyPatch& pp : patches)
+ {
+ if (pp.coupled())
+ {
+ label facei = pp.start();
+ forAll(pp, i)
+ {
+ neiInSet[facei-nInt] = cellLabels.found(own[facei]);
+ ++facei;
+ }
+ }
+ }
+ syncTools::swapBoundaryFaceList(mesh_, neiInSet);
+
+
+ // Check all boundary faces
+ for (const polyPatch& pp : patches)
+ {
+ label facei = pp.start();
+ forAll(pp, i)
+ {
+ if (cellLabels.found(own[facei]) != neiInSet[facei-nInt])
+ {
+ addOrDelete(set, facei, add);
+ }
+ ++facei;
+ }
+ }
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Selected option is not available"
+ << ", option: " << cellActionNames_[option_]
+ << exit(FatalError);
+ }
}
diff --git a/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.H b/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.H
index c9ce52b0fdf40193828da9554a632ab0df61eb45..d1c5638a597043c357a99e9bca83eb15fc339b48 100644
--- a/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.H
+++ b/src/meshTools/topoSet/faceSources/cellToFace/cellToFace.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011 OpenFOAM Foundation
- Copyright (C) 2018-2020 OpenCFD Ltd.
+ Copyright (C) 2018-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -33,8 +33,8 @@ Description
Operands:
\table
Operand | Type | Location
- input | cellSet(s) | $FOAM_CASE/constant/polyMesh/sets/\<set\>
- output | faceSet | $FOAM_CASE/constant/polyMesh/sets/\<set\>
+ input | cellSet(s) | constant/polyMesh/sets/\<set\>
+ output | faceSet | constant/polyMesh/sets/\<set\>
\endtable
Usage
@@ -68,7 +68,7 @@ Usage
where the entries mean:
\table
- Property | Description | Type | Req'd | Dflt
+ Property | Description | Type | Reqd | Deflt
name | Name of faceSet | word | yes | -
type | Type name: faceSet | word | yes | -
action | Action applied on faces - see below | word | yes | -
@@ -87,18 +87,20 @@ Usage
\verbatim
all | All faces of cells in the cellSet
both | Faces where both neighbours are in the cellSet
+ outside | Faces with only one neighbour in the cellSet
\endverbatim
Options for the conditional mandatory entries:
\verbatim
- Entry | Description | Type | Req'd | Dflt
- sets | Names of input cellSets | wordList | cond'l | -
- set | Name of input cellSet | word | cond'l | -
+ Entry | Description | Type | Reqd | Deflt
+ sets | Names of input cellSets | wordList | choice | -
+ set | Name of input cellSet | word | choice | -
\endverbatim
Note
- The order of precedence among the conditional mandatory entries from the
+ - The order of precedence among the conditional mandatory entries from the
highest to the lowest is \c sets, and \c set.
+ - The \c outside option applies to the cellSets individually.
See also
- Foam::topoSetSource
@@ -109,8 +111,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
-#ifndef cellToFace_H
-#define cellToFace_H
+#ifndef Foam_cellToFace_H
+#define Foam_cellToFace_H
#include "topoSetFaceSource.H"
#include "Enum.H"
@@ -121,7 +123,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
- Class cellToFace Declaration
+ Class cellToFace Declaration
\*---------------------------------------------------------------------------*/
class cellToFace
@@ -133,7 +135,8 @@ public:
enum cellAction
{
ALL,
- BOTH
+ BOTH,
+ OUTSIDE
};
diff --git a/src/meshTools/topoSet/faceZoneSources/cellToFaceZone/cellToFaceZone.C b/src/meshTools/topoSet/faceZoneSources/cellToFaceZone/cellToFaceZone.C
new file mode 100644
index 0000000000000000000000000000000000000000..ed6a4c3041d34307478fc5f68da26cf4726d8b94
--- /dev/null
+++ b/src/meshTools/topoSet/faceZoneSources/cellToFaceZone/cellToFaceZone.C
@@ -0,0 +1,287 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "cellToFaceZone.H"
+#include "polyMesh.H"
+#include "faceZoneSet.H"
+#include "cellSet.H"
+#include "syncTools.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(cellToFaceZone, 0);
+ addToRunTimeSelectionTable(topoSetSource, cellToFaceZone, word);
+ addToRunTimeSelectionTable(topoSetSource, cellToFaceZone, istream);
+
+ addToRunTimeSelectionTable(topoSetFaceZoneSource, cellToFaceZone, word);
+ addToRunTimeSelectionTable(topoSetFaceZoneSource, cellToFaceZone, istream);
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::cellToFaceZone::usage_
+(
+ cellToFaceZone::typeName,
+ "\n Usage: cellToFaceZone <slaveCellSet>\n\n"
+ " Select all outside faces in the cellSet."
+ " Orientated so slave side is in cellSet.\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::cellToFaceZone::selectFaces
+(
+ const cellSet& cSet,
+ bitSet& selectedFace,
+ bitSet& doFlip
+) const
+{
+ selectedFace.setSize(mesh_.nFaces());
+ selectedFace = false;
+
+ doFlip.setSize(mesh_.nFaces());
+ doFlip = false;
+
+
+ // Add all faces whose both neighbours are in set.
+
+ const label nInt = mesh_.nInternalFaces();
+ const labelList& own = mesh_.faceOwner();
+ const labelList& nei = mesh_.faceNeighbour();
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+
+ // Check all internal faces
+ for (label facei = 0; facei < nInt; ++facei)
+ {
+ const bool ownFound = cSet.found(own[facei]);
+ const bool neiFound = cSet.found(nei[facei]);
+
+ if (ownFound && !neiFound)
+ {
+ selectedFace.set(facei);
+ doFlip.set(facei, flip_);
+ }
+ else if (!ownFound && neiFound)
+ {
+ selectedFace.set(facei);
+ doFlip.set(facei, !flip_);
+ }
+ }
+
+ // Get coupled cell status
+ boolList neiInSet(mesh_.nBoundaryFaces(), false);
+
+ for (const polyPatch& pp : patches)
+ {
+ if (pp.coupled())
+ {
+ label facei = pp.start();
+ forAll(pp, i)
+ {
+ neiInSet[facei-nInt] = cSet.found(own[facei]);
+ ++facei;
+ }
+ }
+ }
+ syncTools::swapBoundaryFaceList(mesh_, neiInSet);
+
+
+ // Check all boundary faces
+ for (const polyPatch& pp : patches)
+ {
+ label facei = pp.start();
+ forAll(pp, i)
+ {
+ const bool ownFound = cSet.found(own[facei]);
+ const bool neiFound = neiInSet[facei-nInt];
+
+ if (ownFound && !neiFound)
+ {
+ selectedFace.set(facei);
+ doFlip.set(facei, flip_);
+ }
+ else if (!ownFound && neiFound)
+ {
+ selectedFace.set(facei);
+ doFlip.set(facei, !flip_);
+ }
+ ++facei;
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::cellToFaceZone::cellToFaceZone
+(
+ const polyMesh& mesh,
+ const word& setName,
+ const bool flip
+)
+:
+ topoSetFaceZoneSource(mesh),
+ names_(one{}, setName),
+ flip_(flip)
+{}
+
+
+Foam::cellToFaceZone::cellToFaceZone
+(
+ const polyMesh& mesh,
+ const dictionary& dict
+)
+:
+ topoSetFaceZoneSource(mesh),
+ names_(),
+ flip_(dict.getOrDefault("flip", false))
+{
+ // Look for 'sets' or 'set'
+ if (!dict.readIfPresent("sets", names_))
+ {
+ names_.resize(1);
+ dict.readEntry("set", names_.first());
+ }
+}
+
+
+Foam::cellToFaceZone::cellToFaceZone
+(
+ const polyMesh& mesh,
+ Istream& is
+)
+:
+ topoSetFaceZoneSource(mesh),
+ names_(one{}, word(checkIs(is))),
+ flip_(false)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::cellToFaceZone::applyToSet
+(
+ const topoSetSource::setAction action,
+ topoSet& set
+) const
+{
+ if (!isA<faceZoneSet>(set))
+ {
+ WarningInFunction
+ << "Operation only allowed on a faceZoneSet." << endl;
+ return;
+ }
+ else
+ {
+ faceZoneSet& zoneSet = refCast<faceZoneSet>(set);
+
+ if (action == topoSetSource::ADD || action == topoSetSource::NEW)
+ {
+ if (verbose_)
+ {
+ if (flip_)
+ {
+ Info<< " Adding all faces on outside of cellSet "
+ << flatOutput(names_)
+ << "; orientation pointing into cellSet" << endl;
+ }
+ else
+ {
+ Info<< " Adding all faces on outside of cellSet "
+ << flatOutput(names_)
+ << "; orientation pointing away from cellSet" << endl;
+ }
+ }
+
+ bitSet selectedFace(mesh_.nFaces());
+ bitSet doFlip(mesh_.nFaces());
+ for (const word& setName : names_)
+ {
+ // Load the sets
+ cellSet cSet(mesh_, setName);
+ // Select outside faces
+ selectFaces(cSet, selectedFace, doFlip);
+ }
+
+ // Start off from copy
+ DynamicList<label> newAddressing(zoneSet.addressing());
+ DynamicList<bool> newFlipMap(zoneSet.flipMap());
+
+ for (const label facei : selectedFace)
+ {
+ if (!zoneSet.found(facei))
+ {
+ newAddressing.append(facei);
+ newFlipMap.append(doFlip[facei]);
+ }
+ }
+
+ zoneSet.addressing().transfer(newAddressing);
+ zoneSet.flipMap().transfer(newFlipMap);
+ zoneSet.updateSet();
+ }
+ else if (action == topoSetSource::SUBTRACT)
+ {
+ if (verbose_)
+ {
+ Info<< " Removing all faces on outside of cellSet "
+ << flatOutput(names_)
+ << " ..." << endl;
+ }
+
+ bitSet selectedFace(mesh_.nFaces());
+ bitSet doFlip(mesh_.nFaces());
+ for (const word& setName : names_)
+ {
+ // Load the sets
+ cellSet cSet(mesh_, setName);
+ // Select outside faces
+ selectFaces(cSet, selectedFace, doFlip);
+ }
+
+ // Start off empty
+ DynamicList<label> newAddressing(zoneSet.addressing().size());
+ DynamicList<bool> newFlipMap(zoneSet.flipMap().size());
+
+ for (const label facei : selectedFace)
+ {
+ newAddressing.append(facei);
+ newFlipMap.append(doFlip[facei]);
+ }
+ zoneSet.addressing().transfer(newAddressing);
+ zoneSet.flipMap().transfer(newFlipMap);
+ zoneSet.updateSet();
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/topoSet/faceZoneSources/cellToFaceZone/cellToFaceZone.H b/src/meshTools/topoSet/faceZoneSources/cellToFaceZone/cellToFaceZone.H
new file mode 100644
index 0000000000000000000000000000000000000000..338d8f19e320c788dfbfb265e2a55f3dcba60ac0
--- /dev/null
+++ b/src/meshTools/topoSet/faceZoneSources/cellToFaceZone/cellToFaceZone.H
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::cellToFaceZone
+
+Description
+ A \c topoSetSource to select faces with only one
+ neighbour (i.e. outside) in a specified \c cellSet.
+
+ This is just a shortcut for
+ - extracting outside faces as a \c faceSet (\c cellToFace with \c outside).
+ - using \c setsToFaceZone to convert \c faceSet
+ and \c cellSet to oriented \c faceZone.
+
+ Operands:
+ \table
+ Operand | Type | Location
+ input | cellSet(s) | constant/polyMesh/sets/\<set\>
+ output | faceZone | constant/polyMesh/faceZones
+ \endtable
+
+Usage
+ Minimal example by using \c system/topoSetDict.actions:
+ \verbatim
+ {
+ // Mandatory entries
+ name <name>;
+ type faceZoneSet;
+ action <action>;
+ source cellToFaceZone;
+
+ // Select either of the below
+
+ // Option-1
+ sets
+ (
+ <word>
+ <word>
+ ...
+ );
+
+ // Option-2
+ set <word>;
+
+ // Optional entries
+ flip <bool>;
+ }
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Reqd | Deflt
+ name | Name of faceZone | word | yes | -
+ type | Type name: faceZoneSet | word | yes | -
+ action | Action applied on faces - see below | word | yes | -
+ source | Source name: cellToFaceZone | word | yes | -
+ set(s) | Name of input cellSet(s) containing the slave cells <!--
+ --> | word | yes | -
+ flip | Flag to select master/slave cells | bool | no | false
+ \endtable
+
+ Options for the \c action entry:
+ \verbatim
+ new | Create a new faceZone from selected faces
+ add | Add selected faces of a faceZoneSet into this faceZone
+ subtract | Remove selected faces of a faceZoneSet from this faceZone
+ \endverbatim
+
+Notes
+ - \c flip=true sets the orientation of faces
+ pointing into the \c cellSet, and vice versa.
+
+SourceFiles
+ cellToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_cellToFaceZone_H
+#define Foam_cellToFaceZone_H
+
+#include "topoSetFaceZoneSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class cellSet;
+
+/*---------------------------------------------------------------------------*\
+ Class cellToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellToFaceZone
+:
+ public topoSetFaceZoneSource
+{
+ // Private Data
+
+ //- Add usage string
+ static addToUsageTable usage_;
+
+ //- Names of cellSets to use
+ wordList names_;
+
+ //- Whether cellSet is slave cells or master cells
+ const bool flip_;
+
+
+ // Private Member Functions
+
+ //- Select outside faces of cellSet
+ void selectFaces
+ (
+ const cellSet& cSet,
+ bitSet& selectedFace,
+ bitSet& doFlip
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("cellToFaceZone");
+
+
+ // Constructors
+
+ //- Construct from components
+ cellToFaceZone
+ (
+ const polyMesh& mesh,
+ const word& cellSetName,
+ const bool flip
+ );
+
+ //- Construct from dictionary
+ cellToFaceZone(const polyMesh& mesh, const dictionary& dict);
+
+ //- Construct from Istream
+ cellToFaceZone(const polyMesh& mesh, Istream& is);
+
+
+ //- Destructor
+ virtual ~cellToFaceZone() = default;
+
+
+ // Member Functions
+
+ virtual void applyToSet
+ (
+ const topoSetSource::setAction action,
+ topoSet& set
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/Make/files b/src/thermophysicalModels/radiation/Make/files
index acd357afac61656868e98476b0ba57d3d3b540e2..08c27f94bd543ef5fe8c62d9e86df7bab900e069 100644
--- a/src/thermophysicalModels/radiation/Make/files
+++ b/src/thermophysicalModels/radiation/Make/files
@@ -12,6 +12,7 @@ radiationModels/opaqueSolid/opaqueSolid.C
radiationModels/solarLoad/solarLoad.C
radiationModels/solarLoad/faceShading/faceShading.C
radiationModels/solarLoad/faceReflecting/faceReflecting.C
+radiationModels/solarLoad/solarLoadBase.C
/* Scatter model */
submodels/scatterModel/scatterModel/scatterModel.C
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
index ee2cf3a33074aa1c212b58fdbe448b7d76aeb334..fb670bc53eefab0c01b84330578aef45a03dd8f5 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
@@ -31,6 +31,7 @@ License
#include "boundaryRadiationProperties.H"
#include "cyclicAMIPolyPatch.H"
#include "volFields.H"
+#include "surfaceFields.H"
#include "distributedTriSurfaceMesh.H"
#include "OBJstream.H"
@@ -43,153 +44,65 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-void Foam::faceShading::writeRays
-(
- const fileName& fName,
- const DynamicField<point>& endCf,
- const pointField& myFc
-)
-{
- OBJstream os(fName);
-
- Pout<< "Dumping rays to " << os.name() << endl;
-
- forAll(myFc, facei)
- {
- os.writeLine(myFc[facei], endCf[facei]);
- }
-}
-
-
-Foam::triSurface Foam::faceShading::triangulate
-(
- const labelHashSet& includePatches,
- const List<labelHashSet>& includeAllFacesPerPatch
-)
+void Foam::faceShading::calculate()
{
- const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
-
- // Storage for surfaceMesh. Size estimate.
- DynamicList<labelledTri> triangles(mesh_.nBoundaryFaces());
-
- label newPatchI = 0;
-
- for (const label patchI : includePatches)
- {
- const polyPatch& patch = bMesh[patchI];
- const pointField& points = patch.points();
-
- label nTriTotal = 0;
-
- if (includeAllFacesPerPatch[patchI].size())
- {
- for (const label patchFaceI : includeAllFacesPerPatch[patchI])
- {
- const face& f = patch[patchFaceI];
-
- faceList triFaces(f.nTriangles(points));
-
- label nTri = 0;
+ const auto& pbm = mesh_.boundaryMesh();
- f.triangles(points, nTri, triFaces);
-
- for (const face& f : triFaces)
- {
- triangles.append
- (
- labelledTri(f[0], f[1], f[2], newPatchI)
- );
- nTriTotal++;
- }
- }
- newPatchI++;
- }
- }
-
- triangles.shrink();
-
- // Create globally numbered tri surface
- triSurface rawSurface(triangles, mesh_.points());
-
- // Create locally numbered tri surface
- triSurface surface
+ const bitSet isOpaqueFace
(
- rawSurface.localFaces(),
- rawSurface.localPoints()
+ selectOpaqueFaces
+ (
+ radiation::boundaryRadiationProperties::New(mesh_),
+ patchIDs_,
+ zoneIDs_
+ )
);
- // Add patch names to surface
- surface.patches().setSize(newPatchI);
-
- newPatchI = 0;
-
- for (const label patchI : includePatches)
- {
- const polyPatch& patch = bMesh[patchI];
-
- if (includeAllFacesPerPatch[patchI].size())
- {
- surface.patches()[newPatchI].name() = patch.name();
- surface.patches()[newPatchI].geometricType() = patch.type();
-
- newPatchI++;
- }
- }
-
- return surface;
-}
-
-
-void Foam::faceShading::calculate()
-{
- const radiation::boundaryRadiationProperties& boundaryRadiation =
- radiation::boundaryRadiationProperties::New(mesh_);
-
- label nFaces = 0; //total number of direct hit faces
-
- const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
- DynamicList<point> dynCf(nFaces);
- DynamicList<label> dynFacesI;
-
- forAll(patches, patchI)
- {
- const polyPatch& pp = patches[patchI];
- const vectorField::subField cf = pp.faceCentres();
-
- if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
- {
- const tmp<scalarField> tt =
- boundaryRadiation.transmissivity(patchI);
- const scalarField& t = tt();
- const vectorField& n = pp.faceNormals();
+ // Find faces potentially hit by solar rays
+ // - correct normal
+ // - transmissivity 0
+ labelList hitFacesIds;
+ bitSet hitFacesFlips;
+ selectFaces
+ (
+ true, // use normal to do first filtering
+ isOpaqueFace,
+ patchIDs_,
+ zoneIDs_,
- forAll(n, faceI)
- {
- const vector nf(n[faceI]);
- if (((direction_ & nf) > 0) && (t[faceI] == 0.0))
- {
- dynFacesI.append(faceI + pp.start());
- dynCf.append(cf[faceI]);
- nFaces++;
- }
- }
- }
- }
+ hitFacesIds,
+ hitFacesFlips
+ );
Info<< "Number of 'potential' direct hits : "
- << returnReduce(nFaces, sumOp<label>()) << endl;
-
- labelList hitFacesIds(nFaces);
- hitFacesIds.transfer(dynFacesI);
+ << returnReduce(hitFacesIds.size(), sumOp<label>()) << endl;
- pointField Cfs(hitFacesIds.size());
- Cfs.transfer(dynCf);
// * * * * * * * * * * * * * * *
// Create distributedTriSurfaceMesh
Random rndGen(653213);
+ // Find potential obstructions. Include all faces that might potentially
+ // block (so ignore normal)
+ labelList blockingFacesIds;
+ bitSet blockingFacesFlips;
+ selectFaces
+ (
+ false, // use normal to do first filtering
+ isOpaqueFace,
+ patchIDs_,
+ zoneIDs_,
+
+ blockingFacesIds,
+ blockingFacesFlips
+ );
+
+ const triSurface localSurface = triangulate
+ (
+ blockingFacesIds,
+ blockingFacesFlips
+ );
+
// Determine mesh bounding boxes:
List<treeBoundBox> meshBb
(
@@ -210,36 +123,6 @@ void Foam::faceShading::calculate()
);
dict.add("mergeDistance", SMALL);
- labelHashSet includePatches;
- List<labelHashSet> includeAllFacesPerPatch(patches.size());
-
- forAll(patches, patchI)
- {
- const polyPatch& pp = patches[patchI];
- if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
- {
- includePatches.insert(patchI);
-
- const tmp<scalarField> tt =
- boundaryRadiation.transmissivity(patchI);
- const scalarField& tau = tt();
-
- forAll(pp, faceI)
- {
- if (tau[faceI] == 0.0)
- {
- includeAllFacesPerPatch[patchI].insert(faceI);
- }
- }
- }
- }
-
- triSurface localSurface = triangulate
- (
- includePatches,
- includeAllFacesPerPatch
- );
-
distributedTriSurfaceMesh surfacesMesh
(
IOobject
@@ -264,36 +147,28 @@ void Foam::faceShading::calculate()
returnReduce(5.0*mesh_.bounds().mag(), maxOp<scalar>());
// Calculate index of faces which have a direct hit (local)
- DynamicList<label> rayStartFace(nFaces + 0.01*nFaces);
// Shoot Rays
// * * * * * * * * * * * * * * * *
{
- DynamicField<point> start(nFaces);
+ DynamicField<point> start(hitFacesIds.size());
DynamicField<point> end(start.size());
DynamicList<label> startIndex(start.size());
- label i = 0;
- do
- {
- for (; i < Cfs.size(); i++)
- {
- const point& fc = Cfs[i];
-
- const label myFaceId = hitFacesIds[i];
+ const pointField& faceCentres = mesh_.faceCentres();
- const vector d(direction_*maxBounding);
+ const vector d(direction_*maxBounding);
- start.append(fc - 0.001*d);
-
- startIndex.append(myFaceId);
-
- end.append(fc - d);
-
- }
+ forAll(hitFacesIds, i)
+ {
+ const label facei = hitFacesIds[i];
+ const point& fc = faceCentres[facei];
- } while (returnReduceOr(i < Cfs.size()));
+ start.append(fc - 0.001*d);
+ startIndex.append(facei);
+ end.append(fc - d);
+ }
List<pointIndexHit> hitInfo(startIndex.size());
surfacesMesh.findLine(start, end, hitInfo);
@@ -301,11 +176,23 @@ void Foam::faceShading::calculate()
// Collect the rays which has 'only one not wall' obstacle between
// start and end.
// If the ray hit itself get stored in dRayIs
+
+ label nVisible = 0;
+ forAll(hitInfo, rayI)
+ {
+ if (!hitInfo[rayI].hit())
+ {
+ nVisible++;
+ }
+ }
+ rayStartFaces_.setSize(nVisible);
+ nVisible = 0;
+
forAll(hitInfo, rayI)
{
if (!hitInfo[rayI].hit())
{
- rayStartFace.append(startIndex[rayI]);
+ rayStartFaces_[nVisible++] = startIndex[rayI];
}
}
@@ -325,36 +212,39 @@ void Foam::faceShading::calculate()
end.clear();
}
- rayStartFaces_.transfer(rayStartFace);
-
if (debug)
{
- tmp<volScalarField> thitFaces
+ auto thitFaces = tmp<surfaceScalarField>::New
(
- new volScalarField
+ IOobject
(
- IOobject
- (
- "hitFaces",
- mesh_.time().timeName(),
- mesh_,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
+ "hitFaces",
+ mesh_.time().timeName(),
mesh_,
- dimensionedScalar(dimless, Zero)
- )
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh_,
+ dimensionedScalar(dimless, Zero)
);
- volScalarField& hitFaces = thitFaces.ref();
- volScalarField::Boundary& hitFacesBf = hitFaces.boundaryFieldRef();
+ surfaceScalarField& hitFaces = thitFaces.ref();
+ surfaceScalarField::Boundary& hitFacesBf = hitFaces.boundaryFieldRef();
hitFacesBf = 0.0;
- for (const label faceI : rayStartFaces_)
+ for (const label facei : rayStartFaces_)
{
- label patchID = patches.whichPatch(faceI);
- const polyPatch& pp = patches[patchID];
- hitFacesBf[patchID][faceI - pp.start()] = 1.0;
+ const label patchID = pbm.whichPatch(facei);
+ if (patchID == -1)
+ {
+ hitFaces[facei] = 1.0;
+ }
+ else
+ {
+ const polyPatch& pp = pbm[patchID];
+ hitFacesBf[patchID][facei - pp.start()] = 1.0;
+ }
}
hitFaces.write();
}
@@ -364,28 +254,290 @@ void Foam::faceShading::calculate()
}
+Foam::triSurface Foam::faceShading::triangulate
+(
+ const labelUList& faceIDs,
+ const bitSet& flipMap
+) const
+{
+ if (faceIDs.size() != flipMap.size())
+ {
+ FatalErrorInFunction << "Size problem :"
+ << "faceIDs:" << faceIDs.size()
+ << "flipMap:" << flipMap.size()
+ << exit(FatalError);
+ }
+
+ const auto& points = mesh_.points();
+ const auto& faces = mesh_.faces();
+ const auto& bMesh = mesh_.boundaryMesh();
+ const auto& fzs = mesh_.faceZones();
+
+ // Patching of surface:
+ // - non-processor patches
+ // - faceZones
+ // Note: check for faceZones on boundary? Who has priority?
+ geometricSurfacePatchList surfPatches(bMesh.nNonProcessor()+fzs.size());
+ labelList patchID(mesh_.nFaces(), -1);
+ {
+ label newPatchi = 0;
+ for (label patchi = 0; patchi < bMesh.nNonProcessor(); ++patchi)
+ {
+ const auto& pp = bMesh[patchi];
+
+ surfPatches[newPatchi] = geometricSurfacePatch
+ (
+ pp.name(),
+ newPatchi,
+ pp.type()
+ );
+ SubList<label>
+ (
+ patchID,
+ pp.size(),
+ pp.start()
+ ) = newPatchi;
+
+ newPatchi++;
+ }
+ for (const auto& fz : fzs)
+ {
+ surfPatches[newPatchi] = geometricSurfacePatch
+ (
+ fz.name(),
+ newPatchi,
+ fz.type()
+ );
+ UIndirectList<label>(patchID, fz) = newPatchi;
+
+ newPatchi++;
+ }
+ }
+
+
+ // Storage for surfaceMesh. Size estimate.
+ DynamicList<labelledTri> triangles(2*faceIDs.size());
+
+ // Work array
+ faceList triFaces;
+
+ forAll(faceIDs, i)
+ {
+ const label facei = faceIDs[i];
+ const bool flip = flipMap[i];
+ const label patchi = patchID[facei];
+ const face& f = faces[facei];
+
+ // Triangulate face
+ triFaces.setSize(f.nTriangles(points));
+ label nTri = 0;
+ f.triangles(points, nTri, triFaces);
+
+ for (const face& f : triFaces)
+ {
+ if (!flip)
+ {
+ triangles.append(labelledTri(f[0], f[1], f[2], patchi));
+ }
+ else
+ {
+ triangles.append(labelledTri(f[0], f[2], f[1], patchi));
+ }
+ }
+ }
+
+ triangles.shrink();
+
+ // Create globally numbered tri surface
+ triSurface rawSurface(triangles, mesh_.points());
+
+ // Create locally numbered tri surface
+ triSurface surface
+ (
+ rawSurface.localFaces(),
+ rawSurface.localPoints()
+ );
+
+ // Add patch names to surface
+ surface.patches().transfer(surfPatches);
+
+ return surface;
+}
+
+
+Foam::bitSet Foam::faceShading::selectOpaqueFaces
+(
+ const radiation::boundaryRadiationProperties& boundaryRadiation,
+ const labelUList& patchIDs,
+ const labelUList& zoneIDs
+) const
+{
+ const auto& pbm = mesh_.boundaryMesh();
+
+ bitSet isOpaqueFace(mesh_.nFaces(), false);
+
+ // Check selected patches
+ for (const label patchi : patchIDs)
+ {
+ const auto& pp = pbm[patchi];
+ tmp<scalarField> tt = boundaryRadiation.transmissivity(patchi);
+ const scalarField& t = tt.cref();
+
+ forAll(t, i)
+ {
+ isOpaqueFace[i + pp.start()] = (t[i] == 0.0);
+ }
+ }
+
+ // Check selected faceZones
+ const auto& fzs = mesh_.faceZones();
+
+ for (const label zonei : zoneIDs)
+ {
+ const auto& fz = fzs[zonei];
+
+ //- Note: slice mesh face centres preferentially
+ tmp<scalarField> tt = boundaryRadiation.zoneTransmissivity
+ (
+ zonei,
+ fz
+ );
+ const scalarField& t = tt.cref();
+
+ forAll(t, i)
+ {
+ isOpaqueFace[fz[i]] = (t[i] == 0.0);
+ }
+ }
+
+ return isOpaqueFace;
+}
+
+
+void Foam::faceShading::selectFaces
+(
+ const bool useNormal,
+ const bitSet& isCandidateFace,
+ const labelUList& patchIDs,
+ const labelUList& zoneIDs,
+
+ labelList& faceIDs,
+ bitSet& flipMap
+) const
+{
+ const auto& pbm = mesh_.boundaryMesh();
+
+ bitSet isSelected(mesh_.nFaces());
+ DynamicList<label> dynFaces(mesh_.nBoundaryFaces());
+ bitSet isFaceFlipped(mesh_.nFaces());
+
+ // Add patches
+ for (const label patchi : patchIDs)
+ {
+ const auto& pp = pbm[patchi];
+ const vectorField& n = pp.faceNormals();
+
+ forAll(n, i)
+ {
+ const label meshFacei = i + pp.start();
+ if
+ (
+ isCandidateFace[meshFacei]
+ && (
+ !useNormal
+ || ((direction_ & n[i]) > 0)
+ )
+ )
+ {
+ isSelected.set(meshFacei);
+ isFaceFlipped[meshFacei] = false;
+ dynFaces.append(meshFacei);
+ }
+ }
+ }
+
+
+ // Add faceZones
+ const auto& fzs = mesh_.faceZones();
+
+ for (const label zonei : zoneIDs)
+ {
+ const auto& fz = fzs[zonei];
+ const primitiveFacePatch& pp = fz();
+ const vectorField& n = pp.faceNormals();
+
+ forAll(n, i)
+ {
+ const label meshFacei = fz[i];
+
+ if
+ (
+ !isSelected[meshFacei]
+ && isCandidateFace[meshFacei]
+ && (
+ !useNormal
+ || ((direction_ & n[i]) > 0)
+ )
+ )
+ {
+ isSelected.set(meshFacei);
+ dynFaces.append(meshFacei);
+ isFaceFlipped[meshFacei] = fz.flipMap()[i];
+ }
+ }
+ }
+ faceIDs = std::move(dynFaces);
+ flipMap = bitSet(isFaceFlipped, faceIDs);
+}
+
+
+void Foam::faceShading::writeRays
+(
+ const fileName& fName,
+ const DynamicField<point>& endCf,
+ const pointField& myFc
+)
+{
+ OBJstream os(fName);
+
+ Pout<< "Dumping rays to " << os.name() << endl;
+
+ forAll(myFc, facei)
+ {
+ os.writeLine(myFc[facei], endCf[facei]);
+ }
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::faceShading::faceShading
(
const fvMesh& mesh,
- const vector dir,
- const labelList& hitFaceList
+ const vector& dir
)
:
mesh_(mesh),
+ patchIDs_(nonCoupledPatches(mesh)),
+ zoneIDs_(0),
direction_(dir),
- rayStartFaces_(hitFaceList)
-{}
+ rayStartFaces_(0)
+{
+ calculate();
+}
Foam::faceShading::faceShading
(
const fvMesh& mesh,
- const vector dir
+ const labelList& patchIDs,
+ const labelList& zoneIDs,
+ const vector& dir
)
:
mesh_(mesh),
+ patchIDs_(patchIDs),
+ zoneIDs_(zoneIDs),
direction_(dir),
rayStartFaces_(0)
{
@@ -395,6 +547,23 @@ Foam::faceShading::faceShading
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+Foam::labelList Foam::faceShading::nonCoupledPatches(const polyMesh& mesh)
+{
+ const auto& pbm = mesh.boundaryMesh();
+
+ DynamicList<label> ncPatches;
+ forAll(pbm, patchi)
+ {
+ const polyPatch& pp = pbm[patchi];
+ if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
+ {
+ ncPatches.append(patchi);
+ }
+ }
+ return ncPatches;
+}
+
+
void Foam::faceShading::correct()
{
calculate();
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.H b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.H
index 1374ec362c12242cd4c5f58415be46fe6ed4326e..6ef7b1574ff8d0939852739fee7775bdd04d927e 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.H
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
- Copyright (C) 2017-2018 OpenCFD Ltd.
+ Copyright (C) 2017-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -28,6 +28,7 @@ Class
Foam::faceShading
Description
+ Helper class to calculate visible faces for global, sun-like illumination.
faceShading uses the transmissivity value in the boundaryRadiationProperties
in order to evaluate which faces are "hit" by the "direction" vector.
@@ -41,31 +42,36 @@ SourceFiles
#ifndef faceShading_H
#define faceShading_H
-#include "fvMesh.H"
-#include "Time.H"
-#include "meshTools.H"
+#include "boundaryRadiationProperties.H"
#include "DynamicField.H"
-#include "labelIOList.H"
-#include "wallPolyPatch.H"
-#include "triSurfaceTools.H"
-
+#include "triSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
+// Forward Declarations
+class fvMesh;
+class polyMesh;
+
/*---------------------------------------------------------------------------*\
Class faceShading Declaration
\*---------------------------------------------------------------------------*/
class faceShading
{
- // Private data
+ // Private Data
//- Reference to mesh
const fvMesh& mesh_;
+ //- Patches to check for visibility
+ const labelList patchIDs_;
+
+ //- FaceZones to check for visibility
+ const labelList zoneIDs_;
+
//- Direction
vector direction_;
@@ -73,18 +79,42 @@ class faceShading
labelList rayStartFaces_;
- // Private members
+ // Private Member Functions
//- Calculate ray start faces
void calculate();
-
- //- Construct a triSurface from patches and faces on global local index
+ //- Construct a triSurface from selected faces on patches or faceZones.
+ // Resulting surface has regioning with mesh patches first, faceZones
+ // second.
triSurface triangulate
(
- const labelHashSet& includePatches,
- const List<labelHashSet>& includeAllFaces
- );
+ const labelUList& faceIDs,
+ const bitSet& flipMap
+ ) const;
+
+ //- Find opaque faces by looking at the radiation properties.
+ bitSet selectOpaqueFaces
+ (
+ const radiation::boundaryRadiationProperties& boundaryRadiation,
+ const labelUList& patchIDs,
+ const labelUList& zoneIDs
+ ) const;
+
+ //- Pick up candidate faces that might be blocking visibility for
+ // other faces:
+ // - if transmissivity is 0
+ // - (optional) if pointing wrong way
+ void selectFaces
+ (
+ const bool useNormal,
+ const bitSet& isCandidateFace,
+ const labelUList& patchIDs,
+ const labelUList& zoneIDs,
+
+ labelList& faceIDs,
+ bitSet& flipMap
+ ) const;
//- Write rays
void writeRays
@@ -109,19 +139,25 @@ public:
// Constructors
- //- Construct from components
+ //- Helper: return all uncoupled patches
+ static labelList nonCoupledPatches(const polyMesh& mesh);
+
+ //- Construct from mesh and vector to 'sun'. All uncoupled patches
+ // are checked for visibility. faceZones are ignored.
faceShading
(
const fvMesh& mesh,
- const vector dir,
- const labelList& hitFaceList
+ const vector& dir
);
- //- Construct from mesh and vector
+ //- Construct from mesh and vector to 'sun' and selected patches
+ // and faceZones.
faceShading
(
const fvMesh& mesh,
- const vector
+ const labelList& patchIDs,
+ const labelList& zoneIDs,
+ const vector& dir
);
@@ -133,7 +169,7 @@ public:
// Access
- //- const access to direction
+ //- Const access to direction
const vector direction() const
{
return direction_;
@@ -145,7 +181,7 @@ public:
return direction_;
}
- //- Access to rayStartFaces
+ //- Const access to rayStartFaces
const labelList& rayStartFaces() const
{
return rayStartFaces_;
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
index aaf4011c9b2d0842bbc08d1442f8267c194e27be..1b897d7a61e15392dfdd59edc42f993021c3ae12 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
@@ -82,28 +82,28 @@ void Foam::radiation::solarLoad::updateReflectedRays
{
if (includePatches[patchID])
{
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
qrBf[patchID] +=
- reflectedFaces_->qreflective(bandI).boundaryField()[patchID];
+ reflectedFaces_->qreflective(bandi).boundaryField()[patchID];
}
}
else
{
const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
- const labelUList& cellIs = patches[patchID].faceCells();
+ const labelUList& cellis = patches[patchID].faceCells();
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
- forAll(cellIs, i)
+ forAll(cellis, i)
{
- const label cellI = cellIs[i];
+ const label celli = cellis[i];
- Ru_[cellI] +=
+ Ru_[celli] +=
(
- reflectedFaces_->qreflective(bandI).
+ reflectedFaces_->qreflective(bandi).
boundaryField()[patchID][i] * sf[i]
- )/V[cellI];
+ )/V[celli];
}
}
}
@@ -115,8 +115,30 @@ bool Foam::radiation::solarLoad::updateHitFaces()
{
if (!hitFaces_)
{
- hitFaces_.reset(new faceShading(mesh_, solarCalc_.direction()));
- return true;
+ const boundaryRadiationProperties& boundaryRadiation =
+ boundaryRadiationProperties::New(mesh_);
+
+ const labelList activeFaceZoneIds(boundaryRadiation.faceZoneIds());
+
+ if (!activeFaceZoneIds.size())
+ {
+ hitFaces_.reset(new faceShading(mesh_, solarCalc_.direction()));
+ }
+ else
+ {
+ hitFaces_.reset
+ (
+ new faceShading
+ (
+ mesh_,
+ faceShading::nonCoupledPatches(mesh_),
+ activeFaceZoneIds,
+ solarCalc_.direction()
+ )
+ );
+ }
+
+ return true;
}
else
{
@@ -163,10 +185,10 @@ void Foam::radiation::solarLoad::updateAbsorptivity
for (const label patchID : includePatches)
{
absorptivity_[patchID].setSize(nBands_);
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
- absorptivity_[patchID][bandI] =
- boundaryRadiation.absorptivity(patchID, bandI);
+ absorptivity_[patchID][bandi] =
+ boundaryRadiation.absorptivity(patchID, bandi);
}
}
}
@@ -185,29 +207,37 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
// Reset qr and qrPrimary
qrBf = 0.0;
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
volScalarField::Boundary& qprimaryBf =
- qprimaryRad_[bandI].boundaryFieldRef();
+ qprimaryRad_[bandi].boundaryFieldRef();
qprimaryBf = 0.0;
forAll(hitFacesId, i)
{
- const label faceI = hitFacesId[i];
- const label patchID = patches.whichPatch(faceI);
+ const label facei = hitFacesId[i];
+ const label patchID = patches.whichPatch(facei);
+
+ if (patchID == -1)
+ {
+ continue;
+ }
+
const polyPatch& pp = patches[patchID];
- const label localFaceI = faceI - pp.start();
- const vector qPrim =
- solarCalc_.directSolarRad()*solarCalc_.direction();
+ const label localFaceI = facei - pp.start();
+ const vector qPrim
+ (
+ solarCalc_.directSolarRad()*solarCalc_.direction()
+ );
const vectorField& n = pp.faceNormals();
{
qprimaryBf[patchID][localFaceI] +=
(qPrim & n[localFaceI])
- * spectralDistribution_[bandI]
- * absorptivity_[patchID][bandI]()[localFaceI];
+ * spectralDistribution_[bandi]
+ * absorptivity_[patchID][bandi]()[localFaceI];
}
if (includeMappedPatchBasePatches[patchID])
@@ -217,13 +247,13 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
else
{
const vectorField& sf = mesh_.Sf().boundaryField()[patchID];
- const label cellI = pp.faceCells()[localFaceI];
+ const label celli = pp.faceCells()[localFaceI];
- Ru_[cellI] +=
+ Ru_[celli] +=
(qPrim & sf[localFaceI])
- * spectralDistribution_[bandI]
- * absorptivity_[patchID][bandI]()[localFaceI]
- / V[cellI];
+ * spectralDistribution_[bandi]
+ * absorptivity_[patchID][bandi]()[localFaceI]
+ / V[celli];
}
}
}
@@ -251,71 +281,36 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
const vectorField n = pp.faceNormals();
- const labelList& cellIds = pp.faceCells();
+ const labelList& cellids = pp.faceCells();
- forAll(n, faceI)
- {
- const scalar cosEpsilon(verticalDir_ & -n[faceI]);
-
- scalar Ed(0.0);
- scalar Er(0.0);
- const scalar cosTheta(solarCalc_.direction() & -n[faceI]);
+ // Calculate diffusive radiance
+ // contribution from sky and ground
+ tmp<scalarField> tdiffuseSolarRad =
+ solarCalc_.diffuseSolarRad(n);
+ const scalarField& diffuseSolarRad = tdiffuseSolarRad.cref();
- {
- // Above the horizon
- if (cosEpsilon == 0.0)
- {
- // Vertical walls
- scalar Y(0);
-
- if (cosTheta > -0.2)
- {
- Y = 0.55+0.437*cosTheta + 0.313*sqr(cosTheta);
- }
- else
- {
- Y = 0.45;
- }
-
- Ed = solarCalc_.C()*Y*solarCalc_.directSolarRad();
- }
- else
- {
- //Other than vertical walls
- Ed =
- solarCalc_.C()
- * solarCalc_.directSolarRad()
- * (1.0 + cosEpsilon)/2.0;
- }
-
- // Ground reflected
- Er =
- solarCalc_.directSolarRad()
- * (solarCalc_.C() + Foam::sin(solarCalc_.beta()))
- * solarCalc_.groundReflectivity()
- * (1.0 - cosEpsilon)/2.0;
- }
-
- const label cellI = cellIds[faceI];
+ forAll(n, facei)
+ {
+ const label celli = cellids[facei];
if (includeMappedPatchBasePatches[patchID])
{
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
- qrBf[patchID][faceI] +=
- (Ed + Er)
- * spectralDistribution_[bandI]
- * absorptivity_[patchID][bandI]()[faceI];
+ qrBf[patchID][facei] +=
+ diffuseSolarRad[facei]
+ * spectralDistribution_[bandi]
+ * absorptivity_[patchID][bandi]()[facei];
}
}
else
{
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
- Ru_[cellI] +=
- (Ed + Er)
- * spectralDistribution_[bandI]
- * absorptivity_[patchID][bandI]()[faceI]
- * sf[faceI]/V[cellI];
+ Ru_[celli] +=
+ diffuseSolarRad[facei]
+ * spectralDistribution_[bandi]
+ * absorptivity_[patchID][bandi]()[facei]
+ * sf[facei]/V[celli];
}
}
}
@@ -331,30 +326,30 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
const polyPatch& pp = patches[patchID];
const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
- const labelList& cellIds = pp.faceCells();
- forAll(pp, faceI)
+ const labelList& cellids = pp.faceCells();
+ forAll(pp, facei)
{
- const label cellI = cellIds[faceI];
+ const label celli = cellids[facei];
if (includeMappedPatchBasePatches[patchID])
{
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
- qrBf[patchID][faceI] +=
+ qrBf[patchID][facei] +=
solarCalc_.diffuseSolarRad()
- * spectralDistribution_[bandI]
- * absorptivity_[patchID][bandI]()[faceI];
+ * spectralDistribution_[bandi]
+ * absorptivity_[patchID][bandi]()[facei];
}
}
else
{
- for (label bandI = 0; bandI < nBands_; ++bandI)
+ for (label bandi = 0; bandi < nBands_; ++bandi)
{
- Ru_[cellI] +=
+ Ru_[celli] +=
(
- spectralDistribution_[bandI]
- * absorptivity_[patchID][bandI]()[faceI]
+ spectralDistribution_[bandi]
+ * absorptivity_[patchID][bandi]()[facei]
* solarCalc_.diffuseSolarRad()
- )*sf[faceI]/V[cellI];
+ )*sf[facei]/V[celli];
}
}
}
@@ -387,16 +382,16 @@ void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
qprimaryRad_.setSize(nBands_);
- forAll(qprimaryRad_, bandI)
+ forAll(qprimaryRad_, bandi)
{
qprimaryRad_.set
(
- bandI,
+ bandi,
new volScalarField
(
IOobject
(
- "qprimaryRad_" + Foam::name(bandI) ,
+ "qprimaryRad_" + Foam::name(bandi) ,
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
@@ -408,21 +403,12 @@ void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
);
}
- if (coeffs.readIfPresent("gridUp", verticalDir_))
- {
- verticalDir_.normalise();
- }
- else
- {
- const uniformDimensionedVectorField& g =
- meshObjects::gravity::New(mesh_.time());
- verticalDir_ = (-g/mag(g)).value();
- }
-
coeffs.readIfPresent("solidCoupled", solidCoupled_);
coeffs.readIfPresent("wallCoupled", wallCoupled_);
coeffs.readIfPresent("updateAbsorptivity", updateAbsorptivity_);
coeffs.readEntry("useReflectedRays", useReflectedRays_);
+
+ (void) updateHitFaces();
}
/*
@@ -516,11 +502,11 @@ void Foam::radiation::solarLoad::calculateQdiff
// Weight emissivity by spectral distribution
scalarField e(pp.size(), 0.0);
- for (label bandI = 0; bandI < nBands_; bandI++)
+ for (label bandi = 0; bandi < nBands_; bandi++)
{
const tmp<scalarField> te =
- spectralDistribution_[bandI]
- *boundaryRadiation.diffReflectivity(patchID, bandI);
+ spectralDistribution_[bandi]
+ *boundaryRadiation.diffReflectivity(patchID, bandi);
e += te();
}
@@ -543,15 +529,15 @@ void Foam::radiation::solarLoad::calculateQdiff
scalar fullArea = 0.0;
forAll(fineFaces, j)
{
- label faceI = fineFaces[j];
- label globalFaceI = faceI + pp.start();
+ label facei = fineFaces[j];
+ label globalFaceI = facei + pp.start();
if (hitFacesId.found(globalFaceI))
{
- fullArea += sf[faceI];
+ fullArea += sf[facei];
}
- Eave[coarseI] += (e[faceI]*sf[faceI])/fineArea;
- Tave[coarseI] += (pow4(Tf[faceI])*sf[faceI])/fineArea;
+ Eave[coarseI] += (e[facei]*sf[facei])/fineArea;
+ Tave[coarseI] += (pow4(Tf[facei])*sf[facei])/fineArea;
}
localCoarsePartialArea[compactI++] = fullArea/fineArea;
}
@@ -616,11 +602,11 @@ void Foam::radiation::solarLoad::calculateQdiff
scalarField a(ppf.size(), Zero);
- for (label bandI = 0; bandI < nBands_; bandI++)
+ for (label bandi = 0; bandi < nBands_; bandi++)
{
const tmp<scalarField> ta =
- spectralDistribution_[bandI]
- * absorptivity_[patchID][bandI]();
+ spectralDistribution_[bandi]
+ * absorptivity_[patchID][bandi]();
a += ta();
}
@@ -635,8 +621,8 @@ void Foam::radiation::solarLoad::calculateQdiff
scalar aAve = 0.0;
forAll(fineFaces, j)
{
- label faceI = fineFaces[j];
- aAve += (a[faceI]*sf[faceI])/fineArea;
+ label facei = fineFaces[j];
+ aAve += (a[facei]*sf[facei])/fineArea;
}
const scalarList& vf = FmyProc[locaFaceI];
@@ -689,8 +675,8 @@ void Foam::radiation::solarLoad::calculateQdiff
const labelList& fineFaces = coarseToFine_[i][coarseFaceID];
forAll(fineFaces, k)
{
- label faceI = fineFaces[k];
- qrp[faceI] = localqDiffusive[compactId];
+ label facei = fineFaces[k];
+ qrp[facei] = localqDiffusive[compactId];
}
compactId ++;
}
@@ -711,12 +697,12 @@ void Foam::radiation::solarLoad::calculateQdiff
else
{
const polyPatch& pp = patches[patchID];
- const labelList& cellIds = pp.faceCells();
+ const labelList& cellids = pp.faceCells();
const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
- forAll(pp, faceI)
+ forAll(pp, facei)
{
- const label cellI = cellIds[faceI];
- Ru_[cellI] += qSecond[faceI]*sf[faceI]/V[cellI];
+ const label celli = cellids[facei];
+ Ru_[celli] += qSecond[facei]*sf[facei]/V[celli];
}
}
}
@@ -728,6 +714,7 @@ void Foam::radiation::solarLoad::calculateQdiff
Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
:
radiationModel(typeName, T),
+ solarLoadBase(mesh_),
solarCalc_(coeffs_, mesh_),
dict_(coeffs_),
qr_
@@ -762,7 +749,6 @@ Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
spectralDistribution_(),
spectralDistributions_(),
qprimaryRad_(0),
- verticalDir_(Zero),
nBands_(0),
updateTimeIndex_(0),
solidCoupled_(true),
@@ -782,6 +768,7 @@ Foam::radiation::solarLoad::solarLoad
)
:
radiationModel(typeName, dict, T),
+ solarLoadBase(mesh_),
solarCalc_(dict, mesh_),
dict_(dict),
qr_
@@ -816,7 +803,6 @@ Foam::radiation::solarLoad::solarLoad
spectralDistribution_(),
spectralDistributions_(),
qprimaryRad_(0),
- verticalDir_(Zero),
nBands_(0),
updateTimeIndex_(0),
solidCoupled_(true),
@@ -879,7 +865,7 @@ void Foam::radiation::solarLoad::calculate()
const bool timeDependentLoad =
solarCalc_.sunLoadModel() == solarCalculator::mSunLoadTimeDependent;
- if (facesChanged || timeDependentLoad)
+ if (firstIter_ || facesChanged || timeDependentLoad)
{
// Reset Ru
Ru_ = dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), Zero);
@@ -953,4 +939,18 @@ Foam::radiation::solarLoad::Ru() const
}
+const Foam::solarCalculator&
+Foam::radiation::solarLoad::solarCalculatorRef() const
+{
+ return solarCalc_;
+}
+
+
+const Foam::faceShading&
+Foam::radiation::solarLoad::faceShadingRef() const
+{
+ return hitFaces_();
+}
+
+
// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H
index f657fac0bda4ba1c66439336107d396d349b8f4a..2ccb31c3d8d4b101c3dc0d0eccb1afa1f74e6dce 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.H
@@ -99,6 +99,7 @@ SourceFiles
#define radiation_solarLoad_H
#include "radiationModel.H"
+#include "solarLoadBase.H"
#include "volFields.H"
#include "faceShading.H"
#include "faceReflecting.H"
@@ -118,7 +119,8 @@ namespace radiation
class solarLoad
:
- public radiationModel
+ public radiationModel,
+ public solarLoadBase
{
// Private Data
@@ -153,9 +155,6 @@ class solarLoad
//- Primary solar radiative heat flux per band [W/m2]
PtrList<volScalarField> qprimaryRad_;
- //- Vertical direction (default is g vector)
- vector verticalDir_;
-
//- Number of bands
label nBands_;
@@ -262,6 +261,12 @@ public:
{
return qprimaryRad_[bandI];
}
+
+ //- Return const reference to the solar calculator
+ virtual const solarCalculator& solarCalculatorRef() const;
+
+ //- Return const reference to the face shading calculator
+ virtual const faceShading& faceShadingRef() const;
};
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoadBase.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoadBase.C
new file mode 100644
index 0000000000000000000000000000000000000000..617e2778fbb510ef24f078c8cc0eba008dfa808f
--- /dev/null
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoadBase.C
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "solarLoadBase.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+ defineTypeNameAndDebug(solarLoadBase, 0);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::radiation::solarLoadBase::solarLoadBase
+(
+ const fvMesh& mesh
+)
+:
+ regIOobject
+ (
+ IOobject
+ (
+ solarLoadBase::typeName,
+ mesh.polyMesh::instance(),
+ mesh
+ )
+ )
+{}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoadBase.H b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoadBase.H
new file mode 100644
index 0000000000000000000000000000000000000000..5398c5a8b777f140ed2979212e0f8b9781bfb92f
--- /dev/null
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoadBase.H
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::radiation::solarLoadBase
+
+Group
+ grpRadiationModels
+
+Description
+ Base class for solarLoad models.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef radiation_solarLoadBase_H
+#define radiation_solarLoadBase_H
+
+#include "regIOobject.H"
+#include "fvMesh.H"
+#include "solarCalculator.H"
+#include "faceShading.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace radiation
+{
+
+/*---------------------------------------------------------------------------*\
+ Class solarLoadBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class solarLoadBase
+:
+ public regIOobject
+{
+public:
+
+ //- Runtime type information
+ TypeName("solarLoadBase");
+
+
+ // Constructors
+
+ //- Construct
+ solarLoadBase(const fvMesh& mesh);
+
+
+ //- Destructor
+ virtual ~solarLoadBase() = default;
+
+
+ // Member Functions
+
+ // Access
+
+ //- Return const reference to the solar calculator
+ virtual const solarCalculator& solarCalculatorRef() const = 0;
+
+ //- Return const reference to the face shading calculator
+ virtual const faceShading& faceShadingRef() const = 0;
+
+
+ // IO
+
+ bool writeData(Foam::Ostream&) const
+ {
+ return true;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace radiation
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.C b/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.C
index c409ec246701e7709c827f61815471e3cf5a5801..35129a36430a5a2f20793f9aba43a8bbc0ed4b71 100644
--- a/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.C
+++ b/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.C
@@ -52,7 +52,8 @@ Foam::radiation::boundaryRadiationProperties::boundaryRadiationProperties
Foam::GeometricMeshObject,
boundaryRadiationProperties
>(mesh),
- radBoundaryPropertiesPtrList_(mesh.boundary().size())
+ radBoundaryPropertiesPtrList_(mesh.boundary().size()),
+ radZonePropertiesPtrList_(mesh.faceZones().size())
{
IOobject boundaryIO
(
@@ -75,30 +76,110 @@ Foam::radiation::boundaryRadiationProperties::boundaryRadiationProperties
// Model number of bands
label nBands = radiation.nBands();
+ // Load in dictionary
const IOdictionary radiationDict(boundaryIO);
- forAll(mesh.boundary(), patchi)
- {
- const polyPatch& pp = mesh.boundaryMesh()[patchi];
- const dictionary* subDictPtr = radiationDict.findDict(pp.name());
+ wordHashSet matchedEntries;
- if (subDictPtr)
+ // Match patches
+ {
+ const auto& pbm = mesh.boundaryMesh();
+ for (const auto& pp : pbm)
{
- const dictionary& dict = *subDictPtr;
+ const label patchi = pp.index();
+
+ const auto* ePtr = radiationDict.findEntry(pp.name());
+
+ if (ePtr->isDict())
+ {
+ radBoundaryPropertiesPtrList_.set
+ (
+ patchi,
+ boundaryRadiationPropertiesPatch::New(ePtr->dict(), pp)
+ );
+
+ matchedEntries.insert(pp.name());
+
+ if
+ (
+ nBands
+ != radBoundaryPropertiesPtrList_[patchi].nBands()
+ )
+ {
+ FatalErrorInFunction
+ << "Radiation bands : " << nBands << nl
+ << "Bands on patch : " << patchi << " is "
+ << radBoundaryPropertiesPtrList_[patchi].nBands()
+ << abort(FatalError);
+ }
+ }
+ }
+ }
+
- radBoundaryPropertiesPtrList_[patchi].reset
- (
- boundaryRadiationPropertiesPatch::New(dict, pp)
- );
+ // Match faceZones if any dictionary entries have not been used for
+ // patch matching.
+ //
+ // Note: radiation properties are hardcoded to take patch reference.
+ // Supply patch0 for now.
+ {
+ const auto& dummyRef = mesh.boundaryMesh()[0];
+
+ const auto& fzs = mesh.faceZones();
- if (nBands != radBoundaryPropertiesPtrList_[patchi]->nBands())
+ for (const auto& fz : fzs)
+ {
+ const label zonei = fz.index();
+
+ if (!matchedEntries.found(fz.name()))
{
- FatalErrorInFunction
- << "Radiation bands : " << nBands << nl
- << "Bands on patch : " << patchi << " is "
- << radBoundaryPropertiesPtrList_[patchi]->nBands()
- << abort(FatalError);
+ // Note: avoid wildcard matches. Assume user explicitly
+ // provided information for faceZones.
+ const auto* ePtr = radiationDict.findEntry
+ (
+ fz.name(),
+ keyType::LITERAL
+ );
+
+ // Skip face zones that are not used in boundary radiation
+ if (!ePtr)
+ {
+ continue;
+ }
+
+ if (ePtr->isDict())
+ {
+ const dictionary& dict = ePtr->dict();
+
+ radZonePropertiesPtrList_.set
+ (
+ zonei,
+ boundaryRadiationPropertiesPatch::New
+ (
+ dict,
+ dummyRef
+ )
+ );
+
+ matchedEntries.insert(fz.name());
+
+ if
+ (
+ nBands
+ != radZonePropertiesPtrList_[zonei].nBands()
+ )
+ {
+ FatalErrorInFunction
+ << "Radiation bands : " << nBands << nl
+ << "Bands on zone : " << zonei << " is "
+ << radBoundaryPropertiesPtrList_
+ [
+ zonei
+ ].nBands()
+ << abort(FatalError);
+ }
+ }
}
}
}
@@ -117,9 +198,9 @@ Foam::radiation::boundaryRadiationProperties::emissivity
scalarField* T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->e
+ return radBoundaryPropertiesPtrList_[patchi].e
(
bandi,
incomingDirection,
@@ -140,17 +221,17 @@ Foam::radiation::boundaryRadiationProperties::emissivity
Foam::scalar Foam::radiation::boundaryRadiationProperties::faceEmissivity
(
const label patchi,
- const label faceI,
+ const label facei,
const label bandi,
vector incomingDirection,
scalar T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->e
+ return radBoundaryPropertiesPtrList_[patchi].e
(
- faceI,
+ facei,
bandi,
incomingDirection,
T
@@ -176,9 +257,9 @@ Foam::radiation::boundaryRadiationProperties::absorptivity
scalarField* T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->a
+ return radBoundaryPropertiesPtrList_[patchi].a
(
bandi,
incomingDirection,
@@ -199,17 +280,17 @@ Foam::radiation::boundaryRadiationProperties::absorptivity
Foam::scalar Foam::radiation::boundaryRadiationProperties::faceAbsorptivity
(
const label patchi,
- const label faceI,
+ const label facei,
const label bandi,
vector incomingDirection,
scalar T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->a
+ return radBoundaryPropertiesPtrList_[patchi].a
(
- faceI,
+ facei,
bandi,
incomingDirection,
T
@@ -235,9 +316,9 @@ Foam::radiation::boundaryRadiationProperties::transmissivity
scalarField* T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->t
+ return radBoundaryPropertiesPtrList_[patchi].t
(
bandi,
incomingDirection,
@@ -258,17 +339,17 @@ Foam::radiation::boundaryRadiationProperties::transmissivity
Foam::scalar Foam::radiation::boundaryRadiationProperties::faceTransmissivity
(
const label patchi,
- const label faceI,
+ const label facei,
const label bandi,
vector incomingDirection,
scalar T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->t
+ return radBoundaryPropertiesPtrList_[patchi].t
(
- faceI,
+ facei,
bandi,
incomingDirection,
T
@@ -285,6 +366,43 @@ Foam::scalar Foam::radiation::boundaryRadiationProperties::faceTransmissivity
}
+Foam::tmp<Foam::scalarField>
+Foam::radiation::boundaryRadiationProperties::zoneTransmissivity
+(
+ const label zonei,
+ const labelUList& faceIDs,
+ const label bandi,
+ vector incomingDirection,
+ scalar T
+) const
+{
+ if (radZonePropertiesPtrList_.set(zonei))
+ {
+ auto tfld = tmp<scalarField>::New(faceIDs.size());
+ auto& fld = tfld.ref();
+ forAll(fld, i)
+ {
+ fld[i] = radZonePropertiesPtrList_[zonei].t
+ (
+ faceIDs[i],
+ bandi,
+ incomingDirection,
+ T
+ );
+ }
+ return tfld;
+ }
+
+ FatalErrorInFunction
+ << "Zone : " << mesh().faceZones()[zonei].name()
+ << " is not found in the boundaryRadiationProperties. "
+ << "Please add it"
+ << exit(FatalError);
+
+ return tmp<scalarField>::New();
+}
+
+
Foam::tmp<Foam::scalarField>
Foam::radiation::boundaryRadiationProperties::diffReflectivity
(
@@ -294,9 +412,9 @@ Foam::radiation::boundaryRadiationProperties::diffReflectivity
scalarField* T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->rDiff
+ return radBoundaryPropertiesPtrList_[patchi].rDiff
(
bandi,
incomingDirection,
@@ -317,17 +435,17 @@ Foam::radiation::boundaryRadiationProperties::diffReflectivity
Foam::scalar Foam::radiation::boundaryRadiationProperties::faceDiffReflectivity
(
const label patchi,
- const label faceI,
+ const label facei,
const label bandi,
vector incomingDirection,
scalar T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->rDiff
+ return radBoundaryPropertiesPtrList_[patchi].rDiff
(
- faceI,
+ facei,
bandi,
incomingDirection,
T
@@ -353,9 +471,9 @@ Foam::radiation::boundaryRadiationProperties::specReflectivity
scalarField* T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->rSpec
+ return radBoundaryPropertiesPtrList_[patchi].rSpec
(
bandi,
incomingDirection,
@@ -376,17 +494,17 @@ Foam::radiation::boundaryRadiationProperties::specReflectivity
Foam::scalar Foam::radiation::boundaryRadiationProperties::faceSpecReflectivity
(
const label patchi,
- const label faceI,
+ const label facei,
const label bandi,
vector incomingDirection,
scalar T
) const
{
- if (radBoundaryPropertiesPtrList_[patchi])
+ if (radBoundaryPropertiesPtrList_.set(patchi))
{
- return radBoundaryPropertiesPtrList_[patchi]->rSpec
+ return radBoundaryPropertiesPtrList_[patchi].rSpec
(
- faceI,
+ facei,
bandi,
incomingDirection,
T
diff --git a/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.H b/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.H
index 0e65a40e341aa14ccab306b04f1db1fadb1688c3..69dd3973d5b4d2640222e149043adfc330a202a2 100644
--- a/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.H
+++ b/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationProperties.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2015-2018 OpenCFD Ltd.
+ Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -63,12 +63,16 @@ class boundaryRadiationProperties
boundaryRadiationProperties
>
{
- // Private data
+ // Private Data
- //- Ptr list of boundaryRadiationProperties
- List<autoPtr<boundaryRadiationPropertiesPatch>>
+ //- Per patch the boundaryRadiationProperties
+ PtrList<boundaryRadiationPropertiesPatch>
radBoundaryPropertiesPtrList_;
+ //- Per (face)zone the boundaryRadiationProperties
+ PtrList<boundaryRadiationPropertiesPatch>
+ radZonePropertiesPtrList_;
+
public:
@@ -84,6 +88,22 @@ public:
// Member Functions
+ //- Return identifiers of face zones activated for boundary radiation
+ const labelList faceZoneIds() const
+ {
+ DynamicList<label> ncZones;
+
+ forAll(radZonePropertiesPtrList_, i)
+ {
+ if (radZonePropertiesPtrList_.test(i))
+ {
+ ncZones.append(i);
+ }
+ }
+
+ return ncZones;
+ }
+
//- Access boundary emissivity on patch
tmp<scalarField> emissivity
(
@@ -141,6 +161,20 @@ public:
scalar T = 0
) const;
+ // Specialisations for faceZones
+
+ //- Access transmissivity on set of (internal) faces. Zone name only
+ // used to lookup the properties in boundaryRadiationProperties
+ tmp<scalarField> zoneTransmissivity
+ (
+ const label zoneI,
+ const labelUList& faceIDs, // internal faces
+ const label bandI = 0,
+ vector incomingDirection = Zero,
+ scalar T = 0
+ ) const;
+
+
//- Access boundary diffuse reflectivity on patch
tmp<scalarField> diffReflectivity
(
diff --git a/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationPropertiesPatch.H b/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationPropertiesPatch.H
index dc9a057da431c74e3fa363fe75eb0e9d496bfa37..afb1078a43a343705e223363ccc46808e2b28063 100644
--- a/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationPropertiesPatch.H
+++ b/src/thermophysicalModels/radiation/submodels/boundaryRadiationProperties/boundaryRadiationPropertiesPatch.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2015-2018 OpenCFD Ltd.
+ Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -56,9 +56,7 @@ namespace radiation
class boundaryRadiationPropertiesPatch
{
-private:
-
- // Private data
+ // Private Data
//- Dictionary
const dictionary dict_;
@@ -67,7 +65,7 @@ private:
const polyPatch& patch_;
- // Private functions
+ // Private Member Functions
//- Return nbr patch index
label nbrPatchIndex() const;
@@ -128,7 +126,7 @@ public:
virtual ~boundaryRadiationPropertiesPatch() = default;
- // Member functions
+ // Member Functions
//- Return absorptionEmissionModel
const wallAbsorptionEmissionModel& absorptionEmission() const;
diff --git a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
index e14c88bb813d92f442fc3b739166b30a3e24ad73..9cf0a65f6ba2bfaa81e393224db32facad4fe0de 100644
--- a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
+++ b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2015-2021 OpenCFD Ltd.
+ Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -333,4 +333,60 @@ void Foam::solarCalculator::correctDiffuseSolarRad()
}
+Foam::tmp<Foam::scalarField> Foam::solarCalculator::diffuseSolarRad
+(
+ const vectorField& n
+) const
+{
+ auto tload = tmp<scalarField>::New(n.size());
+ auto& load = tload.ref();
+
+ forAll(n, facei)
+ {
+ const scalar cosEpsilon(gridUp_ & -n[facei]);
+
+ scalar Ed = 0;
+ scalar Er = 0;
+ const scalar cosTheta(direction_ & -n[facei]);
+
+ // Above the horizon
+ if (cosEpsilon == 0.0)
+ {
+ // Vertical walls
+ scalar Y = 0;
+
+ if (cosTheta > -0.2)
+ {
+ Y = 0.55+0.437*cosTheta + 0.313*sqr(cosTheta);
+ }
+ else
+ {
+ Y = 0.45;
+ }
+
+ Ed = C_*Y*directSolarRad_;
+ }
+ else
+ {
+ //Other than vertical walls
+ Ed =
+ C_
+ * directSolarRad_
+ * 0.5*(1.0 + cosEpsilon);
+ }
+
+ // Ground reflected
+ Er =
+ directSolarRad_
+ * (C_ + Foam::sin(beta_))
+ * groundReflectivity_
+ * 0.5*(1.0 - cosEpsilon);
+
+ load[facei] = Ed + Er;
+ }
+
+ return tload;
+}
+
+
// ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.H b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.H
index 0ed978ea713d5f41eba30acad5188e994edf806a..e82528d7ef39aa002540a017e057d0b4e924340f 100644
--- a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.H
+++ b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2015-2021 OpenCFD Ltd.
+ Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -395,6 +395,11 @@ public:
//- Correct diffuse solar irradiation
void correctDiffuseSolarRad();
+ //- Return the diffusive solar irradiation
+ //- for sunLoadModel = fairWeather
+ //- or for sunLoadModel = theoreticalMaximum
+ tmp<scalarField> diffuseSolarRad(const vectorField& n) const;
+
// Access
diff --git a/tutorials/incompressible/pisoFoam/RAS/cavity/system/topoSetDict b/tutorials/incompressible/pisoFoam/RAS/cavity/system/topoSetDict
index 801eafd409ce9243ff4dfc1791cc97c2a3249959..e20111f1e33e7566f38a388505f4abef6551855d 100644
--- a/tutorials/incompressible/pisoFoam/RAS/cavity/system/topoSetDict
+++ b/tutorials/incompressible/pisoFoam/RAS/cavity/system/topoSetDict
@@ -303,6 +303,15 @@ actions
set cylinder1;
}
+ {
+ name faceCell1c;
+ type faceSet;
+ action new;
+ source cellToFace;
+ option outside;
+ set cylinder1;
+ }
+
{
name faceCylinder1;
type faceSet;
@@ -430,6 +439,15 @@ actions
// flip true;
// }
+ {
+ name cellToFaceZone1;
+ type faceZoneSet;
+ action new;
+ source cellToFaceZone;
+ set cylinder1;
+ flip true;
+ }
+
{
name pointBox1;
type pointSet;
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/plantCd b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/Cd
similarity index 97%
rename from tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/plantCd
rename to tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/Cd
index 937ded772dac0f8396444ae7a3a327b8a6ca6467..60ea92267704891cb1d0d65717478c807f54dfef 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/plantCd
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/Cd
@@ -10,7 +10,7 @@ FoamFile
version 2.0;
format ascii;
class volScalarField;
- object plantCd;
+ object Cd;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/leafAreaDensity b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/LAD
similarity index 97%
rename from tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/leafAreaDensity
rename to tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/LAD
index 3b7b597b54c7d6efa5f27106c55ce763837eb27a..c6f2009fae11faad7791490fb03fc82c6556fa3d 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/leafAreaDensity
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/common/0.orig/LAD
@@ -10,7 +10,7 @@ FoamFile
version 2.0;
format ascii;
class volScalarField;
- object leafAreaDensity;
+ object LAD;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/plantCd b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/Cd
similarity index 100%
rename from tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/plantCd
rename to tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/Cd
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/leafAreaDensity b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/LAD
similarity index 97%
rename from tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/leafAreaDensity
rename to tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/LAD
index 3a65f8347cc5a2b2ae3b8fd95292a1dde84e4b5a..48b1a763066910d85ce3d6440d79ff3010d61ccc 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/leafAreaDensity
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/common/0.orig/LAD
@@ -10,7 +10,7 @@ FoamFile
version 2.0;
format ascii;
class volScalarField;
- object leafAreaDensity;
+ object LAD;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/neutral/system/setFieldsDict b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/neutral/system/setFieldsDict
index 725d7d2cbd54ce1991a55b73aaa8a09285087de1..4752c8f1a5efadc71d1e6b729bfe30aa5dff6077 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/neutral/system/setFieldsDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/neutral/system/setFieldsDict
@@ -17,8 +17,8 @@ FoamFile
defaultFieldValues
(
volScalarFieldValue qPlant 0
- volScalarFieldValue plantCd 0
- volScalarFieldValue leafAreaDensity 0
+ volScalarFieldValue Cd 0
+ volScalarFieldValue LAD 0
);
regions
@@ -37,8 +37,8 @@ regions
box (0 0 0) (200 200 18.2192);
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
@@ -47,8 +47,8 @@ regions
patch bottom;
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyStable/system/setFieldsDict b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyStable/system/setFieldsDict
index 2b64fb1997275b874c5692ccbe5306fad9ff95a5..aa8aecfb3fefdce0d007f199ce228c0a8135cc6d 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyStable/system/setFieldsDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyStable/system/setFieldsDict
@@ -17,8 +17,8 @@ FoamFile
defaultFieldValues
(
volScalarFieldValue qPlant 0
- volScalarFieldValue plantCd 0
- volScalarFieldValue leafAreaDensity 0
+ volScalarFieldValue Cd 0
+ volScalarFieldValue LAD 0
);
regions
@@ -37,8 +37,8 @@ regions
box (0 0 0) (200 200 18.2192);
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
@@ -47,8 +47,8 @@ regions
patch bottom;
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyUnstable/system/setFieldsDict b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyUnstable/system/setFieldsDict
index 49e7e54314fbf9ab96f2740b127638c36ecbbbf6..906afbfdccbbd2c403043dc107910e7ee7e55c1a 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyUnstable/system/setFieldsDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/slightlyUnstable/system/setFieldsDict
@@ -17,8 +17,8 @@ FoamFile
defaultFieldValues
(
volScalarFieldValue qPlant 0
- volScalarFieldValue plantCd 0
- volScalarFieldValue leafAreaDensity 0
+ volScalarFieldValue Cd 0
+ volScalarFieldValue LAD 0
);
regions
@@ -37,8 +37,8 @@ regions
box (0 0 0) (200 200 18.2192);
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
@@ -47,8 +47,8 @@ regions
patch bottom;
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/stable/system/setFieldsDict b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/stable/system/setFieldsDict
index 10593c0cac5b7c8cf6d33c0e489ab138102ce2a0..07a4f536e3703441ad6dfa9b86224132045675e4 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/stable/system/setFieldsDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/stable/system/setFieldsDict
@@ -17,8 +17,8 @@ FoamFile
defaultFieldValues
(
volScalarFieldValue qPlant 0
- volScalarFieldValue plantCd 0
- volScalarFieldValue leafAreaDensity 0
+ volScalarFieldValue Cd 0
+ volScalarFieldValue LAD 0
);
regions
@@ -37,8 +37,8 @@ regions
box (0 0 0) (200 200 18.2192);
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
@@ -47,8 +47,8 @@ regions
patch bottom;
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/unstable/system/setFieldsDict b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/unstable/system/setFieldsDict
index ed88dd4404ca19bb1ed8f308a5f0c4a3d7933cbd..54debfbeb78cbba21c7208213b304acb8e09ffb4 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/unstable/system/setFieldsDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/unstable/system/setFieldsDict
@@ -17,8 +17,8 @@ FoamFile
defaultFieldValues
(
volScalarFieldValue qPlant 0
- volScalarFieldValue plantCd 0
- volScalarFieldValue leafAreaDensity 0
+ volScalarFieldValue Cd 0
+ volScalarFieldValue LAD 0
);
regions
@@ -37,8 +37,8 @@ regions
box (0 0 0) (200 200 18.2192);
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
@@ -47,8 +47,8 @@ regions
patch bottom;
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
);
diff --git a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/veryStable/system/setFieldsDict b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/veryStable/system/setFieldsDict
index a9a5a0cb5c125fdee1c025d43b28595d15091651..0d41b0c5085172d60f2ddbb53af539b66184851f 100644
--- a/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/veryStable/system/setFieldsDict
+++ b/tutorials/verificationAndValidation/atmosphericModels/atmForestStability/setups.orig/veryStable/system/setFieldsDict
@@ -17,8 +17,8 @@ FoamFile
defaultFieldValues
(
volScalarFieldValue qPlant 0
- volScalarFieldValue plantCd 0
- volScalarFieldValue leafAreaDensity 0
+ volScalarFieldValue Cd 0
+ volScalarFieldValue LAD 0
);
regions
@@ -37,8 +37,8 @@ regions
box (0 0 0) (200 200 18.2192);
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
@@ -47,8 +47,8 @@ regions
patch bottom;
fieldValues
(
- volScalarFieldValue plantCd 0.2
- volScalarFieldValue leafAreaDensity 0.14
+ volScalarFieldValue Cd 0.2
+ volScalarFieldValue LAD 0.14
);
}
);