zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2012-11-02 17:30:58 +00:00
parent e8a129487b f9d4d7f661
commit 96f51164eb
184 changed files with 3440 additions and 5956 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
applications/solvers/combustion/chemFoam/chemFoam.C
View File

@ -69,14 +69,9 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "solveChemistry.H"
{
#include "YEqn.H"
#include "hEqn.H"
#include "pEqn.H"
}
#include "YEqn.H"
#include "hEqn.H"
#include "pEqn.H"
#include "output.H"

2
applications/solvers/combustion/chemFoam/hEqn.H
View File

@ -9,4 +9,6 @@
{
h[0] = h0 + integratedHeat;
}
thermo.correct();
}

1
applications/solvers/combustion/chemFoam/pEqn.H
View File

@ -1,5 +1,4 @@
{
thermo.correct();
rho = thermo.rho();
if (constProp == "volume")
{

23
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
View File

@ -60,14 +60,9 @@ int main(int argc, char *argv[])
#include "createFluidMeshes.H"
#include "createSolidMeshes.H"
#include "createPorousFluidRegions.H"
#include "createPorousSolidMeshes.H"
#include "createFluidFields.H"
#include "createSolidFields.H"
#include "createPorousFluidFields.H"
#include "createPorousSolidFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
@ -116,24 +111,6 @@ int main(int argc, char *argv[])
#include "solveFluid.H"
}
forAll(porousFluidRegions, i)
{
Info<< "\nSolving for fluid porous region "
<< porousFluidRegions[i].name() << endl;
#include "setPorousFluidFields.H"
#include "readPorousFluidRegionPIMPLEControls.H"
#include "solvePorousFluid.H"
}
forAll(porousSolidRegions, i)
{
Info<< "\nSolving for porous solid region "
<< porousSolidRegions[i].name() << endl;
#include "setPorousRegionSolidFields.H"
#include "readPorousSolidMultiRegionPIMPLEControls.H"
#include "solvePorousSolid.H"
}
forAll(solidRegions, i)
{
Info<< "\nSolving for solid region "

3
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/Make/options
View File

@ -1,8 +1,7 @@
EXE_INC = \
-Ifluid \
-Isolid \
-I./porousFluid \
-I./porousSolid \
-I../solid \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \

22
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
View File

@ -50,13 +50,9 @@ int main(int argc, char *argv[])
#include "createFluidMeshes.H"
#include "createSolidMeshes.H"
#include "createPorousFluidRegions.H"
#include "createPorousSolidMeshes.H"
#include "createFluidFields.H"
#include "createSolidFields.H"
#include "createPorousFluidFields.H"
#include "createPorousSolidFields.H"
#include "initContinuityErrs.H"
@ -74,24 +70,6 @@ int main(int argc, char *argv[])
#include "solveFluid.H"
}
forAll(porousFluidRegions, i)
{
Info<< "\nSolving for fluid porous region "
<< porousFluidRegions[i].name() << endl;
#include "setPorousFluidFields.H"
#include "readPorousFluidRegionSIMPLEControls.H"
#include "solvePorousFluid.H"
}
forAll(porousSolidRegions, i)
{
Info<< "\nSolving for porous solid region "
<< porousSolidRegions[i].name() << endl;
#include "setPorousRegionSolidFields.H"
#include "readPorousSolidMultiRegionSIMPLEControls.H"
#include "solvePorousSolid.H"
}
forAll(solidRegions, i)
{
Info<< "\nSolving for solid region "

11
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/UPorousFluidEqn.H
View File

@ -1,11 +0,0 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> porousUEqn
(
fvm::div(porousPhi, porousU)
+ turbPorous.divDevRhoReff(porousU)
+ porousSources(porousRho, porousU)
);
porousUEqn().relax();
solve(porousUEqn() == -fvc::grad(porousP));

181
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/createPorousFluidFields.H
View File

@ -1,181 +0,0 @@
// Initialise porous field pointer lists
PtrList<rhoThermo> thermoPorous(porousFluidRegions.size());
PtrList<volScalarField> rhoPorous(porousFluidRegions.size());
PtrList<volScalarField> kappaPorous(porousFluidRegions.size());
PtrList<volVectorField> UPorous(porousFluidRegions.size());
PtrList<surfaceScalarField> phiPorous(porousFluidRegions.size());
PtrList<compressible::turbulenceModel> turbulencePorous
(
porousFluidRegions.size()
);
PtrList<volScalarField> pPorous(porousFluidRegions.size());
List<scalar> initialMassFluidPorous(porousFluidRegions.size());
List<label> pRefCellFluidPorous(porousFluidRegions.size(),0);
List<scalar> pRefValueFluidPorous(porousFluidRegions.size(),0.0);
PtrList<dimensionedScalar> rhoMaxPorous(fluidRegions.size());
PtrList<dimensionedScalar> rhoMinPorous(fluidRegions.size());
PtrList<IObasicSourceList> heatPorousSources
(
porousFluidRegions.size()
);
forAll(porousFluidRegions, i)
{
Info<< "Reading fluid mesh thermophysical properties for porous "
<< porousFluidRegions[i].name() << nl << endl;
Info<< " Adding to thermoFluid porous\n" << endl;
thermoPorous.set
(
i,
rhoThermo::New(porousFluidRegions[i]).ptr()
);
Info<< " Adding to rhoPorous\n" << endl;
rhoPorous.set
(
i,
new volScalarField
(
IOobject
(
"rho",
runTime.timeName(),
porousFluidRegions[i],
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermoPorous[i].rho()
)
);
Info<< " Adding to UPorous\n" << endl;
UPorous.set
(
i,
new volVectorField
(
IOobject
(
"U",
runTime.timeName(),
porousFluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousFluidRegions[i]
)
);
Info<< " Adding to phiPorous\n" << endl;
phiPorous.set
(
i,
new surfaceScalarField
(
IOobject
(
"phi",
runTime.timeName(),
porousFluidRegions[i],
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rhoPorous[i]*UPorous[i])
& porousFluidRegions[i].Sf()
)
);
Info<< " Adding turbulence to porous\n" << endl;
turbulencePorous.set
(
i,
compressible::turbulenceModel::New
(
rhoPorous[i],
UPorous[i],
phiPorous[i],
thermoPorous[i]
).ptr()
);
Info<< " Adding to kappaFluid\n" << endl;
kappaPorous.set
(
i,
new volScalarField
(
IOobject
(
"kappaPorous",
runTime.timeName(),
porousFluidRegions[i],
IOobject::NO_READ,
IOobject::NO_WRITE
),
thermoPorous[i].Cp()*thermoPorous[i].alpha()
)
);
pPorous.set
(
i,
new volScalarField
(
IOobject
(
"p",
runTime.timeName(),
porousFluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousFluidRegions[i]
)
);
setRefCell
(
thermoPorous[i].p(),
pPorous[i],
porousFluidRegions[i].solutionDict().subDict("SIMPLE"),
pRefCellFluidPorous[i],
pRefValueFluidPorous[i]
);
rhoMaxPorous.set
(
i,
new dimensionedScalar
(
porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
(
"rhoMax"
)
)
);
rhoMinPorous.set
(
i,
new dimensionedScalar
(
porousFluidRegions[i].solutionDict().subDict("SIMPLE").lookup
(
"rhoMin"
)
)
);
heatPorousSources.set
(
i,
new IObasicSourceList(porousFluidRegions[i])
);
}

25
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/createPorousFluidRegions.H
View File

@ -1,25 +0,0 @@
const wordList porousFluidNames(rp["porousFluid"]);
PtrList<fvMesh> porousFluidRegions(porousFluidNames.size());
forAll (porousFluidNames, iPorous)
{
const word porousFluidName = porousFluidNames[iPorous];
Info<< "Create porous fluid region " << porousFluidName
<< nl << endl;
porousFluidRegions.set
(
iPorous,
new fvMesh
(
IOobject
(
porousFluidName,
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}

19
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/hPorousFluidEqn.H
View File

@ -1,19 +0,0 @@
{
fvScalarMatrix hPorousEqn
(
fvm::div(porousPhi, porousH)
- fvm::laplacian(turbPorous.alphaEff(), porousH)
==
- fvc::div(porousPhi, 0.5*magSqr(porousU), "div(phi,K)")
+ porousSources(porousRho, porousH)
);
hPorousEqn.relax();
hPorousEqn.solve();
porousThermo.correct();
Info<< "Min/max in the porous T:"
<< min(porousThermo.T()).value() << ' '
<< max(porousThermo.T()).value() << endl;
}

53
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/pPorousFluidEqn.H
View File

@ -1,53 +0,0 @@
porousRho = porousThermo.rho();
porousRho = max(porousRho, rhoMin);
porousRho = min(porousRho, rhoMax);
porousRho.relax();
volScalarField rAUPorous(1.0/porousUEqn().A());
porousU = rAUPorous*porousUEqn().H();
porousUEqn.clear();
bool closedVolume = false;
porousPhi =
fvc::interpolate(porousRho)
*(fvc::interpolate(porousU) & porousMesh.Sf());
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(porousRho*rAUPorous, porousP) == fvc::div(porousPhi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
porousPhi -= pEqn.flux();
}
}
porousP.relax();
porousU -= rAUPorous*fvc::grad(porousP);
porousU.correctBoundaryConditions();
if (closedVolume)
{
porousP += (initialMass - fvc::domainIntegrate(porousPsi*porousP))
/fvc::domainIntegrate(porousPsi);
}
porousRho = porousThermo.rho();
porousRho = max(porousRho, rhoMin);
porousRho = min(porousRho, rhoMax);
porousRho.relax();
Info<< "rho max/min : "
<< max(porousRho).value() << " "
<< min(porousRho).value() << endl;

4
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/readPorousFluidRegionSIMPLEControls.H
View File

@ -1,4 +0,0 @@
const dictionary& simple = porousMesh.solutionDict().subDict("SIMPLE");
const int nNonOrthCorr =
simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

28
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/setPorousFluidFields.H
View File

@ -1,28 +0,0 @@
const fvMesh& porousMesh = porousFluidRegions[i];
rhoThermo& porousThermo = thermoPorous[i];
volScalarField& porousRho = rhoPorous[i];
volVectorField& porousU = UPorous[i];
surfaceScalarField& porousPhi = phiPorous[i];
compressible::turbulenceModel& turbPorous = turbulencePorous[i];
volScalarField& porousP = porousThermo.p();
const volScalarField& porousPsi = porousThermo.psi();
volScalarField& porousH = porousThermo.he();
const dimensionedScalar initialMass
(
"initialMass",
dimMass,
initialMassFluidPorous[i]
);
IObasicSourceList& porousSources = heatPorousSources[i];
const label pRefCell = pRefCellFluidPorous[i];
const scalar pRefValue = pRefValueFluidPorous[i];
const scalar rhoMax = rhoMaxPorous[i].value();
const scalar rhoMin = rhoMinPorous[i].value();

11
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/solvePorousFluid.H
View File

@ -1,11 +0,0 @@
// Pressure-velocity SIMPLE corrector
porousP.storePrevIter();
porousRho.storePrevIter();
{
#include "UPorousFluidEqn.H"
#include "hPorousFluidEqn.H"
#include "pPorousFluidEqn.H"
}
turbPorous.correct();

43
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/createPorousSolidFields.H
View File

@ -1,43 +0,0 @@
// Initialise solid field pointer lists
PtrList<solidThermo> porousSolidThermos(porousSolidRegions.size());
PtrList<IObasicSourceList> solidHeatSources(porousSolidRegions.size());
PtrList<volScalarField> betavSolid(porousSolidRegions.size());
// Populate solid field pointer lists
forAll(porousSolidRegions, i)
{
Info<< "*** Reading porous solid mesh thermophysical "
<< "properties for region "
<< porousSolidRegions[i].name() << nl << endl;
Info<< " Adding to thermos\n" << endl;
porousSolidThermos.set
(
i,
solidThermo::New(porousSolidRegions[i])
);
Info<< " Adding sources\n" << endl;
solidHeatSources.set
(
i,
new IObasicSourceList(porousSolidRegions[i])
);
Info<< " Adding to betavSolid\n" << endl;
betavSolid.set
(
i,
new volScalarField
(
IOobject
(
"betavSolid",
runTime.timeName(),
porousSolidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousSolidRegions[i]
)
);
}

24
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/createPorousSolidMeshes.H
View File

@ -1,24 +0,0 @@
const wordList porousSolidNames(rp["porousSolid"]);
PtrList<fvMesh> porousSolidRegions(porousSolidNames.size());
forAll(porousSolidNames, i)
{
Info<< "Create solid mesh for region " << porousSolidNames[i]
<< " for time = " << runTime.timeName() << nl << endl;
porousSolidRegions.set
(
i,
new fvMesh
(
IOobject
(
porousSolidNames[i],
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}

4
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/readPorousSolidMultiRegionSIMPLEControls.H
View File

@ -1,4 +0,0 @@
const dictionary& simple = mesh.solutionDict().subDict("SIMPLE");
const int nNonOrthCorr =
simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

23
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/setPorousRegionSolidFields.H
View File

@ -1,23 +0,0 @@
const fvMesh& mesh = porousSolidRegions[i];
solidThermo& thermo = porousSolidThermos[i];
const volScalarField& betav = betavSolid[i];
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
tmp<volScalarField> tcp = thermo.Cp();
const volScalarField& cp = tcp();
tmp<volScalarField> tkappa = thermo.kappa();
//tmp<volSymmTensorField> tkappa = thermo.directionalKappa()*betav;
const volScalarField& kappa = tkappa();
//const volSymmTensorField& K = tK();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
volScalarField& h = thermo.he();
IObasicSourceList& sources = solidHeatSources[i];

18
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solidPorousRegionDiffusionNo.H
View File

@ -1,18 +0,0 @@
scalar DiNum = -GREAT;
forAll(solidRegions, i)
{
# include "setRegionSolidFields.H"
DiNum = max
(
solidRegionDiffNo
(
solidRegions[i],
runTime,
rho*cp,
K
),
DiNum
);
}

17
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousSolid/solvePorousSolid.H
View File

@ -1,17 +0,0 @@
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> hEqn
(
- fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+ sources(rho, h)
);
hEqn().relax();
hEqn().solve();
}
}
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;

20
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/createSolidFields.H
View File

@ -1,20 +0,0 @@
// Initialise solid field pointer lists
PtrList<solidThermo> thermos(solidRegions.size());
PtrList<radiation::radiationModel> radiations(solidRegions.size());
// Populate solid field pointer lists
forAll(solidRegions, i)
{
Info<< "*** Reading solid mesh thermophysical properties for region "
<< solidRegions[i].name() << nl << endl;
Info<< " Adding to thermos\n" << endl;
thermos.set
(
i,
solidThermo::New(solidRegions[i])
);
Info<< " Adding to radiations\n" << endl;
radiations.set(i, radiation::radiationModel::New(thermos[i].T()));
}

24
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/createSolidMeshes.H
View File

@ -1,24 +0,0 @@
const wordList solidsNames(rp["solid"]);
PtrList<fvMesh> solidRegions(solidsNames.size());
forAll(solidsNames, i)
{
Info<< "Create solid mesh for region " << solidsNames[i]
<< " for time = " << runTime.timeName() << nl << endl;
solidRegions.set
(
i,
new fvMesh
(
IOobject
(
solidsNames[i],
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}

18
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H
View File

@ -1,18 +0,0 @@
fvMesh& mesh = solidRegions[i];
solidThermo& thermo = thermos[i];
const radiation::radiationModel& radiation = radiations[i];
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
tmp<volScalarField> tcp = thermo.Cp();
const volScalarField& cp = tcp();
tmp<volScalarField> tkappa = thermo.kappa();
//tmp<volSymmTensorField> tkappa = thermo.directionalkappa();
const volScalarField& kappa = tkappa();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
volScalarField& h = thermo.he();

3
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H
View File

@ -3,7 +3,8 @@
{
fvScalarMatrix hEqn
(
-fvm::laplacian(alpha, h)
- fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+ sources(rho, h)
);
hEqn.relax();
hEqn.solve();

3
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
View File

@ -14,7 +14,7 @@
CoNum
);
}
/*
forAll (porousFluidRegions, porousI)
{
CoNum = max
@ -29,3 +29,4 @@
CoNum
);
}
*/

24
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/UPorousFluidEqn.H
View File

@ -1,24 +0,0 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> porousUEqn
(
fvm::ddt(porousRho, porousU)
+ fvm::div(porousPhi, porousU)
+ turbPorous.divDevRhoReff(porousU)
+ porousSources(porousRho, porousU)
);
porousUEqn().relax();
volScalarField rAUPorous(1.0/porousUEqn().A());
if (momentumPredictor)
{
solve(porousUEqn() == -fvc::grad(porousP));
}
else
{
porousU = rAUPorous*(porousUEqn().H() - fvc::grad(porousP));
porousU.correctBoundaryConditions();
}

147
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/createPorousFluidFields.H
View File

@ -1,147 +0,0 @@
// Initialise porous field pointer lists
PtrList<rhoThermo> thermoPorous(porousFluidRegions.size());
PtrList<volScalarField> rhoPorous(porousFluidRegions.size());
PtrList<volVectorField> UPorous(porousFluidRegions.size());
PtrList<surfaceScalarField> phiPorous(porousFluidRegions.size());
PtrList<volScalarField> KPorous(porousFluidRegions.size());
PtrList<volScalarField> dpdtPorous(fluidRegions.size());
PtrList<compressible::turbulenceModel> turbulencePorous
(
porousFluidRegions.size()
);
PtrList<volScalarField> pPorous(porousFluidRegions.size());
PtrList<IObasicSourceList> heatPorousSources
(
porousFluidRegions.size()
);
forAll(porousFluidRegions, i)
{
Info<< "Reading fluid mesh thermophysical properties for porous "
<< porousFluidRegions[i].name() << nl << endl;
Info<< " Adding to thermoFluid porous\n" << endl;
thermoPorous.set
(
i,
rhoThermo::New(porousFluidRegions[i]).ptr()
);
Info<< " Adding to rhoPorous\n" << endl;
rhoPorous.set
(
i,
new volScalarField
(
IOobject
(
"rho",
runTime.timeName(),
porousFluidRegions[i],
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermoPorous[i].rho()
)
);
Info<< " Adding to UPorous\n" << endl;
UPorous.set
(
i,
new volVectorField
(
IOobject
(
"U",
runTime.timeName(),
porousFluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousFluidRegions[i]
)
);
Info<< " Adding to phiPorous\n" << endl;
phiPorous.set
(
i,
new surfaceScalarField
(
IOobject
(
"phi",
runTime.timeName(),
porousFluidRegions[i],
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rhoPorous[i]*UPorous[i])
& porousFluidRegions[i].Sf()
)
);
Info<< " Adding turbulence to porous\n" << endl;
turbulencePorous.set
(
i,
compressible::turbulenceModel::New
(
rhoPorous[i],
UPorous[i],
phiPorous[i],
thermoPorous[i]
).ptr()
);
Info<< " Adding to KPorous\n" << endl;
KPorous.set
(
i,
new volScalarField
(
"KPorous",
0.5*magSqr(UPorous[i])
)
);
Info<< " Adding to dpdtPorous\n" << endl;
dpdtPorous.set
(
i,
new volScalarField
(
"dpdtPorous",
fvc::ddt(thermoPorous[i].p())
)
);
pPorous.set
(
i,
new volScalarField
(
IOobject
(
"p",
runTime.timeName(),
porousFluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousFluidRegions[i]
)
);
heatPorousSources.set
(
i,
new IObasicSourceList(porousFluidRegions[i])
);
}

25
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/createPorousFluidRegions.H
View File

@ -1,25 +0,0 @@
const wordList porousFluidNames(rp["porousFluid"]);
PtrList<fvMesh> porousFluidRegions(porousFluidNames.size());
forAll (porousFluidNames, iPorous)
{
const word porousFluidName = porousFluidNames[iPorous];
Info<< "Create porous fluid region " << porousFluidName
<< nl << endl;
porousFluidRegions.set
(
iPorous,
new fvMesh
(
IOobject
(
porousFluidName,
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}

21
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/hPorousFluidEqn.H
View File

@ -1,21 +0,0 @@
{
fvScalarMatrix hPorousEqn
(
fvm::ddt(porousRho, porousH)
+ fvm::div(porousPhi, porousH)
- fvm::laplacian(turbPorous.alphaEff(), porousH)
==
porousdpdt
- (fvc::ddt(porousRho, porousK) + fvc::div(porousPhi, porousK))
+ porousSources(porousRho, porousH)
);
hPorousEqn.relax();
hPorousEqn.solve();
porousThermo.correct();
Info<< "Min/max in the porous T:"
<< min(porousThermo.T()).value() << ' '
<< max(porousThermo.T()).value() << endl;
}

63
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/pPorousFluidEqn.H
View File

@ -1,63 +0,0 @@
porousRho = porousThermo.rho();
porousU = rAUPorous*porousUEqn().H();
if (nCorr <= 1)
{
porousUEqn.clear();
}
porousPhi =
fvc::interpolate(porousRho)*
(
(fvc::interpolate(porousU) & porousMesh.Sf())
+ fvc::ddtPhiCorr(rAUPorous, porousRho, porousU, porousPhi)
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::ddt(porousPsi, porousP)
+ fvc::div(porousPhi)
- fvm::laplacian(porousRho*rAUPorous, porousP)
);
pEqn.solve
(
porousMesh.solver
(
porousP.select
(
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
);
if (nonOrth == nNonOrthCorr)
{
porousPhi += pEqn.flux();
}
}
solve(fvm::ddt(porousRho) + fvc::div(porousPhi));
// Explicitly relax pressure for momentum corrector
porousP.relax();
// Recalculate density from the relaxed pressure
porousRho = porousThermo.rho();
porousU -= rAUPorous*fvc::grad(porousP);
porousU.correctBoundaryConditions();
porousK = 0.5*magSqr(porousU);
// Update pressure time derivative if needed
if (porousThermo.dpdt())
{
porousdpdt = fvc::ddt(porousP);
}

11
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/readPorousFluidRegionPIMPLEControls.H
View File

@ -1,11 +0,0 @@
const dictionary& pimple = porousMesh.solutionDict().subDict("PIMPLE");
const int nCorr =
pimple.lookupOrDefault<int>("nCorrectors", 1);
const int nNonOrthCorr =
pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
const bool momentumPredictor =
pimple.lookupOrDefault("momentumPredictor", true);

17
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/setPorousFluidFields.H
View File

@ -1,17 +0,0 @@
fvMesh& porousMesh = porousFluidRegions[i];
rhoThermo& porousThermo = thermoPorous[i];
volScalarField& porousRho = rhoPorous[i];
volVectorField& porousU = UPorous[i];
surfaceScalarField& porousPhi = phiPorous[i];
compressible::turbulenceModel& turbPorous = turbulencePorous[i];
volScalarField& porousK = KPorous[i];
volScalarField& porousdpdt = dpdtPorous[i];
volScalarField& porousP = porousThermo.p();
const volScalarField& porousPsi = porousThermo.psi();
volScalarField& porousH = porousThermo.he();
IObasicSourceList& porousSources = heatPorousSources[i];

28
applications/solvers/heatTransfer/chtMultiRegionFoam/porousFluid/solvePorousFluid.H
View File

@ -1,28 +0,0 @@
if (finalIter)
{
porousMesh.data::add("finalIteration", true);
}
if (oCorr == 0)
{
solve(fvm::ddt(porousRho) + fvc::div(porousPhi));
}
#include "UPorousFluidEqn.H"
#include "hPorousFluidEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pPorousFluidEqn.H"
}
turbPorous.correct();
porousRho = porousThermo.rho();
if (finalIter)
{
porousMesh.data::remove("finalIteration");
}

43
applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/createPorousSolidFields.H
View File

@ -1,43 +0,0 @@
// Initialise solid field pointer lists
PtrList<solidThermo> porousSolidThermos(porousSolidRegions.size());
PtrList<IObasicSourceList> solidHeatSources(porousSolidRegions.size());
PtrList<volScalarField> betavSolid(porousSolidRegions.size());
// Populate solid field pointer lists
forAll(porousSolidRegions, i)
{
Info<< "*** Reading porous solid mesh thermophysical "
<< "properties for region "
<< porousSolidRegions[i].name() << nl << endl;
Info<< " Adding to thermos\n" << endl;
porousSolidThermos.set
(
i,
solidThermo::New(porousSolidRegions[i])
);
Info<< " Adding sources\n" << endl;
solidHeatSources.set
(
i,
new IObasicSourceList(porousSolidRegions[i])
);
Info<< " Adding to betavSolid\n" << endl;
betavSolid.set
(
i,
new volScalarField
(
IOobject
(
"betavSolid",
runTime.timeName(),
porousSolidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
porousSolidRegions[i]
)
);
}

24
applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/createPorousSolidMeshes.H
View File

@ -1,24 +0,0 @@
const wordList porousSolidNames(rp["porousSolid"]);
PtrList<fvMesh> porousSolidRegions(porousSolidNames.size());
forAll(porousSolidNames, i)
{
Info<< "Create solid mesh for region " << porousSolidNames[i]
<< " for time = " << runTime.timeName() << nl << endl;
porousSolidRegions.set
(
i,
new fvMesh
(
IOobject
(
porousSolidNames[i],
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}

4
applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/readPorousSolidMultiRegionPIMPLEControls.H
View File

@ -1,4 +0,0 @@
const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
int nNonOrthCorr =
pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

26
applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/setPorousRegionSolidFields.H
View File

@ -1,26 +0,0 @@
fvMesh& mesh = porousSolidRegions[i];
solidThermo& thermo = porousSolidThermos[i];
const volScalarField& betav = betavSolid[i];
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
tmp<volScalarField> tcp = thermo.Cp();
const volScalarField& cp = tcp();
tmp<volScalarField> tkappa = thermo.kappa();
//tmp<volSymmTensorField> tkappa = thermo.directionalKappa()*betav;
const volScalarField& kappa = tkappa();
//const volSymmTensorField& K = tK();
//tmp<volScalarField> trhoCp = cp*rho;
//const volScalarField& rhoCp = trhoCp();
tmp<volScalarField> talpha = thermo.alpha();
const volScalarField& alpha = talpha();
volScalarField& h = thermo.he();
IObasicSourceList& sources = solidHeatSources[i];

18
applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solidPorousRegionDiffusionNo.H
View File

@ -1,18 +0,0 @@
scalar DiNum = -GREAT;
forAll(solidRegions, i)
{
# include "setRegionSolidFields.H"
DiNum = max
(
solidRegionDiffNo
(
solidRegions[i],
runTime,
rho*cp,
K
),
DiNum
);
}

28
applications/solvers/heatTransfer/chtMultiRegionFoam/porousSolid/solvePorousSolid.H
View File

@ -1,28 +0,0 @@
if (finalIter)
{
mesh.data::add("finalIteration", true);
}
{
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> hEqn
(
fvm::ddt(betav*rho, h)
- fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+ sources(rho, h)
);
hEqn().relax();
hEqn().solve(mesh.solver(h.select(finalIter)));
}
}
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T()) << endl;
if (finalIter)
{
mesh.data::remove("finalIteration");
}

47
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
View File

@ -1,6 +1,8 @@
// Initialise solid field pointer lists
PtrList<solidThermo> thermos(solidRegions.size());
PtrList<radiation::radiationModel> radiations(solidRegions.size());
PtrList<IObasicSourceList> solidHeatSources(solidRegions.size());
PtrList<volScalarField> betavSolid(solidRegions.size());
// Populate solid field pointer lists
forAll(solidRegions, i)
@ -13,4 +15,49 @@
Info<< " Adding to radiations\n" << endl;
radiations.set(i, radiation::radiationModel::New(thermos[i].T()));
Info<< " Adding sources\n" << endl;
solidHeatSources.set
(
i,
new IObasicSourceList(solidRegions[i])
);
IOobject betavSolidIO
(
"betavSolid",
runTime.timeName(),
solidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
if (betavSolidIO.headerOk())
{
betavSolid.set
(
i,
new volScalarField(betavSolidIO, solidRegions[i])
);
}
else
{
betavSolid.set
(
i,
new volScalarField
(
IOobject
(
"betavSolid",
runTime.timeName(),
solidRegions[i],
IOobject::NO_READ,
IOobject::NO_WRITE
),
solidRegions[i],
dimensionedScalar("1", dimless, scalar(1.0))
)
);
}
}

4
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H
View File

@ -14,3 +14,7 @@
const volScalarField& kappa = tkappa();
volScalarField& h = thermo.he();
const volScalarField& betav = betavSolid[i];
IObasicSourceList& sources = solidHeatSources[i];

5
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
View File

@ -8,8 +8,9 @@ if (finalIter)
{
tmp<fvScalarMatrix> hEqn
(
fvm::ddt(rho, h)
- fvm::laplacian(alpha, h)
fvm::ddt(betav*rho, h)
- fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+ sources(rho, h)
);
hEqn().relax();
hEqn().solve(mesh.solver(h.select(finalIter)));

3
applications/test/thermoMixture/Make/files Normal file
View File

@ -0,0 +1,3 @@
Test-thermoMixture.C
EXE = $(FOAM_USER_APPBIN)/Test-thermoMixture

5
applications/test/thermoMixture/Make/options Normal file
View File

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
EXE_LIBS = \
-lspecie

67
src/thermophysicalModels/solidThermo/mixtures/pureSolidMixture/pureSolidMixture.C → applications/test/thermoMixture/Test-thermoMixture.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,48 +21,59 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
ThermoMixture
Description
\*---------------------------------------------------------------------------*/
#include "pureSolidMixture.H"
#include "fvMesh.H"
#include "dictionary.H"
#include "IFstream.H"
#include "specie.H"
#include "perfectGas.H"
#include "hConstThermo.H"
#include "sensibleEnthalpy.H"
#include "thermo.H"
#include "constTransport.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
namespace Foam
int main(int argc, char *argv[])
{
typedef constTransport
<
species::thermo
<
hConstThermo<perfectGas<specie> >,
sensibleEnthalpy
>
> ThermoType;
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
dictionary dict(IFstream("thermoDict")());
template<class ThermoType>
pureSolidMixture<ThermoType>::pureSolidMixture
(
const dictionary& thermoDict,
const fvMesh& mesh
)
:
basicMixture(thermoDict, mesh),
mixture_(thermoDict.subDict("mixture"))
{}
ThermoType t1(dict.subDict("specie1"));
ThermoType t2(dict.subDict("specie2"));
Info<< "Checking Cp of mixture of hConstThermo" << endl;
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Info<< "W 1, 2, (1 + 2) = " << t1.W() << " " << t2.W() << " "
<< (t1 + t2).W() << endl;
template<class ThermoType>
pureSolidMixture<ThermoType>::~pureSolidMixture()
{}
Info<< "Cp 1, 2, 1 + 2 = " << t1.cp(1, 1) << " " << t2.cp(1, 1) << " "
<< (t1 + t2).cp(1, 1) << endl;
ThermoType t3(t1);
t3 += t2;
Info<< "Cp (1 += 2) = " << t3.cp(1, 1) << endl;
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Info<< "\nEnd\n" << endl;
template<class ThermoType>
void pureSolidMixture<ThermoType>::read(const dictionary& thermoDict)
{
mixture_ = ThermoType(thermoDict.subDict("mixture"));
return 0;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

43
applications/test/thermoMixture/thermoDict Normal file
View File

@ -0,0 +1,43 @@
specie1
{
specie
{
nMoles 1;
molWeight 1;
}
thermodynamics
{
Cp 1;
Cv 1;
Hf 0;
}
transport
{
mu 1;
Pr 1;
}
}
specie2
{
specie
{
nMoles 1;
molWeight 0.5;
}
thermodynamics
{
Cp 2;
Cv 2;
Hf 0;
}
transport
{
mu 1;
Pr 1;
}
}

168
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
View File

@ -96,20 +96,29 @@ castellatedMeshControls
// refinement.
nCellsBetweenLevels 1;
// Explicit feature edge refinement
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Specifies a level for any cell intersected by explicitly provided
// edges.
// This is a featureEdgeMesh, read from constant/triSurface for now.
// Specify 'levels' in the same way as the 'distance' mode in the
// refinementRegions (see below). The old specification
// level 2;
// is equivalent to
// levels ((0 2));
features
(
//{
// file "someLine.eMesh";
// level 2;
// //level 2;
// levels ((0.0 2) (1.0 3));
//}
);
// Surface based refinement
// ~~~~~~~~~~~~~~~~~~~~~~~~
@ -178,7 +187,7 @@ castellatedMeshControls
// three modes
// - distance. 'levels' specifies per distance to the surface the
// wanted refinement level. The distances need to be specified in
// descending order.
// increasing order.
// - inside. 'levels' is only one entry and only the level is used. All
// cells inside the surface get refined up to the level. The surface
// needs to be closed for this to be possible.
@ -217,37 +226,37 @@ castellatedMeshControls
// Settings for the snapping.
snapControls
{
//- Number of patch smoothing iterations before finding correspondence
// to surface
// Number of patch smoothing iterations before finding correspondence
// to surface
nSmoothPatch 3;
//- Maximum relative distance for points to be attracted by surface.
// True distance is this factor times local maximum edge length.
// Maximum relative distance for points to be attracted by surface.
// True distance is this factor times local maximum edge length.
// Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
tolerance 2.0;
//- Number of mesh displacement relaxation iterations.
// Number of mesh displacement relaxation iterations.
nSolveIter 30;
//- Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
// Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
nRelaxIter 5;
// Feature snapping
//- Number of feature edge snapping iterations.
// Leave out altogether to disable.
// Number of feature edge snapping iterations.
// Leave out altogether to disable.
nFeatureSnapIter 10;
//- Detect (geometric only) features by sampling the surface
// (default=false).
// Detect (geometric only) features by sampling the surface
// (default=false).
implicitFeatureSnap false;
//- Use castellatedMeshControls::features (default = true)
// Use castellatedMeshControls::features (default = true)
explicitFeatureSnap true;
//- Detect features between multiple surfaces
// (only for explicitFeatureSnap, default = false)
// Detect features between multiple surfaces
// (only for explicitFeatureSnap, default = false)
multiRegionFeatureSnap false;
}
@ -258,9 +267,43 @@ addLayersControls
// size of the refined cell outside layer (true) or absolute sizes (false).
relativeSizes true;
// Layer thickness specification. This can be specified in one of four ways
// - expansionRatio and finalLayerThickness (cell nearest internal mesh)
// - expansionRatio and firstLayerThickness (cell on surface)
// - overall thickness and firstLayerThickness
// - overall thickness and finalLayerThickness
// Expansion factor for layer mesh
expansionRatio 1.0;
// Wanted thickness of the layer furthest away from the wall.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
finalLayerThickness 0.3;
// Wanted thickness of the layer next to the wall.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
//firstLayerThickness 0.3;
// Wanted overall thickness of layers.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
//thickness 0.5
// Minimum overall thickness of total layers. If for any reason layer
// cannot be above minThickness do not add layer.
// If relativeSizes this is relative to undistorted size of cell
// outside layer..
minThickness 0.25;
// Per final patch (so not geometry!) the layer information
// Note: This behaviour changed after 21x. Any non-mentioned patches
// now slide unless nSurfaceLayers is explicitly mentioned to be 0.
// now slide unless:
// - nSurfaceLayers is explicitly mentioned to be 0.
// - angle to nearest surface < slipFeatureAngle (see below)
layers
{
sphere.stl_firstSolid
@ -272,9 +315,9 @@ addLayersControls
{
nSurfaceLayers 1;
// Per patch layer data
expansionRatio 1.3;
expansionRatio 1.3;
finalLayerThickness 0.3;
minThickness 0.1;
minThickness 0.1;
}
// Disable any mesh shrinking and layer addition on any point of
@ -285,41 +328,24 @@ addLayersControls
}
}
// Expansion factor for layer mesh
expansionRatio 1.0;
//- Wanted thickness of final added cell layer. If multiple layers
// is the
// thickness of the layer furthest away from the wall.
// Relative to undistorted size of cell outside layer.
// is the thickness of the layer furthest away from the wall.
// See relativeSizes parameter.
finalLayerThickness 0.3;
//- Minimum thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// Relative to undistorted size of cell outside layer.
// See relativeSizes parameter.
minThickness 0.25;
//- If points get not extruded do nGrow layers of connected faces that are
// also not grown. This helps convergence of the layer addition process
// close to features.
// If points get not extruded do nGrow layers of connected faces that are
// also not grown. This helps convergence of the layer addition process
// close to features.
// Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
nGrow 0;
// Advanced settings
//- When not to extrude surface. 0 is flat surface, 90 is when two faces
// make straight angle.
// When not to extrude surface. 0 is flat surface, 90 is when two faces
// make straight angle.
featureAngle 60;
//- At non-patched sides allow mesh to slip if extrusion direction makes
// angle larger than slipFeatureAngle.
// At non-patched sides allow mesh to slip if extrusion direction makes
// angle larger than slipFeatureAngle.
slipFeatureAngle 30;
//- Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
// Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
nRelaxIter 5;
// Number of smoothing iterations of surface normals
@ -365,51 +391,51 @@ addLayersControls
// where to undo.
meshQualityControls
{
//- Maximum non-orthogonality allowed. Set to 180 to disable.
// Maximum non-orthogonality allowed. Set to 180 to disable.
maxNonOrtho 65;
//- Max skewness allowed. Set to <0 to disable.
// Max skewness allowed. Set to <0 to disable.
maxBoundarySkewness 20;
maxInternalSkewness 4;
//- Max concaveness allowed. Is angle (in degrees) below which concavity
// is allowed. 0 is straight face, <0 would be convex face.
// Set to 180 to disable.
// Max concaveness allowed. Is angle (in degrees) below which concavity
// is allowed. 0 is straight face, <0 would be convex face.
// Set to 180 to disable.
maxConcave 80;
//- Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
// Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;
//- Minimum quality of the tet formed by the face-centre
// and variable base point minimum decomposition triangles and
// the cell centre. This has to be a positive number for tracking
// to work. Set to very negative number (e.g. -1E30) to
// disable.
// <0 = inside out tet,
// 0 = flat tet
// 1 = regular tet
// Minimum quality of the tet formed by the face-centre
// and variable base point minimum decomposition triangles and
// the cell centre. This has to be a positive number for tracking
// to work. Set to very negative number (e.g. -1E30) to
// disable.
// <0 = inside out tet,
// 0 = flat tet
// 1 = regular tet
minTetQuality 1e-9;
//- Minimum face area. Set to <0 to disable.
// Minimum face area. Set to <0 to disable.
minArea -1;
//- Minimum face twist. Set to <-1 to disable. dot product of face normal
//- and face centre triangles normal
// Minimum face twist. Set to <-1 to disable. dot product of face normal
// and face centre triangles normal
minTwist 0.05;
//- minimum normalised cell determinant
//- 1 = hex, <= 0 = folded or flattened illegal cell
// minimum normalised cell determinant
// 1 = hex, <= 0 = folded or flattened illegal cell
minDeterminant 0.001;
//- minFaceWeight (0 -> 0.5)
// minFaceWeight (0 -> 0.5)
minFaceWeight 0.05;
//- minVolRatio (0 -> 1)
// minVolRatio (0 -> 1)
minVolRatio 0.01;
//must be >0 for Fluent compatibility
// must be >0 for Fluent compatibility
minTriangleTwist -1;
//- if >0 : preserve single cells with all points on the surface if the
@ -420,9 +446,9 @@ meshQualityControls
// Advanced
//- Number of error distribution iterations
// Number of error distribution iterations
nSmoothScale 4;
//- amount to scale back displacement at error points
// amount to scale back displacement at error points
errorReduction 0.75;
// Optional : some meshing phases allow usage of relaxed rules.

11
applications/utilities/postProcessing/sampling/sample/sampleDict
View File

@ -90,6 +90,8 @@ fields
// uniform: extra number of sampling points
// polyLine, cloud: list of coordinates
// patchCloud: list of coordinates and set of patches to look for nearest
// patchSeed: random sampling on set of patches. Points slightly off
// face centre.
sets
(
lineX1
@ -135,6 +137,15 @@ sets
maxDistance 0.1; // maximum distance to search
patches (".*Wall.*");
}
patchSeed
{
patches (".*Wall.*");
// Number of points to seed. Divided amongst all processors according
// to fraction of patches they hold.
maxPoints 100;
}
);

414
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View File

@ -318,7 +318,6 @@ void Foam::autoLayerDriver::handleNonManifolds
}
}
Info<< "Set displacement to zero for all " << nNonManif
<< " non-manifold points" << endl;
}
@ -443,86 +442,6 @@ void Foam::autoLayerDriver::handleFeatureAngle
}
//Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::undistortedEdgeLength
//(
// const indirectPrimitivePatch& pp,
// const bool relativeSizes,
// const bool faceSize
//)
//{
// const fvMesh& mesh = meshRefiner_.mesh();
//
// tmp<scalarField> tfld(new scalarField());
// scalarField& fld = tfld();
//
// if (faceSize)
// {
// fld.setSize(pp.size());
// }
// else
// {
// fld.setSize(pp.nPoints());
// }
//
//
// if (relativeSizes)
// {
// const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
// const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
//
// if (faceSize)
// {
// forAll(pp, i)
// {
// label faceI = pp.addressing()[i];
// label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
// fld[i] = edge0Len/(1<<ownLevel);
// }
// }
// else
// {
// // Determine per point the max cell level of connected cells
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// labelList maxPointLevel(pp.nPoints(), labelMin);
//
// forAll(pp, i)
// {
// label faceI = pp.addressing()[i];
// label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
//
// const face& f = pp.localFaces()[i];
// forAll(f, fp)
// {
// maxPointLevel[f[fp]] =
// max(maxPointLevel[f[fp]], ownLevel);
// }
// }
//
// syncTools::syncPointList
// (
// mesh,
// pp.meshPoints(),
// maxPointLevel,
// maxEqOp<label>(),
// labelMin // null value
// );
//
//
// forAll(maxPointLevel, pointI)
// {
// // Find undistorted edge size for this level.
// fld[i] = edge0Len/(1<<maxPointLevel[pointI]);
// }
// }
// }
// else
// {
// // Use actual cell size
// }
//}
// No extrusion on cells with warped faces. Calculates the thickness of the
// layer and compares it to the space the warped face takes up. Disables
// extrusion if layer thickness is more than faceRatio of the thickness of
@ -702,7 +621,6 @@ void Foam::autoLayerDriver::handleWarpedFaces
//}
// No extrusion on faces with differing number of layers for points
void Foam::autoLayerDriver::setNumLayers
(
const labelList& patchToNLayers,
@ -865,15 +783,6 @@ Foam::autoLayerDriver::makeLayerDisplacementField
}
//Pout<< "*** makeLayerDisplacementField : boundary conditions:" << endl;
//forAll(patchFieldTypes, patchI)
//{
// Pout<< "\t" << patchI << " name:" << pointPatches[patchI].name()
// << " type:" << patchFieldTypes[patchI]
// << " nLayers:" << numLayers[patchI]
// << endl;
//}
const polyMesh& mesh = pMesh();
// Note: time().timeName() instead of meshRefinement::timeName() since
@ -1053,11 +962,10 @@ void Foam::autoLayerDriver::determineSidePatches
patchDict.add("nFaces", 0);
patchDict.add("startFace", mesh.nFaces());
Pout<< "Adding patch " << patchI
<< " name:" << name
<< " between " << Pstream::myProcNo()
<< " and " << nbrProcI
<< endl;
//Pout<< "Adding patch " << patchI
// << " name:" << name
// << " between " << Pstream::myProcNo()
// << " and " << nbrProcI << endl;
label procPatchI = meshRefiner_.appendPatch
(
@ -1090,12 +998,7 @@ void Foam::autoLayerDriver::calculateLayerThickness
(
const indirectPrimitivePatch& pp,
const labelList& patchIDs,
const scalarField& patchExpansionRatio,
const bool relativeSizes,
const scalarField& patchFinalLayerThickness,
const scalarField& patchMinThickness,
const layerParameters& layerParams,
const labelList& cellLevel,
const labelList& patchNLayers,
const scalar edge0Len,
@ -1111,12 +1014,13 @@ void Foam::autoLayerDriver::calculateLayerThickness
// Rework patch-wise layer parameters into minimum per point
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Note: only layer parameters consistent with layer specification
// method (see layerParameters) will be correct.
scalarField firstLayerThickness(pp.nPoints(), GREAT);
scalarField finalLayerThickness(pp.nPoints(), GREAT);
scalarField totalThickness(pp.nPoints(), GREAT);
scalarField expRatio(pp.nPoints(), GREAT);
// Reuse input fields
expansionRatio.setSize(pp.nPoints());
expansionRatio = GREAT;
thickness.setSize(pp.nPoints());
thickness = GREAT;
minThickness.setSize(pp.nPoints());
minThickness = GREAT;
@ -1130,20 +1034,30 @@ void Foam::autoLayerDriver::calculateLayerThickness
{
label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
expansionRatio[ppPointI] = min
firstLayerThickness[ppPointI] = min
(
expansionRatio[ppPointI],
patchExpansionRatio[patchI]
firstLayerThickness[ppPointI],
layerParams.firstLayerThickness()[patchI]
);
thickness[ppPointI] = min
finalLayerThickness[ppPointI] = min
(
thickness[ppPointI],
patchFinalLayerThickness[patchI]
finalLayerThickness[ppPointI],
layerParams.finalLayerThickness()[patchI]
);
totalThickness[ppPointI] = min
(
totalThickness[ppPointI],
layerParams.thickness()[patchI]
);
expRatio[ppPointI] = min
(
expRatio[ppPointI],
layerParams.expansionRatio()[patchI]
);
minThickness[ppPointI] = min
(
minThickness[ppPointI],
patchMinThickness[patchI]
layerParams.minThickness()[patchI]
);
}
}
@ -1152,7 +1066,7 @@ void Foam::autoLayerDriver::calculateLayerThickness
(
mesh,
pp.meshPoints(),
expansionRatio,
firstLayerThickness,
minEqOp<scalar>(),
GREAT // null value
);
@ -1160,7 +1074,23 @@ void Foam::autoLayerDriver::calculateLayerThickness
(
mesh,
pp.meshPoints(),
thickness,
finalLayerThickness,
minEqOp<scalar>(),
GREAT // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
totalThickness,
minEqOp<scalar>(),
GREAT // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
expRatio,
minEqOp<scalar>(),
GREAT // null value
);
@ -1182,14 +1112,18 @@ void Foam::autoLayerDriver::calculateLayerThickness
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// by multiplying with the internal cell size.
if (relativeSizes)
if (layerParams.relativeSizes())
{
if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2)
if
(
min(layerParams.minThickness()) < 0
|| max(layerParams.minThickness()) > 2
)
{
FatalErrorIn("calculateLayerThickness(..)")
<< "Thickness should be factor of local undistorted cell size."
<< " Valid values are [0..2]." << nl
<< " minThickness:" << patchMinThickness
<< " minThickness:" << layerParams.minThickness()
<< exit(FatalError);
}
@ -1225,38 +1159,114 @@ void Foam::autoLayerDriver::calculateLayerThickness
{
// Find undistorted edge size for this level.
scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
thickness[pointI] *= edgeLen;
firstLayerThickness[pointI] *= edgeLen;
finalLayerThickness[pointI] *= edgeLen;
totalThickness[pointI] *= edgeLen;
minThickness[pointI] *= edgeLen;
}
}
// Rework thickness (of final layer) into overall thickness of all layers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Rework thickness parameters into overall thickness
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(thickness, pointI)
forAll(firstLayerThickness, pointI)
{
// Calculate layer thickness based on expansion ratio
// and final layer height
if (expansionRatio[pointI] == 1)
{
thickness[pointI] *= patchNLayers[pointI];
}
else
{
thickness[pointI] = layerParams.layerThickness
(
patchNLayers[pointI],
firstLayerThickness[pointI],
finalLayerThickness[pointI],
totalThickness[pointI],
expRatio[pointI]
);
scalar invExpansion = 1.0 / expansionRatio[pointI];
label nLay = patchNLayers[pointI];
thickness[pointI] *=
(1.0 - pow(invExpansion, nLay))
/ (1.0 - invExpansion);
}
expansionRatio[pointI] = layerParams.layerExpansionRatio
(
patchNLayers[pointI],
firstLayerThickness[pointI],
finalLayerThickness[pointI],
totalThickness[pointI],
expRatio[pointI]
);
}
//Info<< "calculateLayerThickness : min:" << gMin(thickness)
// << " max:" << gMax(thickness) << endl;
// Print a bit
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Find maximum length of a patch name, for a nicer output
label maxPatchNameLen = 0;
forAll(patchIDs, i)
{
label patchI = patchIDs[i];
word patchName = patches[patchI].name();
maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
}
Info<< nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
<< setw(0) << " faces layers avg thickness[m]" << nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << " "
<< setw(0) << " near-wall overall" << nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
<< setw(0) << " ----- ------ --------- -------" << endl;
forAll(patchIDs, i)
{
label patchI = patchIDs[i];
const labelList& meshPoints = patches[patchI].meshPoints();
scalar sumThickness = 0;
scalar sumNearWallThickness = 0;
forAll(meshPoints, patchPointI)
{
label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
sumThickness += thickness[ppPointI];
sumNearWallThickness += layerParams.firstLayerThickness
(
patchNLayers[ppPointI],
firstLayerThickness[ppPointI],
finalLayerThickness[ppPointI],
thickness[ppPointI],
expansionRatio[ppPointI]
);
}
label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
// For empty patches, totNPoints is 0.
scalar avgThickness = 0;
scalar avgNearWallThickness = 0;
if (totNPoints > 0)
{
avgThickness =
returnReduce(sumThickness, sumOp<scalar>())
/ totNPoints;
avgNearWallThickness =
returnReduce(sumNearWallThickness, sumOp<scalar>())
/ totNPoints;
}
Info<< setf(ios_base::left) << setw(maxPatchNameLen)
<< patches[patchI].name() << setprecision(3)
<< " " << setw(8)
<< returnReduce(patches[patchI].size(), sumOp<scalar>())
<< " " << setw(6) << layerParams.numLayers()[patchI]
<< " " << setw(8) << avgNearWallThickness
<< " " << setw(8) << avgThickness
<< endl;
}
Info<< endl;
}
}
@ -2618,7 +2628,8 @@ void Foam::autoLayerDriver::addLayers
const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
// Determine (wanted) point-wise layer thickness and expansion ratio
// Determine (wanted) point-wise overall layer thickness and expansion
// ratio
scalarField thickness(pp().nPoints());
scalarField minThickness(pp().nPoints());
scalarField expansionRatio(pp().nPoints());
@ -2626,12 +2637,7 @@ void Foam::autoLayerDriver::addLayers
(
pp,
meshMover().adaptPatchIDs(),
layerParams.expansionRatio(),
layerParams.relativeSizes(), // thickness relative to cellsize?
layerParams.finalLayerThickness(), // wanted thicknes
layerParams.minThickness(), // minimum thickness
layerParams,
cellLevel,
patchNLayers,
edge0Len,
@ -2642,87 +2648,6 @@ void Foam::autoLayerDriver::addLayers
);
// Print a bit
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Find maximum length of a patch name, for a nicer output
label maxPatchNameLen = 0;
forAll(meshMover().adaptPatchIDs(), i)
{
label patchI = meshMover().adaptPatchIDs()[i];
word patchName = patches[patchI].name();
maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
}
Info<< nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
<< setw(0) << " faces layers avg thickness[m]" << nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << " "
<< setw(0) << " near-wall overall" << nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
<< setw(0) << " ----- ------ --------- -------" << endl;
forAll(meshMover().adaptPatchIDs(), i)
{
label patchI = meshMover().adaptPatchIDs()[i];
const labelList& meshPoints = patches[patchI].meshPoints();
scalar sumThickness = 0;
scalar sumNearWallThickness = 0;
forAll(meshPoints, patchPointI)
{
label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
sumThickness += thickness[ppPointI];
label nLay = patchNLayers[ppPointI];
if (nLay > 0)
{
if (expansionRatio[ppPointI] == 1)
{
sumNearWallThickness += thickness[ppPointI]/nLay;
}
else
{
scalar s =
(1.0-pow(expansionRatio[ppPointI], nLay))
/ (1.0-expansionRatio[ppPointI]);
sumNearWallThickness += thickness[ppPointI]/s;
}
}
}
label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
// For empty patches, totNPoints is 0.
scalar avgThickness = 0;
scalar avgNearWallThickness = 0;
if (totNPoints > 0)
{
avgThickness =
returnReduce(sumThickness, sumOp<scalar>())
/ totNPoints;
avgNearWallThickness =
returnReduce(sumNearWallThickness, sumOp<scalar>())
/ totNPoints;
}
Info<< setf(ios_base::left) << setw(maxPatchNameLen)
<< patches[patchI].name() << setprecision(3)
<< " " << setw(8)
<< returnReduce(patches[patchI].size(), sumOp<scalar>())
<< " " << setw(6) << layerParams.numLayers()[patchI]
<< " " << setw(8) << avgNearWallThickness
<< " " << setw(8) << avgThickness
<< endl;
}
Info<< endl;
}
// Calculate wall to medial axis distance for smoothing displacement
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2958,8 +2883,8 @@ void Foam::autoLayerDriver::addLayers
// Determine per point/per face number of layers to extrude. Also
// handles the slow termination of layers when going switching layers
labelList nPatchPointLayers(pp().nPoints(),-1);
labelList nPatchFaceLayers(pp().localFaces().size(),-1);
labelList nPatchPointLayers(pp().nPoints(), -1);
labelList nPatchFaceLayers(pp().size(), -1);
setupLayerInfoTruncation
(
meshMover(),
@ -2970,31 +2895,22 @@ void Foam::autoLayerDriver::addLayers
nPatchFaceLayers
);
// Calculate displacement for first layer for addPatchLayer.
// (first layer = layer of cells next to the original mesh)
vectorField firstDisp(patchNLayers.size(), vector::zero);
// Calculate displacement for final layer for addPatchLayer.
// (layer of cells next to the original mesh)
vectorField finalDisp(patchNLayers.size(), vector::zero);
forAll(nPatchPointLayers, i)
{
if (nPatchPointLayers[i] > 0)
{
if (expansionRatio[i] == 1.0)
{
firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
}
else
{
label nLay = nPatchPointLayers[i];
scalar h =
pow(expansionRatio[i], nLay - 1)
* (1.0 - expansionRatio[i])
/ (1.0 - pow(expansionRatio[i], nLay));
firstDisp[i] = h*patchDisp[i];
}
}
scalar ratio = layerParams.finalLayerThicknessRatio
(
nPatchPointLayers[i],
expansionRatio[i]
);
finalDisp[i] = ratio*patchDisp[i];
}
const scalarField invExpansionRatio(1.0 / expansionRatio);
const scalarField invExpansionRatio(1.0/expansionRatio);
// Add topo regardless of whether extrudeStatus is extruderemove.
// Not add layer if patchDisp is zero.
@ -3009,7 +2925,7 @@ void Foam::autoLayerDriver::addLayers
labelList(0), // exposed patchIDs, not used for adding layers
nPatchFaceLayers, // layers per face
nPatchPointLayers, // layers per point
firstDisp, // thickness of layer nearest internal mesh
finalDisp, // thickness of layer nearest internal mesh
meshMod
);

7
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
View File

@ -237,12 +237,7 @@ class autoLayerDriver
(
const indirectPrimitivePatch& pp,
const labelList& patchIDs,
const scalarField& patchExpansionRatio,
const bool relativeSizes,
const scalarField& patchFinalLayerThickness,
const scalarField& patchMinThickness,
const layerParameters& layerParams,
const labelList& cellLevel,
const labelList& patchNLayers,
const scalar edge0Len,

3
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
View File

@ -97,6 +97,7 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
refineParams.curvature(),
true, // featureRefinement
false, // featureDistanceRefinement
false, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
@ -207,6 +208,7 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
refineParams.curvature(),
false, // featureRefinement
false, // featureDistanceRefinement
false, // internalRefinement
true, // surfaceRefinement
true, // curvatureRefinement
@ -368,6 +370,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
refineParams.curvature(),
false, // featureRefinement
true, // featureDistanceRefinement
true, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement

11
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
View File

@ -1406,9 +1406,10 @@ void Foam::autoSnapDriver::doSnap
adaptPatchIDs
)
);
indirectPrimitivePatch& pp = ppPtr();
// Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
const scalarField snapDist(calcSnapDistance(snapParams, pp));
// Construct iterative mesh mover.
@ -1420,7 +1421,7 @@ void Foam::autoSnapDriver::doSnap
motionSmoother meshMover
(
mesh,
ppPtr(),
pp,
adaptPatchIDs,
meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
motionDict
@ -1475,7 +1476,7 @@ void Foam::autoSnapDriver::doSnap
}
// Check for displacement being outwards.
outwardsDisplacement(ppPtr(), disp);
outwardsDisplacement(pp, disp);
// Set initial distribution of displacement field (on patches)
// from patchDisp and make displacement consistent with b.c.
@ -1489,8 +1490,8 @@ void Foam::autoSnapDriver::doSnap
(
mesh.time().path()
/ "patchDisplacement_" + name(iter) + ".obj",
ppPtr().localPoints(),
ppPtr().localPoints() + disp
pp.localPoints(),
pp.localPoints() + disp
);
}

26
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
View File

@ -141,14 +141,6 @@ class autoSnapDriver
vectorField& pointSurfaceNormal,
vectorField& pointRotation
) const;
//void calcNearestFace
//(
// const label iter,
// const indirectPrimitivePatch& pp,
// vectorField& faceDisp,
// vectorField& faceSurfaceNormal,
// vectorField& faceRotation
//) const;
void calcNearestFace
(
const label iter,
@ -158,15 +150,19 @@ class autoSnapDriver
labelList& faceSurfaceRegion,
vectorField& faceRotation
) const;
void interpolateFaceToPoint
void calcNearestFacePointProperties
(
const label iter,
const indirectPrimitivePatch& pp,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceRotation,
const List<List<point> >& pointFaceCentres,
vectorField& patchDisp
//vectorField& patchRotationDisp
const vectorField& faceDisp,
const vectorField& faceSurfaceNormal,
const labelList& faceSurfaceRegion,
List<List<point> >& pointFaceSurfNormals,
List<List<point> >& pointFaceDisp,
List<List<point> >& pointFaceCentres,
List<labelList>& pointFacePatchID
) const;
void correctAttraction
(
@ -276,6 +272,7 @@ class autoSnapDriver
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch&,
const scalarField&,
@ -339,6 +336,7 @@ class autoSnapDriver
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,

510
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
View File

@ -321,7 +321,7 @@ void Foam::autoSnapDriver::calcNearestFace
const indirectPrimitivePatch& pp,
vectorField& faceDisp,
vectorField& faceSurfaceNormal,
labelList& faceSurfaceRegion,
labelList& faceSurfaceGlobalRegion,
vectorField& faceRotation
) const
{
@ -333,8 +333,8 @@ void Foam::autoSnapDriver::calcNearestFace
faceDisp = vector::zero;
faceSurfaceNormal.setSize(pp.size());
faceSurfaceNormal = vector::zero;
faceSurfaceRegion.setSize(pp.size());
faceSurfaceRegion = -1;
faceSurfaceGlobalRegion.setSize(pp.size());
faceSurfaceGlobalRegion = -1;
// Divide surfaces into zoned and unzoned
labelList zonedSurfaces = surfaces.getNamedSurfaces();
@ -419,7 +419,7 @@ void Foam::autoSnapDriver::calcNearestFace
label faceI = ppFaces[hitI];
faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
faceSurfaceNormal[faceI] = hitNormal[hitI];
faceSurfaceRegion[faceI] = surfaces.globalRegion
faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
(
hitSurface[hitI],
hitRegion[hitI]
@ -477,7 +477,11 @@ void Foam::autoSnapDriver::calcNearestFace
label faceI = ppFaces[hitI];
faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
faceSurfaceNormal[faceI] = hitNormal[hitI];
faceSurfaceRegion[faceI] = hitRegion[hitI];
faceSurfaceGlobalRegion[faceI] = surfaces.globalRegion
(
hitSurface[hitI],
hitRegion[hitI]
);
}
}
@ -517,6 +521,169 @@ void Foam::autoSnapDriver::calcNearestFace
}
// Collect (possibly remote) per point data of all surrounding faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - faceSurfaceNormal
// - faceDisp
// - faceCentres&faceNormal
void Foam::autoSnapDriver::calcNearestFacePointProperties
(
const label iter,
const indirectPrimitivePatch& pp,
const vectorField& faceDisp,
const vectorField& faceSurfaceNormal,
const labelList& faceSurfaceGlobalRegion,
List<List<point> >& pointFaceSurfNormals,
List<List<point> >& pointFaceDisp,
List<List<point> >& pointFaceCentres,
List<labelList>& pointFacePatchID
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
// For now just get all surrounding face data. Expensive - should just
// store and sync data on coupled points only
// (see e.g PatchToolsNormals.C)
pointFaceSurfNormals.setSize(pp.nPoints());
pointFaceDisp.setSize(pp.nPoints());
pointFaceCentres.setSize(pp.nPoints());
pointFacePatchID.setSize(pp.nPoints());
// Fill local data
forAll(pp.pointFaces(), pointI)
{
const labelList& pFaces = pp.pointFaces()[pointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
pNormals.setSize(pFaces.size());
List<point>& pDisp = pointFaceDisp[pointI];
pDisp.setSize(pFaces.size());
List<point>& pFc = pointFaceCentres[pointI];
pFc.setSize(pFaces.size());
labelList& pFid = pointFacePatchID[pointI];
pFid.setSize(pFaces.size());
forAll(pFaces, i)
{
label faceI = pFaces[i];
pNormals[i] = faceSurfaceNormal[faceI];
pDisp[i] = faceDisp[faceI];
pFc[i] = pp.faceCentres()[faceI];
//label meshFaceI = pp.addressing()[faceI];
//pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
pFid[i] = globalToPatch_[faceSurfaceGlobalRegion[faceI]];
}
}
// Collect additionally 'normal' boundary faces for boundaryPoints of pp
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// points on the boundary of pp should pick up non-pp normals
// as well for the feature-reconstruction to behave correctly.
// (the movement is already constrained outside correctly so it
// is only that the unconstrained attraction vector is calculated
// correctly)
{
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
labelList patchID(pbm.patchID());
// Unmark all non-coupled boundary faces
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (pp.coupled() || isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
label meshFaceI = pp.start()+i;
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
}
}
// Remove any meshed faces
forAll(pp.addressing(), i)
{
label meshFaceI = pp.addressing()[i];
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
// See if pp point uses any non-meshed boundary faces
const labelList& boundaryPoints = pp.boundaryPoints();
forAll(boundaryPoints, i)
{
label pointI = boundaryPoints[i];
label meshPointI = pp.meshPoints()[pointI];
const point& pt = mesh.points()[meshPointI];
const labelList& pFaces = mesh.pointFaces()[meshPointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
List<point>& pDisp = pointFaceDisp[pointI];
List<point>& pFc = pointFaceCentres[pointI];
labelList& pFid = pointFacePatchID[pointI];
forAll(pFaces, i)
{
label meshFaceI = pFaces[i];
if (!mesh.isInternalFace(meshFaceI))
{
label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
if (patchI != -1)
{
vector fn = mesh.faceAreas()[meshFaceI];
pNormals.append(fn/mag(fn));
pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
pFc.append(mesh.faceCentres()[meshFaceI]);
pFid.append(patchI);
}
}
}
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceSurfNormals,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceDisp,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceCentres,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFacePatchID,
listPlusEqOp<label>(),
List<label>()
);
}
// Gets passed in offset to nearest point on feature edge. Calculates
// if the point has a different number of faces on either side of the feature
// and if so attracts the point to that non-dominant plane.
@ -580,56 +747,7 @@ Foam::pointIndexHit Foam::autoSnapDriver::findMultiPatchPoint
}
return pointIndexHit(false, vector::zero, labelMax);
}
////XXXXXXXX
//void Foam::autoSnapDriver::attractMultiPatchPoint
//(
// const label iter,
// const scalar featureCos,
//
// const indirectPrimitivePatch& pp,
// const scalarField& snapDist,
// const label pointI,
//
// const List<List<point> >& pointFaceSurfNormals,
// const labelListList& pointFaceSurfaceRegion,
// const List<List<point> >& pointFaceDisp,
// const List<List<point> >& pointFaceCentres,
// const labelListList& pointFacePatchID,
//
// vector& patchAttraction,
// pointConstraint& patchConstraint
//) const
//{
// // Collect
//
// );
//
// if
// (
// (constraint.first() > patchConstraints[pointI].first())
// || (magSqr(attraction) < magSqr(patchAttraction[pointI]))
// )
// {
// patchAttraction[pointI] = attraction;
// patchConstraints[pointI] = constraint;
//
// // Check the number of directions
// if (patchConstraints[pointI].first() == 1)
// {
// // Flat surface. Check for different patchIDs
// pointIndexHit multiPatchPt
// (
// findMultiPatchPoint
// (
// pt,
// pointFacePatchID[pointI],
// pointFaceCentres[pointI]
// )
// );
// if (multiPatchPt.hit())
// {
// // Behave like when having two surface normals so
////XXXXXXXX
void Foam::autoSnapDriver::binFeatureFace
(
@ -1361,6 +1479,7 @@ void Foam::autoSnapDriver::determineFeatures
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
@ -1464,69 +1583,72 @@ void Foam::autoSnapDriver::determineFeatures
if (patchConstraints[pointI].first() == 1)
{
// Flat surface. Check for different patchIDs
pointIndexHit multiPatchPt
(
findMultiPatchPoint
(
pt,
pointFacePatchID[pointI],
pointFaceCentres[pointI]
)
);
if (multiPatchPt.hit())
if (multiRegionFeatureSnap)
{
// Behave like when having two surface normals so
// attract to nearest feature edge (with a guess for
// the multipatch point as starting point)
label featI = -1;
pointIndexHit nearInfo = findNearFeatureEdge
pointIndexHit multiPatchPt
(
pp,
snapDist,
pointI,
multiPatchPt.hitPoint(), //estimatedPt
featI,
edgeAttractors,
edgeConstraints,
patchAttraction,
patchConstraints
findMultiPatchPoint
(
pt,
pointFacePatchID[pointI],
pointFaceCentres[pointI]
)
);
if (multiPatchPt.hit())
{
// Behave like when having two surface normals so
// attract to nearest feature edge (with a guess for
// the multipatch point as starting point)
label featI = -1;
pointIndexHit nearInfo = findNearFeatureEdge
(
pp,
snapDist,
pointI,
multiPatchPt.hitPoint(), //estimatedPt
if (nearInfo.hit())
{
// Dump
if (featureEdgeStr.valid())
featI,
edgeAttractors,
edgeConstraints,
patchAttraction,
patchConstraints
);
if (nearInfo.hit())
{
meshTools::writeOBJ(featureEdgeStr(), pt);
featureEdgeVertI++;
meshTools::writeOBJ
(
featureEdgeStr(),
nearInfo.hitPoint()
);
featureEdgeVertI++;
featureEdgeStr()
<< "l " << featureEdgeVertI-1 << ' '
<< featureEdgeVertI << nl;
// Dump
if (featureEdgeStr.valid())
{
meshTools::writeOBJ(featureEdgeStr(), pt);
featureEdgeVertI++;
meshTools::writeOBJ
(
featureEdgeStr(),
nearInfo.hitPoint()
);
featureEdgeVertI++;
featureEdgeStr()
<< "l " << featureEdgeVertI-1 << ' '
<< featureEdgeVertI << nl;
}
}
}
else
{
if (missedEdgeStr.valid())
else
{
meshTools::writeOBJ(missedEdgeStr(), pt);
missedVertI++;
meshTools::writeOBJ
(
missedEdgeStr(),
nearInfo.missPoint()
);
missedVertI++;
missedEdgeStr()
<< "l " << missedVertI-1 << ' '
<< missedVertI << nl;
if (missedEdgeStr.valid())
{
meshTools::writeOBJ(missedEdgeStr(), pt);
missedVertI++;
meshTools::writeOBJ
(
missedEdgeStr(),
nearInfo.missPoint()
);
missedVertI++;
missedEdgeStr()
<< "l " << missedVertI-1 << ' '
<< missedVertI << nl;
}
}
}
}
@ -1645,6 +1767,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
(
const label iter,
const scalar featureCos,
const bool multiRegionFeatureSnap,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
@ -1697,6 +1820,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
(
iter,
featureCos,
multiRegionFeatureSnap,
pp,
snapDist,
@ -2131,6 +2255,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
//MEJ: any faces that have multi-patch points only keep the multi-patch
// points. The other points on the face will be dragged along
// (hopefully)
if (multiRegionFeatureSnap)
{
autoPtr<OFstream> multiPatchStr;
if (debug&meshRefinement::OBJINTERSECTIONS)
@ -2536,10 +2661,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
{
const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
const Switch explicitFeatureAttraction = snapParams.explicitFeatureSnap();
const Switch multiRegionFeatureSnap = snapParams.multiRegionFeatureSnap();
Info<< "Overriding displacement on features :" << nl
<< " implicit features : " << implicitFeatureAttraction << nl
<< " explicit features : " << explicitFeatureAttraction << nl
<< " implicit features : " << implicitFeatureAttraction << nl
<< " explicit features : " << explicitFeatureAttraction << nl
<< " multi-patch features : " << multiRegionFeatureSnap << nl
<< endl;
@ -2547,6 +2674,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
const pointField& localPoints = pp.localPoints();
const fvMesh& mesh = meshRefiner_.mesh();
// Displacement and orientation per pp face
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2554,7 +2682,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
vectorField faceDisp(pp.size(), vector::zero);
// normal of surface at point on surface
vectorField faceSurfaceNormal(pp.size(), vector::zero);
labelList faceSurfaceRegion(pp.size(), -1);
labelList faceSurfaceGlobalRegion(pp.size(), -1);
vectorField faceRotation(pp.size(), vector::zero);
calcNearestFace
@ -2563,7 +2691,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
pp,
faceDisp,
faceSurfaceNormal,
faceSurfaceRegion,
faceSurfaceGlobalRegion,
faceRotation
);
@ -2587,158 +2715,25 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - faceSurfaceNormal
// - faceDisp
// - faceRotation
// - faceCentres&faceNormal
// For now just get all surrounding face data. Expensive - should just
// store and sync data on coupled points only
// (see e.g PatchToolsNormals.C)
List<List<point> > pointFaceSurfNormals(pp.nPoints());
List<List<point> > pointFaceDisp(pp.nPoints());
//List<List<point> > pointFaceRotation(pp.nPoints());
List<List<point> > pointFaceCentres(pp.nPoints());
List<labelList> pointFacePatchID(pp.nPoints());
// Fill local data
forAll(pp.pointFaces(), pointI)
{
const labelList& pFaces = pp.pointFaces()[pointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
pNormals.setSize(pFaces.size());
List<point>& pDisp = pointFaceDisp[pointI];
pDisp.setSize(pFaces.size());
//List<point>& pRot = pointFaceRotation[pointI];
//pRot.setSize(pFaces.size());
List<point>& pFc = pointFaceCentres[pointI];
pFc.setSize(pFaces.size());
labelList& pFid = pointFacePatchID[pointI];
pFid.setSize(pFaces.size());
forAll(pFaces, i)
{
label faceI = pFaces[i];
pNormals[i] = faceSurfaceNormal[faceI];
pDisp[i] = faceDisp[faceI];
//pRot[i] = faceRotation[faceI];
pFc[i] = pp.faceCentres()[faceI];
label meshFaceI = pp.addressing()[faceI];
pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI);
}
}
// Collect additionally 'normal' boundary faces for boundaryPoints of pp
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// points on the boundary of pp should pick up non-pp normals
// as well for the feature-reconstruction to behave correctly.
// (the movement is already constrained outside correctly so it
// is only that the unconstrained attraction vector is calculated
// correctly)
{
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
labelList patchID(pbm.patchID());
// Unmark all non-coupled boundary faces
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (pp.coupled() || isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
label meshFaceI = pp.start()+i;
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
}
}
// Remove any meshed faces
forAll(pp.addressing(), i)
{
label meshFaceI = pp.addressing()[i];
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
// See if pp point uses any non-meshed boundary faces
const labelList& boundaryPoints = pp.boundaryPoints();
forAll(boundaryPoints, i)
{
label pointI = boundaryPoints[i];
label meshPointI = pp.meshPoints()[pointI];
const point& pt = mesh.points()[meshPointI];
const labelList& pFaces = mesh.pointFaces()[meshPointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
List<point>& pDisp = pointFaceDisp[pointI];
List<point>& pFc = pointFaceCentres[pointI];
labelList& pFid = pointFacePatchID[pointI];
forAll(pFaces, i)
{
label meshFaceI = pFaces[i];
if (!mesh.isInternalFace(meshFaceI))
{
label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
if (patchI != -1)
{
vector fn = mesh.faceAreas()[meshFaceI];
pNormals.append(fn/mag(fn));
pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
pFc.append(mesh.faceCentres()[meshFaceI]);
pFid.append(patchI);
}
}
}
}
}
syncTools::syncPointList
calcNearestFacePointProperties
(
mesh,
pp.meshPoints(),
iter,
pp,
faceDisp,
faceSurfaceNormal,
faceSurfaceGlobalRegion,
pointFaceSurfNormals,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceDisp,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
//syncTools::syncPointList
//(
// mesh,
// pp.meshPoints(),
// pointFaceRotation,
// listPlusEqOp<point>(),
// List<point>(),
// listTransform()
//);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceCentres,
listPlusEqOp<point>(),
List<point>(),
listTransform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFacePatchID,
listPlusEqOp<label>(),
List<label>()
pointFacePatchID
);
@ -2793,6 +2788,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
(
iter,
featureCos,
multiRegionFeatureSnap,
pp,
snapDist,

440
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
View File

@ -34,6 +34,72 @@ License
const Foam::scalar Foam::layerParameters::defaultConcaveAngle = 90;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::scalar Foam::layerParameters::layerExpansionRatio
(
const label n,
const scalar totalOverFirst
) const
{
if (n <= 1)
{
return 1.0;
}
//scalar totalOverFirst = totalThickness/firstLayerThickess;
const label maxIters = 10;
const scalar tol = 1e-8;
if (mag(n-totalOverFirst) < tol)
{
return 1.0;
}
// Calculate the bounds of the solution
scalar minR;
scalar maxR;
if (totalOverFirst < n)
{
minR = 0.0;
maxR = pow(totalOverFirst/n, 1/(n-1));
}
else
{
minR = pow(totalOverFirst/n, 1/(n-1));
maxR = totalOverFirst/(n - 1);
}
//Info<< "Solution bounds = (" << minR << ", " << maxR << ")" << nl << endl;
// Starting guess
scalar r = 0.5*(minR + maxR);
for (label i = 0; i < maxIters; ++i)
{
const scalar prevr = r;
const scalar fx = pow(r, n) - totalOverFirst*r - (1 - totalOverFirst);
const scalar dfx = n*pow(r, n - 1) - totalOverFirst;
r -= fx/dfx;
const scalar error = mag(r - prevr);
//Info<< i << " " << r << " Error = " << error << endl;
if (error < tol)
{
break;
}
}
return r;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
@ -44,17 +110,12 @@ Foam::layerParameters::layerParameters
)
:
numLayers_(boundaryMesh.size(), -1),
expansionRatio_
(
boundaryMesh.size(),
readScalar(dict.lookup("expansionRatio"))
),
relativeSizes_(dict.lookup("relativeSizes")),
finalLayerThickness_
(
boundaryMesh.size(),
readScalar(dict.lookup("finalLayerThickness"))
),
layerSpec_(ILLEGAL),
firstLayerThickness_(boundaryMesh.size(), -123),
finalLayerThickness_(boundaryMesh.size(), -123),
thickness_(boundaryMesh.size(), -123),
expansionRatio_(boundaryMesh.size(), -123),
minThickness_
(
boundaryMesh.size(),
@ -103,24 +164,108 @@ Foam::layerParameters::layerParameters
nRelaxedIter_(labelMax),
additionalReporting_(dict.lookupOrDefault("additionalReporting", false))
{
if (nGrow_ > 0)
// Detect layer specification mode
label nSpec = 0;
bool haveFirst = dict.found("firstLayerThickness");
if (haveFirst)
{
WarningIn("layerParameters::layerParameters(..)")
<< "The nGrow parameter effect has changed with respect to 1.6.x."
<< endl
<< "Please set nGrow=0 for 1.6.x behaviour."
firstLayerThickness_ = scalarField
(
boundaryMesh.size(),
readScalar(dict.lookup("firstLayerThickness"))
);
nSpec++;
}
bool haveFinal = dict.found("finalLayerThickness");
if (haveFinal)
{
finalLayerThickness_ = scalarField
(
boundaryMesh.size(),
readScalar(dict.lookup("finalLayerThickness"))
);
nSpec++;
}
bool haveTotal = dict.found("thickness");
if (haveTotal)
{
thickness_ = scalarField
(
boundaryMesh.size(),
readScalar(dict.lookup("thickness"))
);
nSpec++;
}
bool haveExp = dict.found("expansionRatio");
if (haveExp)
{
expansionRatio_ = scalarField
(
boundaryMesh.size(),
readScalar(dict.lookup("expansionRatio"))
);
nSpec++;
}
if (haveFirst && haveTotal)
{
layerSpec_ = FIRST_AND_TOTAL;
Info<< "Layer thickness specified as first layer and overall thickness."
<< endl;
}
else if (haveFirst && haveExp)
{
layerSpec_ = FIRST_AND_EXPANSION;
Info<< "Layer thickness specified as first layer and expansion ratio."
<< endl;
}
else if (haveFinal && haveTotal)
{
layerSpec_ = FINAL_AND_TOTAL;
Info<< "Layer thickness specified as final layer and overall thickness."
<< endl;
}
else if (haveFinal && haveExp)
{
layerSpec_ = FINAL_AND_EXPANSION;
Info<< "Layer thickness specified as final layer and expansion ratio."
<< endl;
}
if (layerSpec_ == ILLEGAL || nSpec != 2)
{
FatalIOErrorIn
(
"layerParameters::layerParameters"
"(const dictionary&, const polyBoundaryMesh&)",
dict
) << "Over- or underspecified layer thickness."
<< " Please specify" << nl
<< " first layer thickness ('firstLayerThickness')"
<< " and overall thickness ('thickness') or" << nl
<< " first layer thickness ('firstLayerThickness')"
<< " and expansion ratio ('expansionRatio') or" << nl
<< " final layer thickness ('finalLayerThickness')"
<< " and expansion ratio ('expansionRatio') or" << nl
<< " final layer thickness ('finalLayerThickness')"
<< " and overall thickness ('thickness')"
<< exit(FatalIOError);
}
dict.readIfPresent("nRelaxedIter", nRelaxedIter_);
if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
{
FatalErrorIn("layerParameters::layerParameters(..)")
FatalIOErrorIn("layerParameters::layerParameters(..)", dict)
<< "Layer iterations should be >= 0." << endl
<< "nLayerIter:" << nLayerIter_
<< " nRelaxedIter:" << nRelaxedIter_
<< exit(FatalError);
<< exit(FatalIOError);
}
@ -130,7 +275,7 @@ Foam::layerParameters::layerParameters
{
if (iter().isDict())
{
const word& key = iter().keyword();
const keyType& key = iter().keyword();
const labelHashSet patchIDs
(
boundaryMesh.patchSet(List<wordRe>(1, key))
@ -154,16 +299,69 @@ Foam::layerParameters::layerParameters
numLayers_[patchI] =
readLabel(layerDict.lookup("nSurfaceLayers"));
layerDict.readIfPresent
(
"expansionRatio",
expansionRatio_[patchI]
);
layerDict.readIfPresent
(
"finalLayerThickness",
finalLayerThickness_[patchI]
);
switch (layerSpec_)
{
case FIRST_AND_TOTAL:
layerDict.readIfPresent
(
"firstLayerThickness",
firstLayerThickness_[patchI]
);
layerDict.readIfPresent
(
"thickness",
thickness_[patchI]
);
break;
case FIRST_AND_EXPANSION:
layerDict.readIfPresent
(
"firstLayerThickness",
firstLayerThickness_[patchI]
);
layerDict.readIfPresent
(
"expansionRatio",
expansionRatio_[patchI]
);
break;
case FINAL_AND_TOTAL:
layerDict.readIfPresent
(
"finalLayerThickness",
finalLayerThickness_[patchI]
);
layerDict.readIfPresent
(
"thickness",
thickness_[patchI]
);
break;
case FINAL_AND_EXPANSION:
layerDict.readIfPresent
(
"finalLayerThickness",
finalLayerThickness_[patchI]
);
layerDict.readIfPresent
(
"expansionRatio",
expansionRatio_[patchI]
);
break;
default:
FatalIOErrorIn
(
"layerParameters::layerParameters(..)",
dict
) << "problem." << exit(FatalIOError);
break;
}
layerDict.readIfPresent
(
"minThickness",
@ -176,4 +374,190 @@ Foam::layerParameters::layerParameters
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::layerParameters::layerThickness
(
const label nLayers,
const scalar firstLayerThickess,
const scalar finalLayerThickess,
const scalar totalThickness,
const scalar expansionRatio
) const
{
switch (layerSpec_)
{
case FIRST_AND_TOTAL:
case FINAL_AND_TOTAL:
{
return totalThickness;
}
break;
case FIRST_AND_EXPANSION:
{
if (mag(expansionRatio-1) < SMALL)
{
return firstLayerThickess * nLayers;
}
else
{
return firstLayerThickess *
(1.0 - pow(expansionRatio, nLayers))
/ (1.0 - expansionRatio);
}
}
break;
case FINAL_AND_EXPANSION:
{
if (mag(expansionRatio-1) < SMALL)
{
return finalLayerThickess * nLayers;
}
else
{
scalar invExpansion = 1.0 / expansionRatio;
return finalLayerThickess *
(1.0 - pow(invExpansion, nLayers))
/ (1.0 - invExpansion);
}
}
break;
default:
{
FatalErrorIn("layerParameters::layerThickness(..)")
<< "Illegal thickness specification " << layerSpec_
<< exit(FatalError);
return -VGREAT;
}
}
}
Foam::scalar Foam::layerParameters::layerExpansionRatio
(
const label nLayers,
const scalar firstLayerThickess,
const scalar finalLayerThickess,
const scalar totalThickness,
const scalar expansionRatio
) const
{
switch (layerSpec_)
{
case FIRST_AND_EXPANSION:
case FINAL_AND_EXPANSION:
{
return expansionRatio;
}
break;
case FIRST_AND_TOTAL:
{
return layerExpansionRatio
(
nLayers,
totalThickness/firstLayerThickess
);
}
break;
case FINAL_AND_TOTAL:
{
return
1.0
/ layerExpansionRatio
(
nLayers,
totalThickness/finalLayerThickess
);
}
break;
default:
{
FatalErrorIn("layerParameters::layerThickness(..)")
<< "Illegal thickness specification" << exit(FatalError);
return -VGREAT;
}
}
}
Foam::scalar Foam::layerParameters::firstLayerThickness
(
const label nLayers,
const scalar firstLayerThickess,
const scalar finalLayerThickess,
const scalar totalThickness,
const scalar expansionRatio
) const
{
switch (layerSpec_)
{
case FIRST_AND_EXPANSION:
case FIRST_AND_TOTAL:
{
return firstLayerThickess;
}
case FINAL_AND_EXPANSION:
{
return finalLayerThickess*pow(1.0/expansionRatio, nLayers-1);
}
break;
case FINAL_AND_TOTAL:
{
scalar r = layerExpansionRatio
(
nLayers,
firstLayerThickess,
finalLayerThickess,
totalThickness,
expansionRatio
);
return finalLayerThickess/pow(r, nLayers-1);
}
break;
default:
{
FatalErrorIn("layerParameters::layerThickness(..)")
<< "Illegal thickness specification" << exit(FatalError);
return -VGREAT;
}
}
}
Foam::scalar Foam::layerParameters::finalLayerThicknessRatio
(
const label nLayers,
const scalar expansionRatio
) const
{
if (nLayers > 0)
{
if (mag(expansionRatio-1) < SMALL)
{
return 1.0/nLayers;
}
else
{
return
pow(expansionRatio, nLayers - 1)
* (1.0 - expansionRatio)
/ (1.0 - pow(expansionRatio, nLayers));
}
}
else
{
return 0.0;
}
}
// ************************************************************************* //

341
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
View File

@ -55,6 +55,27 @@ class refinementSurfaces;
class layerParameters
{
public:
// Public data types
//- Enumeration defining the layer specification:
// - first and total thickness specified
// - first and expansion ratio specified
// - final and total thickness specified
// - final and expansion ratio specified
enum layerSpecification
{
ILLEGAL,
FIRST_AND_TOTAL,
FIRST_AND_EXPANSION,
FINAL_AND_TOTAL,
FINAL_AND_EXPANSION
};
private:
// Static data members
//- Default angle for faces to be convcave
@ -68,12 +89,21 @@ class layerParameters
//- How many layers to add.
labelList numLayers_;
scalarField expansionRatio_;
//- Are sizes relative to local cell size
Switch relativeSizes_;
scalarField finalLayerThickness_;
//- How thickness is specified.
layerSpecification layerSpec_;
scalarField firstLayerThickness_;
scalarField finalLayerThickness_;
scalarField thickness_;
scalarField expansionRatio_;
//- Minimum total thickness
scalarField minThickness_;
@ -111,6 +141,14 @@ class layerParameters
// Private Member Functions
//- Calculate expansion ratio from overall size v.s. thickness of
// first layer.
scalar layerExpansionRatio
(
const label n,
const scalar totalOverFirst
) const;
//- Disallow default bitwise copy construct
layerParameters(const layerParameters&);
@ -128,154 +166,209 @@ public:
// Member Functions
// Access
// Per patch information
// Per patch information
//- How many layers to add.
// -1 : no specification. Assume 0 layers but allow sliding
// to make layers
// 0 : specified to have 0 layers. No sliding allowed.
// >0 : number of layers
const labelList& numLayers() const
{
return numLayers_;
}
// Expansion factor for layer mesh
const scalarField& expansionRatio() const
{
return expansionRatio_;
}
//- Are size parameters relative to inner cell size or
// absolute distances.
bool relativeSizes() const
{
return relativeSizes_;
}
//- Wanted thickness of final added cell layer. If multiple
// layers is the thickness of the layer furthest away
// from the wall (i.e. nearest the original mesh)
// If relativeSize() this number is relative to undistorted
// size of the cell outside layer.
const scalarField& finalLayerThickness() const
{
return finalLayerThickness_;
}
//- Minimum thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// If relativeSize() this number is relative to undistorted
// size of the cell outside layer.
const scalarField& minThickness() const
{
return minThickness_;
}
scalar featureAngle() const
//- How many layers to add.
// -1 : no specification. Assume 0 layers but allow sliding
// to make layers
// 0 : specified to have 0 layers. No sliding allowed.
// >0 : number of layers
const labelList& numLayers() const
{
return featureAngle_;
return numLayers_;
}
//- At non-patched sides allow mesh to slip if extrusion
// direction makes angle larger than slipFeatureAngle.
scalar slipFeatureAngle() const
//- Are size parameters relative to inner cell size or
// absolute distances.
bool relativeSizes() const
{
return slipFeatureAngle_;
return relativeSizes_;
}
scalar concaveAngle() const
// Expansion factor for layer mesh
const scalarField& expansionRatio() const
{
return concaveAngle_;
return expansionRatio_;
}
//- If points get not extruded do nGrow layers of connected faces
// that are not grown. Is used to not do layers at all close to
// features.
label nGrow() const
//- Wanted thickness of the layer furthest away
// from the wall (i.e. nearest the original mesh).
// If relativeSize() this number is relative to undistorted
// size of the cell outside layer.
const scalarField& finalLayerThickness() const
{
return nGrow_;
return finalLayerThickness_;
}
//- Number of smoothing iterations of surface normals
label nSmoothSurfaceNormals() const
//- Wanted thickness of the layer nearest to the wall.
// If relativeSize() this number is relative to undistorted
// size of the cell outside layer.
const scalarField& firstLayerThickness() const
{
return nSmoothSurfaceNormals_;
return firstLayerThickness_;
}
//- Number of smoothing iterations of interior mesh movement
// direction
label nSmoothNormals() const
//- Wanted overall thickness of all layers.
// If relativeSize() this number is relative to undistorted
// size of the cell outside layer.
const scalarField& thickness() const
{
return nSmoothNormals_;
return thickness_;
}
//- Stop layer growth on highly warped cells
scalar maxFaceThicknessRatio() const
//- Minimum overall thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// If relativeSize() this number is relative to undistorted
// size of the cell outside layer.
const scalarField& minThickness() const
{
return maxFaceThicknessRatio_;
}
scalar layerTerminationCos() const
{
return layerTerminationCos_;
}
//- Smooth layer thickness over surface patches
label nSmoothThickness() const
{
return nSmoothThickness_;
}
//- Reduce layer growth where ratio thickness to medial
// distance is large
scalar maxThicknessToMedialRatio() const
{
return maxThicknessToMedialRatio_;
}
//- Angle used to pick up medial axis points
scalar minMedianAxisAngleCos() const
{
return minMedianAxisAngleCos_;
}
//- Create buffer region for new layer terminations
label nBufferCellsNoExtrude() const
{
return nBufferCellsNoExtrude_;
}
label nSnap() const
{
return nSnap_;
}
const Switch& additionalReporting() const
{
return additionalReporting_;
return minThickness_;
}
// Overall
scalar featureAngle() const
{
return featureAngle_;
}
//- Number of overall layer addition iterations
label nLayerIter() const
{
return nLayerIter_;
}
//- At non-patched sides allow mesh to slip if extrusion
// direction makes angle larger than slipFeatureAngle.
scalar slipFeatureAngle() const
{
return slipFeatureAngle_;
}
//- Number of iterations after which relaxed motion rules
// are to be used.
label nRelaxedIter() const
{
return nRelaxedIter_;
}
scalar concaveAngle() const
{
return concaveAngle_;
}
//- If points get not extruded do nGrow layers of connected faces
// that are not grown. Is used to not do layers at all close to
// features.
label nGrow() const
{
return nGrow_;
}
//- Number of smoothing iterations of surface normals
label nSmoothSurfaceNormals() const
{
return nSmoothSurfaceNormals_;
}
//- Number of smoothing iterations of interior mesh movement
// direction
label nSmoothNormals() const
{
return nSmoothNormals_;
}
//- Stop layer growth on highly warped cells
scalar maxFaceThicknessRatio() const
{
return maxFaceThicknessRatio_;
}
scalar layerTerminationCos() const
{
return layerTerminationCos_;
}
//- Smooth layer thickness over surface patches
label nSmoothThickness() const
{
return nSmoothThickness_;
}
//- Reduce layer growth where ratio thickness to medial
// distance is large
scalar maxThicknessToMedialRatio() const
{
return maxThicknessToMedialRatio_;
}
//- Angle used to pick up medial axis points
scalar minMedianAxisAngleCos() const
{
return minMedianAxisAngleCos_;
}
//- Create buffer region for new layer terminations
label nBufferCellsNoExtrude() const
{
return nBufferCellsNoExtrude_;
}
label nSnap() const
{
return nSnap_;
}
const Switch& additionalReporting() const
{
return additionalReporting_;
}
// Overall
//- Number of overall layer addition iterations
label nLayerIter() const
{
return nLayerIter_;
}
//- Number of iterations after which relaxed motion rules
// are to be used.
label nRelaxedIter() const
{
return nRelaxedIter_;
}
// Helper
//- Determine overall thickness. Uses two of the four parameters
// according to the layerSpecification
scalar layerThickness
(
const label nLayers,
const scalar firstLayerThickess,
const scalar finalLayerThickess,
const scalar totalThickness,
const scalar expansionRatio
) const;
//- Determine expansion ratio. Uses two of the four parameters
// according to the layerSpecification
scalar layerExpansionRatio
(
const label nLayers,
const scalar firstLayerThickess,
const scalar finalLayerThickess,
const scalar totalThickness,
const scalar expansionRatio
) const;
//- Determine first layer (near-wall) thickness. Uses two of the
// four parameters according to the layerSpecification
scalar firstLayerThickness
(
const label nLayers,
const scalar firstLayerThickess,
const scalar finalLayerThickess,
const scalar totalThickness,
const scalar expansionRatio
) const;
//- Determine ratio of final layer thickness to
// overall layer thickness
scalar finalLayerThicknessRatio
(
const label nLayers,
const scalar expansionRatio
) const;
};

6
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
View File

@ -36,7 +36,11 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
nSnap_(readLabel(dict.lookup("nRelaxIter"))),
nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false))
implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false)),
multiRegionFeatureSnap_
(
dict.lookupOrDefault("multiRegionFeatureSnap", false)
)
{}

7
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
View File

@ -68,6 +68,8 @@ class snapParameters
const Switch implicitFeatureSnap_;
const Switch multiRegionFeatureSnap_;
// Private Member Functions
@ -134,6 +136,11 @@ public:
return implicitFeatureSnap_;
}
Switch multiRegionFeatureSnap() const
{
return multiRegionFeatureSnap_;
}
};

9
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -264,6 +264,14 @@ private:
label& nRefine
) const;
//- Mark cells for distance-to-feature based refinement.
label markInternalDistanceToFeatureRefinement
(
const label nAllowRefine,
labelList& refineCell,
label& nRefine
) const;
//- Mark cells for refinement-shells based refinement.
label markInternalRefinement
(
@ -656,6 +664,7 @@ public:
const scalar curvature,
const bool featureRefinement,
const bool featureDistanceRefinement,
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,

119
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -290,7 +290,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
forAll(features_, featI)
{
const featureEdgeMesh& featureMesh = features_[featI];
const label featureLevel = features_.levels()[featI];
const label featureLevel = features_.levels()[featI][0];
const labelListList& pointEdges = featureMesh.pointEdges();
// Find regions on edgeMesh
@ -511,6 +511,90 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
}
// Mark cells for distance-to-feature based refinement.
Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
(
const label nAllowRefine,
labelList& refineCell,
label& nRefine
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
const pointField& cellCentres = mesh_.cellCentres();
// Detect if there are any distance shells
if (features_.maxDistance() <= 0.0)
{
return 0;
}
label oldNRefine = nRefine;
// Collect cells to test
pointField testCc(cellLevel.size()-nRefine);
labelList testLevels(cellLevel.size()-nRefine);
label testI = 0;
forAll(cellLevel, cellI)
{
if (refineCell[cellI] == -1)
{
testCc[testI] = cellCentres[cellI];
testLevels[testI] = cellLevel[cellI];
testI++;
}
}
// Do test to see whether cells is inside/outside shell with higher level
labelList maxLevel;
features_.findHigherLevel(testCc, testLevels, maxLevel);
// Mark for refinement. Note that we didn't store the original cellID so
// now just reloop in same order.
testI = 0;
forAll(cellLevel, cellI)
{
if (refineCell[cellI] == -1)
{
if (maxLevel[testI] > testLevels[testI])
{
bool reachedLimit = !markForRefine
(
maxLevel[testI], // mark with any positive value
nAllowRefine,
refineCell[cellI],
nRefine
);
if (reachedLimit)
{
if (debug)
{
Pout<< "Stopped refining internal cells"
<< " since reaching my cell limit of "
<< mesh_.nCells()+7*nRefine << endl;
}
break;
}
}
testI++;
}
}
if
(
returnReduce(nRefine, sumOp<label>())
> returnReduce(nAllowRefine, sumOp<label>())
)
{
Info<< "Reached refinement limit." << endl;
}
return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
// Mark cells for non-surface intersection based refinement.
Foam::label Foam::meshRefinement::markInternalRefinement
(
@ -1128,6 +1212,7 @@ Foam::labelList Foam::meshRefinement::refineCandidates
const scalar curvature,
const bool featureRefinement,
const bool featureDistanceRefinement,
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,
@ -1196,8 +1281,24 @@ Foam::labelList Foam::meshRefinement::refineCandidates
nRefine
);
Info<< "Marked for refinement due to explicit features : "
<< nFeatures << " cells." << endl;
Info<< "Marked for refinement due to explicit features "
<< ": " << nFeatures << " cells." << endl;
}
// Inside distance-to-feature shells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (featureDistanceRefinement)
{
label nShell = markInternalDistanceToFeatureRefinement
(
nAllowRefine,
refineCell,
nRefine
);
Info<< "Marked for refinement due to distance to explicit features "
": " << nShell << " cells." << endl;
}
// Inside refinement shells
@ -1212,8 +1313,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
refineCell,
nRefine
);
Info<< "Marked for refinement due to refinement shells : "
<< nShell << " cells." << endl;
Info<< "Marked for refinement due to refinement shells "
<< ": " << nShell << " cells." << endl;
}
// Refinement based on intersection of surface
@ -1230,8 +1331,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
refineCell,
nRefine
);
Info<< "Marked for refinement due to surface intersection : "
<< nSurf << " cells." << endl;
Info<< "Marked for refinement due to surface intersection "
<< ": " << nSurf << " cells." << endl;
}
// Refinement based on curvature of surface
@ -1254,8 +1355,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
refineCell,
nRefine
);
Info<< "Marked for refinement due to curvature/regions : "
<< nCurv << " cells." << endl;
Info<< "Marked for refinement due to curvature/regions "
<< ": " << nCurv << " cells." << endl;
}
// Pack cells-to-refine

180
src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C
View File

@ -25,6 +25,7 @@ License
#include "refinementFeatures.H"
#include "Time.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -34,15 +35,15 @@ void Foam::refinementFeatures::read
const PtrList<dictionary>& featDicts
)
{
forAll(featDicts, i)
forAll(featDicts, featI)
{
const dictionary& dict = featDicts[i];
const dictionary& dict = featDicts[featI];
fileName featFileName(dict.lookup("file"));
set
(
i,
featI,
new featureEdgeMesh
(
IOobject
@ -58,15 +59,74 @@ void Foam::refinementFeatures::read
)
);
const featureEdgeMesh& eMesh = operator[](i);
const featureEdgeMesh& eMesh = operator[](featI);
//eMesh.mergePoints(meshRefiner_.mergeDistance());
levels_[i] = readLabel(dict.lookup("level"));
Info<< "Refinement level " << levels_[i]
<< " for all cells crossed by feature " << featFileName
<< " (" << eMesh.points().size() << " points, "
if (dict.found("levels"))
{
List<Tuple2<scalar, label> > distLevels(dict["levels"]);
if (dict.size() < 1)
{
FatalErrorIn
(
"refinementFeatures::read"
"(const objectRegistry&"
", const PtrList<dictionary>&)"
) << " : levels should be at least size 1" << endl
<< "levels : " << dict["levels"]
<< exit(FatalError);
}
distances_[featI].setSize(distLevels.size());
levels_[featI].setSize(distLevels.size());
forAll(distLevels, j)
{
distances_[featI][j] = distLevels[j].first();
levels_[featI][j] = distLevels[j].second();
// Check in incremental order
if (j > 0)
{
if
(
(distances_[featI][j] <= distances_[featI][j-1])
|| (levels_[featI][j] > levels_[featI][j-1])
)
{
FatalErrorIn
(
"refinementFeatures::read"
"(const objectRegistry&"
", const PtrList<dictionary>&)"
) << " : Refinement should be specified in order"
<< " of increasing distance"
<< " (and decreasing refinement level)." << endl
<< "Distance:" << distances_[featI][j]
<< " refinementLevel:" << levels_[featI][j]
<< exit(FatalError);
}
}
}
}
else
{
// Look up 'level' for single level
levels_[featI] = labelList(1, readLabel(dict.lookup("level")));
distances_[featI] = scalarField(1, 0.0);
}
Info<< "Refinement level according to distance to "
<< featFileName << " (" << eMesh.points().size() << " points, "
<< eMesh.edges().size() << " edges)." << endl;
forAll(levels_[featI], j)
{
Info<< " level " << levels_[featI][j]
<< " for all cells within " << distances_[featI][j]
<< " meter." << endl;
}
}
}
@ -127,6 +187,80 @@ void Foam::refinementFeatures::buildTrees
}
// Find maximum level of a shell.
void Foam::refinementFeatures::findHigherLevel
(
const pointField& pt,
const label featI,
labelList& maxLevel
) const
{
const labelList& levels = levels_[featI];
const scalarField& distances = distances_[featI];
// Collect all those points that have a current maxLevel less than
// (any of) the shell. Also collect the furthest distance allowable
// to any shell with a higher level.
pointField candidates(pt.size());
labelList candidateMap(pt.size());
scalarField candidateDistSqr(pt.size());
label candidateI = 0;
forAll(maxLevel, pointI)
{
forAllReverse(levels, levelI)
{
if (levels[levelI] > maxLevel[pointI])
{
candidates[candidateI] = pt[pointI];
candidateMap[candidateI] = pointI;
candidateDistSqr[candidateI] = sqr(distances[levelI]);
candidateI++;
break;
}
}
}
candidates.setSize(candidateI);
candidateMap.setSize(candidateI);
candidateDistSqr.setSize(candidateI);
// Do the expensive nearest test only for the candidate points.
const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
List<pointIndexHit> nearInfo(candidates.size());
forAll(candidates, candidateI)
{
nearInfo[candidateI] = tree.findNearest
(
candidates[candidateI],
candidateDistSqr[candidateI]
);
}
// Update maxLevel
forAll(nearInfo, candidateI)
{
if (nearInfo[candidateI].hit())
{
// Check which level it actually is in.
label minDistI = findLower
(
distances,
mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
);
label pointI = candidateMap[candidateI];
// pt is inbetween shell[minDistI] and shell[minDistI+1]
maxLevel[pointI] = levels[minDistI+1];
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::refinementFeatures::refinementFeatures
@ -136,6 +270,7 @@ Foam::refinementFeatures::refinementFeatures
)
:
PtrList<featureEdgeMesh>(featDicts.size()),
distances_(featDicts.size()),
levels_(featDicts.size()),
edgeTrees_(featDicts.size()),
pointTrees_(featDicts.size())
@ -175,6 +310,7 @@ Foam::refinementFeatures::refinementFeatures
)
:
PtrList<featureEdgeMesh>(featDicts.size()),
distances_(featDicts.size()),
levels_(featDicts.size()),
edgeTrees_(featDicts.size()),
pointTrees_(featDicts.size())
@ -336,4 +472,32 @@ void Foam::refinementFeatures::findNearestPoint
}
void Foam::refinementFeatures::findHigherLevel
(
const pointField& pt,
const labelList& ptLevel,
labelList& maxLevel
) const
{
// Maximum level of any shell. Start off with level of point.
maxLevel = ptLevel;
forAll(*this, featI)
{
findHigherLevel(pt, featI, maxLevel);
}
}
Foam::scalar Foam::refinementFeatures::maxDistance() const
{
scalar overallMax = -GREAT;
forAll(distances_, featI)
{
overallMax = max(overallMax, max(distances_[featI]));
}
return overallMax;
}
// ************************************************************************* //

34
src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H
View File

@ -57,8 +57,11 @@ private:
// Private data
//- Refinement levels
labelList levels_;
//- Per shell the list of ranges
List<scalarField> distances_;
//- Per shell per distance the refinement level
labelListList levels_;
//- Edge
PtrList<indexedOctree<treeDataEdge> > edgeTrees_;
@ -75,6 +78,13 @@ private:
//- Build edge tree and feature point tree
void buildTrees(const label, const labelList&);
//- Find shell level higher than ptLevel
void findHigherLevel
(
const pointField& pt,
const label featI,
labelList& maxLevel
) const;
public:
@ -101,11 +111,18 @@ public:
// Access
const labelList& levels() const
//- Per featureEdgeMesh the list of level
const labelListList& levels() const
{
return levels_;
}
//- Per featureEdgeMesh the list of ranges
const List<scalarField>& distances() const
{
return distances_;
}
const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
{
return edgeTrees_;
@ -119,6 +136,9 @@ public:
// Query
//- Highest distance of all features
scalar maxDistance() const;
//- Find nearest point on nearest feature edge
void findNearestEdge
(
@ -141,6 +161,14 @@ public:
labelList& nearIndex
) const;
//- Find shell level higher than ptLevel
void findHigherLevel
(
const pointField& pt,
const labelList& ptLevel,
labelList& maxLevel
) const;
};

7
src/regionModels/pyrolysisModels/Make/options
View File

@ -3,10 +3,10 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
@ -20,12 +20,13 @@ LIB_LIBS = \
-lmeshTools \
-lchemistryModel \
-lspecie \
-lsolidSpecie \
-lfluidThermophysicalModels \
-lsolidChemistryModel \
-lsolidThermo \
-lcompressibleTurbulenceModel \
-lcompressibleRASModels \
-lcompressibleLESModels \
-lLESdeltas \
-lregionModels \
-lradiationModels
-lradiationModels \
-lreactionThermophysicalModels

2
src/regionModels/thermoBaffleModels/Make/options
View File

@ -1,7 +1,6 @@
EXE_INC = \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solid/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
@ -18,5 +17,4 @@ LIB_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lOpenFOAM \
-lsolidSpecie \
-lradiationModels

1
src/sampling/Make/files
View File

@ -18,6 +18,7 @@ sampledSet/sampledSets/sampledSetsGrouping.C
sampledSet/sampledSetsFunctionObject/sampledSetsFunctionObject.C
sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
sampledSet/uniform/uniformSet.C
sampledSet/array/arraySet.C
setWriters = sampledSet/writers

195
src/sampling/sampledSet/array/arraySet.C Normal file
View File

@ -0,0 +1,195 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "arraySet.H"
#include "sampledSet.H"
#include "meshSearch.H"
#include "DynamicList.H"
#include "polyMesh.H"
#include "addToRunTimeSelectionTable.H"
#include "word.H"
#include "transform.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(arraySet, 0);
addToRunTimeSelectionTable(sampledSet, arraySet, word);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::arraySet::calcSamples
(
DynamicList<point>& samplingPts,
DynamicList<label>& samplingCells,
DynamicList<label>& samplingFaces,
DynamicList<label>& samplingSegments,
DynamicList<scalar>& samplingCurveDist
) const
{
const meshSearch& queryMesh = searchEngine();
label nTotalSamples
(
pointsDensity_.x()
*pointsDensity_.y()
*pointsDensity_.z()
);
List<point> sampleCoords(nTotalSamples);
const scalar deltax = spanBox_.x()/(pointsDensity_.x() + 1);
const scalar deltay = spanBox_.y()/(pointsDensity_.y() + 1);
const scalar deltaz = spanBox_.z()/(pointsDensity_.z() + 1);
label p(0);
for (label k=1; k<=pointsDensity_.z(); k++)
{
for (label j=1; j<=pointsDensity_.y(); j++)
{
for (label i=1; i<=pointsDensity_.x(); i++)
{
vector t(deltax*i , deltay*j, deltaz*k);
sampleCoords[p] = coordSys_.origin() + t;
p++;
}
}
}
forAll(sampleCoords, i)
{
sampleCoords[i] = transform(coordSys_.R(), sampleCoords[i]);
}
forAll(sampleCoords, sampleI)
{
label cellI = queryMesh.findCell(sampleCoords[sampleI]);
if (cellI != -1)
{
samplingPts.append(sampleCoords[sampleI]);
samplingCells.append(cellI);
samplingFaces.append(-1);
samplingSegments.append(0);
samplingCurveDist.append(1.0 * sampleI);
}
}
}
void Foam::arraySet::genSamples()
{
// Storage for sample points
DynamicList<point> samplingPts;
DynamicList<label> samplingCells;
DynamicList<label> samplingFaces;
DynamicList<label> samplingSegments;
DynamicList<scalar> samplingCurveDist;
calcSamples
(
samplingPts,
samplingCells,
samplingFaces,
samplingSegments,
samplingCurveDist
);
samplingPts.shrink();
samplingCells.shrink();
samplingFaces.shrink();
samplingSegments.shrink();
samplingCurveDist.shrink();
setSamples
(
samplingPts,
samplingCells,
samplingFaces,
samplingSegments,
samplingCurveDist
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::arraySet::arraySet
(
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const word& axis,
const coordinateSystem& origin,
const Vector<label>& pointsDensity,
const Vector<scalar>& spanBox
)
:
sampledSet(name, mesh, searchEngine, axis),
coordSys_(origin),
pointsDensity_(pointsDensity),
spanBox_(spanBox)
{
genSamples();
if (debug)
{
write(Info);
}
}
Foam::arraySet::arraySet
(
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const dictionary& dict
)
:
sampledSet(name, mesh, searchEngine, dict),
coordSys_(dict),
pointsDensity_(dict.lookup("pointsDensity")),
spanBox_(dict.lookup("spanBox"))
{
genSamples();
if (debug)
{
write(Info);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::arraySet::~arraySet()
{}
// ************************************************************************* //

108
src/thermophysicalModels/solidSpecie/reaction/Reactions/IrreversibleSolidReaction/IrreversibleSolidReaction.H → src/sampling/sampledSet/array/arraySet.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,103 +22,99 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::IrreversibleSolidReaction
Foam::arraySet
Description
Simple extension of Reaction to handle reversible reactions
SourceFiles
IrreversibleSolidReaction.C
arraySet.C
\*---------------------------------------------------------------------------*/
#ifndef IrreversibleSolidReaction_H
#define IrreversibleSolidReaction_H
#ifndef arraySet_H
#define arraySet_H
#include "solidReaction.H"
#include "sampledSet.H"
#include "DynamicList.H"
#include "coordinateSystem.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class passiveParticle;
template<class Type> class particle;
/*---------------------------------------------------------------------------*\
Class IrreversibleSolidReaction Declaration
Class arraySet Declaration
\*---------------------------------------------------------------------------*/
template<class ReactionRate>
class IrreversibleSolidReaction
class arraySet
:
public solidReaction
public sampledSet
{
// Private data
// Reaction rate
ReactionRate k_;
//- Coordinate syste
coordinateSystem coordSys_;
// Reaction order
scalar nReact_;
//- Point density vector
Vector<label> pointsDensity_;
//- Span box
Vector<scalar> spanBox_;
// Private Member Functions
//- Disallow default bitwise assignment
void operator=
//- Samples all points in sampleCoords.
void calcSamples
(
const IrreversibleSolidReaction<ReactionRate>&
);
DynamicList<point>& samplingPts,
DynamicList<label>& samplingCells,
DynamicList<label>& samplingFaces,
DynamicList<label>& samplingSegments,
DynamicList<scalar>& samplingCurveDist
) const;
//- Uses calcSamples to obtain samples. Copies them into *this.
void genSamples();
public:
//- Runtime type information
TypeName("irreversible");
TypeName("array");
// Constructors
//- Construct from components
IrreversibleSolidReaction
arraySet
(
const solidReaction& reaction,
const ReactionRate& reactionRate,
const scalar nReact
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const word& axis,
const coordinateSystem& coordSys,
const Vector<label>& pointsDensity,
const Vector<scalar>& spanBox
);
//- Construct from Istream
IrreversibleSolidReaction
//- Construct from dictionary
arraySet
(
const speciesTable& components,
Istream& is,
const speciesTable& pyrolysisGases
const word& name,
const polyMesh& mesh,
const meshSearch& searchEngine,
const dictionary& dict
);
//- Destructor
virtual ~IrreversibleSolidReaction()
{}
// Member Functions
// IrreversibleSolidReaction rate coefficients
//- Forward rate constant
virtual scalar kf
(
const scalar p,
const scalar T,
const scalarField& c
) const;
//- Reaction order
virtual scalar nReact() const;
//- Write
virtual void write(Ostream&) const;
virtual ~arraySet();
};
@ -128,12 +124,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "IrreversibleSolidReaction.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

22
src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
View File

@ -164,17 +164,19 @@ void Foam::sampledSurfaces::write()
writeGeometry();
}
sampleAndWrite<volScalarField>(mesh_);
sampleAndWrite<volVectorField>(mesh_);
sampleAndWrite<volSphericalTensorField>(mesh_);
sampleAndWrite<volSymmTensorField>(mesh_);
sampleAndWrite<volTensorField>(mesh_);
const IOobjectList objects(mesh_, mesh_.time().timeName());
sampleAndWrite<surfaceScalarField>(mesh_);
sampleAndWrite<surfaceVectorField>(mesh_);
sampleAndWrite<surfaceSphericalTensorField>(mesh_);
sampleAndWrite<surfaceSymmTensorField>(mesh_);
sampleAndWrite<surfaceTensorField>(mesh_);
sampleAndWrite<volScalarField>(objects);
sampleAndWrite<volVectorField>(objects);
sampleAndWrite<volSphericalTensorField>(objects);
sampleAndWrite<volSymmTensorField>(objects);
sampleAndWrite<volTensorField>(objects);
sampleAndWrite<surfaceScalarField>(objects);
sampleAndWrite<surfaceVectorField>(objects);
sampleAndWrite<surfaceSphericalTensorField>(objects);
sampleAndWrite<surfaceSymmTensorField>(objects);
sampleAndWrite<surfaceTensorField>(objects);
}
}

2
src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
View File

@ -160,7 +160,7 @@ class sampledSurfaces
//- Sample and write all sampled fields
template<class Type>
void sampleAndWrite(const fvMesh&);
void sampleAndWrite(const IOobjectList&);
//- Disallow default bitwise copy construct and assignment
sampledSurfaces(const sampledSurfaces&);

25
src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C
View File

@ -166,13 +166,18 @@ void Foam::sampledSurfaces::sampleAndWrite
template<class GeoField>
void Foam::sampledSurfaces::sampleAndWrite(const fvMesh& mesh)
void Foam::sampledSurfaces::sampleAndWrite(const IOobjectList& allObjects)
{
forAll (fieldSelection_, fieldI)
{
const wordRe field = fieldSelection_[fieldI];
IOobject* fieldIOPtr = allObjects.lookup(field);
if (mesh.thisDb().foundObject<GeoField>(field))
if
(
fieldIOPtr != NULL
&& fieldIOPtr->headerClassName() == GeoField::typeName
)
{
if (Pstream::master() && verbose_)
{
@ -181,17 +186,25 @@ void Foam::sampledSurfaces::sampleAndWrite(const fvMesh& mesh)
if (loadFromFiles_)
{
const GeoField& geoField =
mesh.thisDb().lookupObject<GeoField>(field);
const GeoField geoField
(
IOobject
(
field,
mesh_.time().timeName(),
mesh_,
IOobject::MUST_READ
),
mesh_
);
const_cast<GeoField&>(geoField).readOpt() = IOobject::MUST_READ;
sampleAndWrite(geoField);
}
else
{
sampleAndWrite
(
mesh.thisDb().lookupObject<GeoField>(field)
mesh_.thisDb().lookupObject<GeoField>(field)
);
}
}

3
src/thermophysicalModels/Allwmake
View File

@ -4,6 +4,7 @@ makeType=${1:-libso}
set -x
wmake $makeType specie
wmake $makeType solidSpecie
wmake $makeType thermophysicalFunctions
./properties/Allwmake $*
@ -14,7 +15,7 @@ wmake $makeType chemistryModel
wmake $makeType barotropicCompressibilityModel
wmake $makeType SLGThermo
wmake $makeType solidSpecie
wmake $makeType solidThermo
wmake $makeType solidChemistryModel

21
src/thermophysicalModels/basic/mixtures/pureMixture/pureMixture.H
View File

@ -87,6 +87,27 @@ public:
return mixture_;
}
const ThermoType& cellVolMixture
(
const scalar,
const scalar,
const label
) const
{
return mixture_;
}
const ThermoType& patchFaceVolMixture
(
const scalar,
const scalar,
const label,
const label
) const
{
return mixture_;
}
//- Read dictionary
void read(const dictionary&);
};

8
src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
View File

@ -27,7 +27,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "unitConversion.H"
#include "zeroGradientFvPatchFields.H"
#include "basicSolidMixture.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -71,7 +71,7 @@ greyMeanSolidAbsorptionEmission::X(const word specie) const
}
}
const scalarField& Yj = mixture_.Y(specie);
const label mySpecieI = mixture_.components()[specie];
const label mySpecieI = mixture_.species()[specie];
forAll(Xj, iCell)
{
Xj[iCell] = Yj[iCell]/mixture_.rho(mySpecieI, p[iCell], T[iCell]);
@ -93,10 +93,10 @@ greyMeanSolidAbsorptionEmission
coeffsDict_((dict.subDict(typeName + "Coeffs"))),
thermo_(mesh.lookupObject<solidThermo>("thermophysicalProperties")),
speciesNames_(0),
mixture_(dynamic_cast<const basicSolidMixture&>(thermo_)),
mixture_(dynamic_cast<const basicMultiComponentMixture&>(thermo_)),
solidData_(mixture_.Y().size())
{
if (!isA<basicSolidMixture>(thermo_))
if (!isA<basicMultiComponentMixture>(thermo_))
{
FatalErrorIn
(

2
src/thermophysicalModels/radiationModels/submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.H
View File

@ -90,7 +90,7 @@ private:
HashTable<label> speciesNames_;
//- Basic multicomponent mixture
const basicSolidMixture& mixture_;
const basicMultiComponentMixture& mixture_;
//- List of solid species data
List<FixedList<scalar, 2> > solidData_;

3
src/thermophysicalModels/reactionThermo/Make/options
View File

@ -2,10 +2,13 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
LIB_LIBS = \
-lfiniteVolume \
-lfluidThermophysicalModels \
-lspecie \
-lsolidSpecie \
-lmeshTools

5
src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/makeChemistryReaders.C
View File

@ -25,6 +25,7 @@ License
#include "makeReactionThermo.H"
#include "thermoPhysicsTypes.H"
#include "solidThermoPhysicsTypes.H"
#include "chemistryReader.H"
#include "foamChemistryReader.H"
@ -41,6 +42,8 @@ makeChemistryReader(gasThermoPhysics);
makeChemistryReader(constIncompressibleGasThermoPhysics);
makeChemistryReader(incompressibleGasThermoPhysics);
makeChemistryReader(icoPoly8ThermoPhysics);
makeChemistryReader(hConstSolidThermoPhysics);
makeChemistryReader(hExponentialSolidThermoPhysics);
makeChemistryReaderType(foamChemistryReader, constGasThermoPhysics);
makeChemistryReaderType(foamChemistryReader, gasThermoPhysics);
@ -51,6 +54,8 @@ makeChemistryReaderType
);
makeChemistryReaderType(foamChemistryReader, incompressibleGasThermoPhysics);
makeChemistryReaderType(foamChemistryReader, icoPoly8ThermoPhysics);
makeChemistryReaderType(foamChemistryReader, hConstSolidThermoPhysics);
makeChemistryReaderType(foamChemistryReader, hExponentialSolidThermoPhysics);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

17
src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
View File

@ -41,6 +41,7 @@ License
#include "powerSeriesReactionRate.H"
#include "addToRunTimeSelectionTable.H"
/* * * * * * * * * * * * * * * * * Static data * * * * * * * * * * * * * * * */
namespace Foam
@ -178,7 +179,8 @@ void Foam::chemkinReader::addReactionType
{
reactions_.append
(
new IrreversibleReaction<gasThermoPhysics, ReactionRateType>
new IrreversibleReaction
<Reaction, gasThermoPhysics, ReactionRateType>
(
Reaction<gasThermoPhysics>
(
@ -197,7 +199,8 @@ void Foam::chemkinReader::addReactionType
{
reactions_.append
(
new ReversibleReaction<gasThermoPhysics, ReactionRateType>
new ReversibleReaction
<Reaction, gasThermoPhysics, ReactionRateType>
(
Reaction<gasThermoPhysics>
(
@ -496,7 +499,7 @@ void Foam::chemkinReader::addReaction
reactions_.append
(
new NonEquilibriumReversibleReaction
<gasThermoPhysics, ArrheniusReactionRate>
<Reaction, gasThermoPhysics, ArrheniusReactionRate>
(
Reaction<gasThermoPhysics>
(
@ -549,7 +552,11 @@ void Foam::chemkinReader::addReaction
reactions_.append
(
new NonEquilibriumReversibleReaction
<gasThermoPhysics, thirdBodyArrheniusReactionRate>
<
Reaction,
gasThermoPhysics,
thirdBodyArrheniusReactionRate
>
(
Reaction<gasThermoPhysics>
(
@ -654,7 +661,7 @@ void Foam::chemkinReader::addReaction
reactions_.append
(
new NonEquilibriumReversibleReaction
<gasThermoPhysics, LandauTellerReactionRate>
<Reaction, gasThermoPhysics, LandauTellerReactionRate>
(
Reaction<gasThermoPhysics>
(

65
src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
View File

@ -77,7 +77,8 @@ Foam::multiComponentMixture<ThermoType>::multiComponentMixture
:
basicMultiComponentMixture(thermoDict, specieNames, mesh),
speciesData_(species_.size()),
mixture_("mixture", *thermoData[specieNames[0]])
mixture_("mixture", *thermoData[specieNames[0]]),
mixtureVol_("volMixture", *thermoData[specieNames[0]])
{
forAll(species_, i)
{
@ -101,7 +102,8 @@ Foam::multiComponentMixture<ThermoType>::multiComponentMixture
:
basicMultiComponentMixture(thermoDict, thermoDict.lookup("species"), mesh),
speciesData_(species_.size()),
mixture_("mixture", constructSpeciesData(thermoDict))
mixture_("mixture", constructSpeciesData(thermoDict)),
mixtureVol_("volMixture", speciesData_[0])
{
correctMassFractions();
}
@ -148,6 +150,65 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::patchFaceMixture
}
template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellVolMixture
(
const scalar p,
const scalar T,
const label celli
) const
{
scalar rhoInv = 0.0;
forAll(speciesData_, i)
{
rhoInv += Y_[i][celli]/speciesData_[i].rho(p, T);
}
mixtureVol_ =
Y_[0][celli]/speciesData_[0].rho(p, T)/rhoInv*speciesData_[0];
for (label n=1; n<Y_.size(); n++)
{
mixtureVol_ +=
Y_[n][celli]/speciesData_[n].rho(p, T)/rhoInv*speciesData_[n];
}
return mixtureVol_;
}
template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::
patchFaceVolMixture
(
const scalar p,
const scalar T,
const label patchi,
const label facei
) const
{
scalar rhoInv = 0.0;
forAll(speciesData_, i)
{
rhoInv +=
Y_[i].boundaryField()[patchi][facei]/speciesData_[i].rho(p, T);
}
mixtureVol_ =
Y_[0].boundaryField()[patchi][facei]/speciesData_[0].rho(p, T)/rhoInv
* speciesData_[0];
for (label n=1; n<Y_.size(); n++)
{
mixtureVol_ +=
Y_[n].boundaryField()[patchi][facei]/speciesData_[n].rho(p,T)
/ rhoInv*speciesData_[n];
}
return mixtureVol_;
}
template<class ThermoType>
void Foam::multiComponentMixture<ThermoType>::read
(

19
src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.H
View File

@ -60,6 +60,10 @@ class multiComponentMixture
//- Temporary storage for the cell/face mixture thermo data
mutable ThermoType mixture_;
//- Temporary storage for the volume weighted
// cell/face mixture thermo data
mutable ThermoType mixtureVol_;
// Private Member Functions
@ -110,6 +114,21 @@ public:
const label facei
) const;
const ThermoType& cellVolMixture
(
const scalar p,
const scalar T,
const label celli
) const;
const ThermoType& patchFaceVolMixture
(
const scalar p,
const scalar T,
const label patchi,
const label facei
) const;
//- Return the raw specie thermodynamic data
const PtrList<ThermoType>& speciesData() const
{

6
src/thermophysicalModels/solidChemistryModel/Make/options
View File

@ -1,10 +1,11 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude
@ -12,4 +13,5 @@ EXE_INC = \
LIB_LIBS = \
-lchemistryModel \
-lfiniteVolume \
-lODE
-lODE\
-lreactionThermophysicalModels

70
src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.C
View File

@ -24,7 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "ODESolidChemistryModel.H"
#include "reactingSolidMixture.H"
#include "reactingMixture.H"
#include "solidReaction.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -38,24 +39,20 @@ ODESolidChemistryModel
CompType(mesh),
ODE(),
Ys_(this->solidThermo().composition().Y()),
pyrolisisGases_
(
mesh.lookupObject<dictionary>
("thermophysicalProperties").lookup("gaseousSpecies")
),
reactions_
(
dynamic_cast<const reactingSolidMixture<SolidThermo>& >
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
)
),
pyrolisisGases_(reactions_[0].gasSpecies()),
solidThermo_
(
dynamic_cast<const reactingSolidMixture<SolidThermo>& >
dynamic_cast<const reactingMixture<SolidThermo>& >
(
this->solidThermo()
).solidData()
).speciesData()
),
gasThermo_(pyrolisisGases_.size()),
nGases_(pyrolisisGases_.size()),
@ -184,7 +181,7 @@ ODESolidChemistryModel
mesh.lookupObject<dictionary>
(
"thermophysicalProperties"
).subDict(pyrolisisGases_[gasI] + "Coeffs");
).subDict(pyrolisisGases_[gasI]);
gasThermo_.set
(
@ -193,14 +190,13 @@ ODESolidChemistryModel
);
}
Info<< "ODESolidChemistryModel: Number of solids = " << nSolids_
<< " and reactions = " << nReaction_ << endl;
Info<< "solidChemistryModel: Number of solids = " << nSolids_ << endl;
Info<< "Number of gases from pyrolysis = " << nGases_ << endl;
Info<< "solidChemistryModel: Number of gases = " << nGases_ << endl;
forAll(reactions_, i)
{
Info<< indent << "Reaction " << i << nl << reactions_[i] << nl;
Info<< indent << reactions_[i] << nl;
}
}
@ -235,23 +231,23 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
forAll(reactions_, i)
{
const solidReaction& R = reactions_[i];
const Reaction<SolidThermo>& R = reactions_[i];
scalar omegai = omega
(
R, c, T, p, pf, cf, lRef, pr, cr, rRef
);
scalar rhoL = 0.0;
forAll(R.slhs(), s)
forAll(R.lhs(), s)
{
label si = R.slhs()[s];
label si = R.lhs()[s].index;
om[si] -= omegai;
rhoL = solidThermo_[si].rho(p, T);
}
scalar sr = 0.0;
forAll(R.srhs(), s)
forAll(R.rhs(), s)
{
label si = R.srhs()[s];
label si = R.rhs()[s].index;
scalar rhoR = solidThermo_[si].rho(p, T);
sr = rhoR/rhoL;
om[si] += sr*omegai;
@ -263,7 +259,7 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
}
forAll(R.grhs(), g)
{
label gi = R.grhs()[g];
label gi = R.grhs()[g].index;
om[gi + nSolids_] += (1.0 - sr)*omegai;
}
}
@ -276,7 +272,7 @@ template<class CompType, class SolidThermo, class GasThermo>
Foam::scalar
Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
(
const solidReaction& R,
const Reaction<SolidThermo>& R,
const scalarField& c,
const scalar T,
const scalar p,
@ -299,16 +295,17 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::omega
scalar kf = R.kf(p, T, c1);
const scalar exponent = R.nReact();
//const scalar exponent = R.nReact();
const label Nl = R.slhs().size();
const label Nl = R.lhs().size();
for (label s=0; s<Nl; s++)
{
label si = R.slhs()[s];
label si = R.lhs()[s].index;
const scalar exp = R.lhs()[si].exponent;
kf *=
pow(c1[si]/Ys0_[si][cellI], exponent)
pow(c1[si]/Ys0_[si][cellI], exp)
*(Ys0_[si][cellI]);
}
@ -390,18 +387,18 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
for (label ri=0; ri<reactions_.size(); ri++)
{
const solidReaction& R = reactions_[ri];
const Reaction<SolidThermo>& R = reactions_[ri];
scalar kf0 = R.kf(p, T, c2);
forAll(R.slhs(), j)
forAll(R.lhs(), j)
{
label sj = R.slhs()[j];
label sj = R.lhs()[j].index;
scalar kf = kf0;
forAll(R.slhs(), i)
forAll(R.lhs(), i)
{
label si = R.slhs()[i];
scalar exp = R.nReact();
label si = R.lhs()[i].index;
scalar exp = R.lhs()[i].exponent;
if (i == j)
{
if (exp < 1.0)
@ -428,14 +425,14 @@ void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::jacobian
}
}
forAll(R.slhs(), i)
forAll(R.lhs(), i)
{
label si = R.slhs()[i];
label si = R.lhs()[i].index;
dfdc[si][sj] -= kf;
}
forAll(R.srhs(), i)
forAll(R.rhs(), i)
{
label si = R.srhs()[i];
label si = R.rhs()[i].index;
dfdc[si][sj] += kf;
}
}
@ -551,7 +548,8 @@ ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::nEqns() const
template<class CompType, class SolidThermo, class GasThermo>
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::calculate()
void Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::
calculate()
{
if (!this->chemistry_)
{

12
src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModel.H
View File

@ -38,7 +38,7 @@ SourceFiles
#ifndef ODESolidChemistryModel_H
#define ODESolidChemistryModel_H
#include "solidReaction.H"
#include "Reaction.H"
#include "ODE.H"
#include "volFieldsFwd.H"
#include "DimensionedField.H"
@ -72,12 +72,12 @@ protected:
//- Reference to solid mass fractions
PtrList<volScalarField>& Ys_;
//- Reactions
const PtrList<Reaction<SolidThermo> >& reactions_;
//- List of gas species present in reaction system
speciesTable pyrolisisGases_;
//- Reactions
const PtrList<solidReaction>& reactions_;
//- Thermodynamic data of solids
const PtrList<SolidThermo>& solidThermo_;
@ -149,7 +149,7 @@ public:
// Member Functions
//- The reactions
inline const PtrList<solidReaction>& reactions() const;
inline const PtrList<Reaction<SolidThermo> >& reactions() const;
//- Thermodynamic data of gases
inline const PtrList<GasThermo>& gasThermo() const;
@ -180,7 +180,7 @@ public:
// species and charateristic times
virtual scalar omega
(
const solidReaction& r,
const Reaction<SolidThermo>& r,
const scalarField& c,
const scalar T,
const scalar p,

2
src/thermophysicalModels/solidChemistryModel/ODESolidChemistryModel/ODESolidChemistryModelI.H
View File

@ -45,7 +45,7 @@ Foam::ODESolidChemistryModel<CompType, SolidThermo, GasThermo>::RRg()
template<class CompType, class SolidThermo, class GasThermo>
inline const Foam::PtrList<Foam::solidReaction>&
inline const Foam::PtrList<Foam::Reaction<SolidThermo> >&
Foam::ODESolidChemistryModel<CompType, SolidThermo,GasThermo>::reactions() const
{
return reactions_;

1
src/thermophysicalModels/solidChemistryModel/solidChemistrySolver/makeSolidChemistrySolvers.C
View File

@ -51,6 +51,7 @@ namespace Foam
hExponentialSolidThermoPhysics,
gasThermoPhysics
)
}

1
src/thermophysicalModels/solidSpecie/Make/files
View File

@ -1,4 +1,3 @@
reaction/Reactions/solidReaction/solidReaction.C
reaction/reactions/makeSolidReactions.C
LIB = $(FOAM_LIBBIN)/libsolidSpecie

3
src/thermophysicalModels/solidSpecie/Make/options
View File

@ -1,2 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
LIB_LIBS = \
-lspecie

92
src/thermophysicalModels/solidSpecie/reaction/Reactions/IrreversibleSolidReaction/IrreversibleSolidReaction.C
View File

@ -1,92 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "IrreversibleSolidReaction.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class ReactionRate>
Foam::IrreversibleSolidReaction<ReactionRate>::IrreversibleSolidReaction
(
const solidReaction& reaction,
const ReactionRate& k,
const scalar nReact
)
:
solidReaction(reaction),
k_(k),
nReact_(nReact)
{}
template<class ReactionRate>
Foam::IrreversibleSolidReaction<ReactionRate>::IrreversibleSolidReaction
(
const speciesTable& components,
Istream& is,
const speciesTable& pyrolysisGases
)
:
solidReaction(components, is, pyrolysisGases),
k_(components, is),
nReact_(readScalar(is))
{
is.readEnd("solidArrheniusReactionRate(Istream&)");
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ReactionRate>
Foam::scalar Foam::IrreversibleSolidReaction<ReactionRate>::kf
(
const scalar p,
const scalar T,
const scalarField& c
) const
{
return k_(p, T, c);
}
template<class ReactionRate>
Foam::scalar Foam::IrreversibleSolidReaction<ReactionRate>::nReact() const
{
return nReact_;
}
template<class ReactionRate>
void Foam::IrreversibleSolidReaction<ReactionRate>::write
(
Ostream& os
) const
{
solidReaction::write(os);
os << token::SPACE << "Reaction order: " << nReact_ << nl << k_;
}
// ************************************************************************* //

323
src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.C
View File

@ -26,299 +26,118 @@ License
#include "solidReaction.H"
#include "DynamicList.H"
/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
namespace Foam
{
defineTypeNameAndDebug(solidReaction, 0);
defineRunTimeSelectionTable(solidReaction, Istream);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solidReaction::solidReaction
template<class ReactionThermo>
Foam::solidReaction<ReactionThermo>::solidReaction
(
const speciesTable& componets,
const Reaction<ReactionThermo>& reaction,
const speciesTable& pyrolisisGases,
const List<label>& slhs,
const List<label>& srhs,
const List<label>& grhs
const List<specieCoeffs>& glhs,
const List<specieCoeffs>& grhs
)
:
components_(componets),
Reaction<ReactionThermo>(reaction),
pyrolisisGases_(pyrolisisGases),
slhs_(slhs),
srhs_(srhs),
glhs_(glhs),
grhs_(grhs)
{}
Foam::solidReaction::solidReaction
template<class ReactionThermo>
Foam::solidReaction<ReactionThermo>::solidReaction
(
const solidReaction& r,
const speciesTable& componets,
const solidReaction<ReactionThermo>& r,
const speciesTable& pyrolisisGases
)
:
components_(componets),
Reaction<ReactionThermo>(r),
pyrolisisGases_(pyrolisisGases),
slhs_(r.slhs_),
srhs_(r.srhs_),
glhs_(r.glhs_),
grhs_(r.grhs_)
{}
Foam::solidReaction::solidReaction
(
const speciesTable& components,
Istream& is,
const speciesTable& pyrolisisGases
)
:
components_(components),
pyrolisisGases_(pyrolisisGases)
{
setLRhs(is);
}
Foam::label Foam::solidReaction::componentIndex
(
bool& isGas,
Istream& is
)
{
token t(is);
if (t.isWord())
{
word componentName = t.wordToken();
size_t i = componentName.find('=');
if (i != word::npos)
{
string exponentStr = componentName
(
i + 1,
componentName.size() - i - 1
);
componentName = componentName(0, i);
}
if (components_.contains(componentName))
{
isGas = false;
return (components_[componentName]);
}
else if (pyrolisisGases_.contains(componentName))
{
isGas = true;
return (pyrolisisGases_[componentName]);
}
else
{
FatalIOErrorIn
(
"solidReaction::componentIndex(bool&, Istream& is)",
is
)
<< "Cannot find component" << componentName
<< "in tables :" << pyrolisisGases_ << " or "
<< components_
<< exit(FatalIOError);
return -1;
}
}
else
{
FatalIOErrorIn("solidReaction::componentIndex(bool&, Istream& is)", is)
<< "Expected a word but found " << t.info()
<< exit(FatalIOError);
return -1;
}
}
void Foam::solidReaction::setLRhs(Istream& is)
{
DynamicList<label> dlsrhs;
label index = 0;
while (is)
{
bool isGas = false;
index = componentIndex(isGas, is);
if (index != -1)
{
dlsrhs.append(index);
token t(is);
if (t.isPunctuation())
{
if (t == token::ADD)
{
if (isGas)
{
grhs_ = dlsrhs.shrink();
dlsrhs.clear();
//is.putBack(t);
//return;
}
else
{
srhs_ = dlsrhs.shrink();
dlsrhs.clear(); //is.putBack(t);
//return;
}
}
else if (t == token::ASSIGN)
{
if (isGas)
{
Info << "Pyrolysis Gases should appear on lhs of the"
"reaction" << endl;
}
else
{
slhs_ = dlsrhs.shrink();
dlsrhs.clear();
}
}
else if (isGas)
{
grhs_ = dlsrhs.shrink();
is.putBack(t);
return;
}
else
{
srhs_ = dlsrhs.shrink();
is.putBack(t);
return;
}
}
else if (isGas)
{
grhs_ = dlsrhs.shrink();
is.putBack(t);
return;
}
else
{
srhs_ = dlsrhs.shrink();
is.putBack(t);
return;
}
}
else
{
FatalIOErrorIn("solidReaction::lsrhs(Istream& is)", is)
<< "Cannot find component in tables"
<< exit(FatalIOError);
}
}
FatalIOErrorIn("solidReaction::lsrhs(Istream& is)", is)
<< "Cannot continue reading reaction data from stream"
<< exit(FatalIOError);
}
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::solidReaction> Foam::solidReaction::New
template<class ReactionThermo>
Foam::solidReaction<ReactionThermo>::solidReaction
(
const speciesTable& species,
Istream& is,
const speciesTable& pyrolisisGases
const HashPtrTable<ReactionThermo>& thermoDatabase,
Istream& is
)
:
Reaction<ReactionThermo>(species, thermoDatabase, is),
pyrolisisGases_(),
glhs_(),
grhs_()
{
if (is.eof())
{
FatalIOErrorIn
(
"solidReaction::New(const speciesTable& species,"
" const HashPtrTable& thermoDatabase, Istream&)",
is
) << "solidReaction type not specified" << nl << nl
<< "Valid solidReaction types are :" << endl
<< IstreamConstructorTablePtr_->sortedToc()
<< exit(FatalIOError);
}
const word reactionTypeName(is);
IstreamConstructorTable::iterator cstrIter
= IstreamConstructorTablePtr_->find(reactionTypeName);
if (cstrIter == IstreamConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"solidReaction::New(const speciesTable& species,"
" const HashPtrTable& thermoDatabase, Istream&)",
is
) << "Unknown reaction type "
<< reactionTypeName << nl << nl
<< "Valid reaction types are :" << endl
<< IstreamConstructorTablePtr_->sortedToc()
<< exit(FatalIOError);
}
return autoPtr<solidReaction>
notImplemented
(
cstrIter()(species, is, pyrolisisGases)
"template<class ReactionThermo>"
"Foam::solidReaction<ReactionThermo>::solidReaction"
"("
" const speciesTable& species,"
" const HashPtrTable<ReactionThermo>& thermoDatabase,"
" Istream& is"
")"
);
}
template<class ReactionThermo>
Foam::solidReaction<ReactionThermo>::solidReaction
(
const speciesTable& species,
const HashPtrTable<ReactionThermo>& thermoDatabase,
const dictionary& dict
)
:
Reaction<ReactionThermo>(species, thermoDatabase, dict),
pyrolisisGases_(dict.parent().parent().lookup("gaseousSpecies")),
glhs_(),
grhs_()
{
this->setLRhs
(
IStringStream(dict.lookup("reaction"))(),
pyrolisisGases_,
glhs_,
grhs_
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::solidReaction::write(Ostream& os) const
template<class ReactionThermo>
const Foam::List<typename Foam::solidReaction<ReactionThermo>::specieCoeffs>&
Foam::solidReaction<ReactionThermo>::glhs() const
{
os << type() << nl << " ";
forAll(slhs_, i)
{
os << components_[slhs_[i]];
}
os << " = ";
forAll(srhs_, i)
{
os << components_[srhs_[i]];
}
os << " + ";
forAll(grhs_, i)
{
os << pyrolisisGases_[grhs_[i]];
}
os << endl << " ";
return glhs_;
}
Foam::scalar Foam::solidReaction::kf
(
const scalar p,
const scalar T,
const scalarField& c
) const
template<class ReactionThermo>
const Foam::List<typename Foam::Reaction<ReactionThermo>::specieCoeffs>&
Foam::solidReaction<ReactionThermo>::grhs() const
{
return 0.0;
return grhs_;
}
Foam::scalar Foam::solidReaction::nReact() const
template<class ReactionThermo>
const Foam::speciesTable& Foam::solidReaction<ReactionThermo>::
gasSpecies() const
{
return 1.0;
return pyrolisisGases_;
}
template<class ReactionThermo>
void Foam::solidReaction<ReactionThermo>::write(Ostream& os) const
{
Reaction<ReactionThermo>::write(os);
}

166
src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H
View File

@ -38,9 +38,7 @@ SourceFiles
#define solidReaction_H
#include "speciesTable.H"
#include "scalarField.H"
#include "typeInfo.H"
#include "runTimeSelectionTables.H"
#include "Reaction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -48,44 +46,40 @@ namespace Foam
{
// Forward declaration of friend functions and operators
template<class ReactionThermo>
class solidReaction;
inline Ostream& operator<<(Ostream&, const solidReaction&);
template<class ReactionThermo>
inline Ostream& operator<<(Ostream&, const solidReaction<ReactionThermo>&);
/*---------------------------------------------------------------------------*\
Class solidReaction Declaration
\*---------------------------------------------------------------------------*/
template<class ReactionThermo>
class solidReaction
:
public Reaction<ReactionThermo>
{
private:
// Private data
//- List of solid names present in reaction system
const speciesTable& components_;
typedef typename Reaction<ReactionThermo>::specieCoeffs specieCoeffs;
//- List of gas species present in reaction system
speciesTable pyrolisisGases_;
//- Solid components index for the left-hand-side of the reaction
List<label> slhs_;
//- Gas specie index for the left-hand-side of the reaction
List<specieCoeffs> glhs_;
//- Solid components index for the right-hand-side of the reaction
List<label> srhs_;
//- Specie index for the right-hand-side of the reaction
List<label> grhs_;
//- Gas specie index for the right-hand-side of the reaction
List<specieCoeffs> grhs_;
// Private Member Functions
//- Set rhs and lhs of the reaction
void setLRhs(Istream&);
//- Look for the component index in the reaction
label componentIndex(bool& isGas, Istream& is);
//- Disallow default bitwise assignment
void operator=(const solidReaction&);
@ -94,58 +88,7 @@ private:
public:
//- Runtime type information
TypeName("Reaction");
// Declare run-time constructor selection tables
declareRunTimeSelectionTable
(
autoPtr,
solidReaction,
Istream,
(
const speciesTable& components,
Istream& is,
const speciesTable& pyrolysisGases
),
(components, is, pyrolysisGases)
);
// Public classes
//- Class used for the read-construction of PtrLists of reaction
class iNew
{
const speciesTable& components_;
speciesTable pyrolisisGases_;
public:
iNew
(
const speciesTable& components,
Istream& pyrolisisGases
)
:
components_(components),
pyrolisisGases_(pyrolisisGases)
{}
autoPtr<solidReaction> operator()(Istream& is) const
{
return autoPtr<solidReaction>
(
solidReaction::New
(
components_,
is,
pyrolisisGases_
)
);
}
};
TypeName("SolidReaction");
// Constructors
@ -153,48 +96,59 @@ public:
//- Construct from components
solidReaction
(
const speciesTable& components,
const Reaction<ReactionThermo>& reaction,
const speciesTable& pyrolisisGases,
const List<label>& slhs,
const List<label>& srhs,
const List<label>& grhs
const List<specieCoeffs>& glhs,
const List<specieCoeffs>& grhs
);
//- Construct as copy given new speciesTable
solidReaction
(
const solidReaction&,
const speciesTable& components,
const solidReaction<ReactionThermo>&,
const speciesTable& pyrolisisGases
);
//- Construct from Istream
solidReaction
(
const speciesTable& components,
Istream& is,
const speciesTable& pyrolisisGases
const speciesTable& pyrolisisGases,
const HashPtrTable<ReactionThermo>& thermoDatabase,
Istream& is
);
//- Construct from dictionary
solidReaction
(
const speciesTable& species,
const HashPtrTable<ReactionThermo>& thermoDatabase,
const dictionary& dict
);
//- Construct and return a clone
virtual autoPtr<solidReaction > clone() const
virtual autoPtr<Reaction<ReactionThermo> > clone() const
{
return autoPtr<solidReaction >
return autoPtr<Reaction<ReactionThermo> >
(
new solidReaction(*this)
new solidReaction<ReactionThermo>(*this)
);
}
// Selectors
//- Return a pointer to a new patchField created on freestore from input
static autoPtr<solidReaction > New
//- Construct and return a clone with new speciesTable
virtual autoPtr<Reaction<ReactionThermo> > clone
(
const speciesTable& components,
Istream&,
const speciesTable& pyrolisisGases
);
const speciesTable& species
) const
{
return autoPtr<Reaction<ReactionThermo> >
(
new solidReaction<ReactionThermo>(*this, species)
);
}
//- Destructor
@ -206,23 +160,13 @@ public:
// Access
inline const List<label>& slhs() const;
inline const List<label>& srhs() const;
inline const List<label>& grhs() const;
inline const speciesTable& pyrolisisGases() const;
// - Acces to gas components of the reaction
virtual const List<specieCoeffs>& grhs() const;
virtual const List<specieCoeffs>& glhs() const;
// solidReaction rate coefficients
virtual scalar kf
(
const scalar p,
const scalar T,
const scalarField& c
) const;
virtual scalar nReact() const;
// - Access to gas specie list
virtual const speciesTable& gasSpecies() const;
//- Write
@ -231,10 +175,10 @@ public:
// Ostream Operator
friend Ostream& operator<<
friend Ostream& operator<< <ReactionThermo>
(
Ostream&,
const solidReaction&
const solidReaction<ReactionThermo>&
);
};
@ -249,6 +193,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "solidReaction.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

29
src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReactionI.H
View File

@ -30,38 +30,13 @@ License
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const List<label>& solidReaction::slhs() const
{
return slhs_;
}
inline const List<label>& solidReaction::srhs() const
{
return srhs_;
}
inline const List<label>& solidReaction::grhs() const
{
return grhs_;
}
inline const speciesTable& solidReaction::pyrolisisGases() const
{
return pyrolisisGases_;
}
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class ReactionThermo>
inline Ostream& operator<<
(
Ostream& os,
const solidReaction& r
const solidReaction<ReactionThermo>& r
)
{
r.write(os);

18
src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H
View File

@ -75,13 +75,25 @@ public:
Istream& is
);
//- Construct from dictionary
inline solidArrheniusReactionRate
(
const speciesTable& species,
const dictionary& dict
);
//- Destructor
virtual ~solidArrheniusReactionRate()
{}
// Member Functions
//- Return the type name
static word type()
{
return "SolidArrhenius";
return "Arrhenius";
}
inline scalar operator()
@ -92,6 +104,10 @@ public:
) const;
//- Write to stream
inline void write(Ostream& os) const;
// Ostream Operator
inline friend Ostream& operator<<

21
src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H
View File

@ -30,6 +30,7 @@ inline Foam::solidArrheniusReactionRate::solidArrheniusReactionRate
const scalar A,
const scalar Ta,
const scalar Tcrit
//const scalar nReact
)
:
A_(A),
@ -50,6 +51,18 @@ inline Foam::solidArrheniusReactionRate::solidArrheniusReactionRate
{}
inline Foam::solidArrheniusReactionRate::solidArrheniusReactionRate
(
const speciesTable&,
const dictionary& dict
)
:
A_(readScalar(dict.lookup("A"))),
Ta_(readScalar(dict.lookup("Ta"))),
Tcrit_(readScalar(dict.lookup("Tcrit")))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::scalar Foam::solidArrheniusReactionRate::operator()
@ -74,6 +87,14 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::operator()
}
inline void Foam::solidArrheniusReactionRate::write(Ostream& os) const
{
os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;
os.writeKeyword("Ta") << Ta_ << token::END_STATEMENT << nl;
os.writeKeyword("Tcrit") << Tcrit_ << token::END_STATEMENT << nl;
}
inline Foam::Ostream& Foam::operator<<
(
Ostream& os,

40
src/thermophysicalModels/solidSpecie/reaction/reactions/makeSolidReaction.H
View File

@ -33,7 +33,9 @@ Description
#define makeSolidReactionThermo_H
#include "solidReaction.H"
#include "IrreversibleSolidReaction.H"
#include "IrreversibleReaction.H"
#include "Reaction.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -43,32 +45,46 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#define makeReaction(ReactionType, ReactionRate) \
#define makeSolidReaction(ReactionType, Thermo, ReactionRate) \
\
typedef solidReaction Reaction; \
typedef solidReaction<Thermo> solidReaction##Thermo; \
\
typedef ReactionType<ReactionRate> \
ReactionType##ReactionRate; \
typedef Reaction<Thermo> Reaction##Thermo; \
\
typedef ReactionType<solidReaction, Thermo, ReactionRate> \
ReactionType##Thermo##ReactionRate; \
\
defineTemplateRunTimeSelectionTable(Reaction##Thermo, Istream); \
defineTemplateRunTimeSelectionTable(Reaction##Thermo, dictionary); \
\
defineTemplateTypeNameAndDebug(solidReaction##Thermo, 0); \
defineTemplateTypeNameAndDebug(Reaction##Thermo, 0); \
\
template<> \
const word ReactionType##ReactionRate::typeName \
const word ReactionType##Thermo##ReactionRate::typeName \
( \
ReactionType::typeName_() \
+ ReactionRate::type() \
+ Reaction::typeName_() \
+ solidReaction##Thermo::typeName_() \
); \
\
addToRunTimeSelectionTable \
( \
Reaction, \
ReactionType##ReactionRate, \
Reaction##Thermo, \
ReactionType##Thermo##ReactionRate, \
Istream \
); \
\
addToRunTimeSelectionTable \
( \
Reaction##Thermo, \
ReactionType##Thermo##ReactionRate, \
dictionary \
);
#define makeIRReactions(ReactionRate) \
#define makeSolidIRReactions(Thermo, ReactionRate) \
\
makeReaction(IrreversibleSolidReaction, ReactionRate)
makeSolidReaction(IrreversibleReaction, Thermo, ReactionRate)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

9
src/thermophysicalModels/solidSpecie/reaction/reactions/makeSolidReactions.C
View File

@ -25,6 +25,7 @@ License
#include "makeSolidReaction.H"
#include "solidArrheniusReactionRate.H"
#include "solidThermoPhysicsTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -33,7 +34,13 @@ namespace Foam
// * * * * * * * * * * * * * Make Solid reactions * * * * * * * * * * * * //
makeIRReactions(solidArrheniusReactionRate)
makeSolidIRReactions(hConstSolidThermoPhysics, solidArrheniusReactionRate)
makeSolidIRReactions
(
hExponentialSolidThermoPhysics,
solidArrheniusReactionRate
)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

50
src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransport.H
View File

@ -26,7 +26,7 @@ Class
Description
Constant properties Transport package.
Templated into a given thermodynamics package (needed for thermal
Templated into a given Thermodynamics package (needed for thermal
conductivity).
SourceFiles
@ -44,20 +44,20 @@ SourceFiles
namespace Foam
{
template<class thermo> class constAnIsoSolidTransport;
template<class Thermo> class constAnIsoSolidTransport;
template<class thermo>
inline constAnIsoSolidTransport<thermo> operator*
template<class Thermo>
inline constAnIsoSolidTransport<Thermo> operator*
(
const scalar,
const constAnIsoSolidTransport<thermo>&
const constAnIsoSolidTransport<Thermo>&
);
template<class thermo>
template<class Thermo>
Ostream& operator<<
(
Ostream&,
const constAnIsoSolidTransport<thermo>&
const constAnIsoSolidTransport<Thermo>&
);
@ -65,10 +65,10 @@ Ostream& operator<<
Class constAnIsoSolidTransport Declaration
\*---------------------------------------------------------------------------*/
template<class thermo>
template<class Thermo>
class constAnIsoSolidTransport
:
public thermo
public Thermo
{
// Private data
@ -79,7 +79,7 @@ class constAnIsoSolidTransport
// Private Member Functions
//- Construct from components
inline constAnIsoSolidTransport(const thermo& t, const vector kappa);
inline constAnIsoSolidTransport(const Thermo& t, const vector kappa);
public:
@ -93,18 +93,36 @@ public:
const constAnIsoSolidTransport&
);
//- Construct from Istream
//constAnIsoSolidTransport(Istream&);
//- Construct from dictionary
constAnIsoSolidTransport(const dictionary&);
// Selector from dictionary
inline static autoPtr<constAnIsoSolidTransport> New
(
const dictionary& dict
);
// Member functions
//- Return the instantiated type name
static word typeName()
{
return "constAnIso<" + Thermo::typeName() + '>';
}
//- Isotropic thermal conductivity [W/mK]
inline scalar kappa(const scalar T) const;
inline scalar kappa(const scalar p, const scalar T) const;
//- Un-isotropic thermal conductivity [W/mK]
inline vector Kappa(const scalar T) const;
inline vector Kappa(const scalar p, const scalar T) const;
//- Dynamic viscosity [kg/ms]
inline scalar mu(const scalar p, const scalar T) const;
//- Thermal diffusivity of enthalpy [kg/ms]
inline vector alphah(const scalar p, const scalar T) const;
//- Write to Ostream
void write(Ostream& os) const;
@ -121,7 +139,7 @@ public:
// Friend operators
friend constAnIsoSolidTransport operator* <thermo>
friend constAnIsoSolidTransport operator* <Thermo>
(
const scalar,
const constAnIsoSolidTransport&
@ -130,7 +148,7 @@ public:
// Ostream Operator
friend Ostream& operator<< <thermo>
friend Ostream& operator<< <Thermo>
(
Ostream&,
const constAnIsoSolidTransport&

92
src/thermophysicalModels/solidSpecie/transport/const/constAnIsoSolidTransportI.H
View File

@ -25,53 +25,89 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class thermo>
inline Foam::constAnIsoSolidTransport<thermo>::constAnIsoSolidTransport
template<class Thermo>
inline Foam::constAnIsoSolidTransport<Thermo>::constAnIsoSolidTransport
(
const thermo& t,
const Thermo& t,
const vector kappa
)
:
thermo(t),
Thermo(t),
kappa_(kappa)
{}
template<class thermo>
inline Foam::constAnIsoSolidTransport<thermo>::constAnIsoSolidTransport
template<class Thermo>
inline Foam::constAnIsoSolidTransport<Thermo>::constAnIsoSolidTransport
(
const word& name,
const constAnIsoSolidTransport& ct
)
:
thermo(name, ct),
Thermo(name, ct),
kappa_(ct.kappa_)
{}
template<class Thermo>
inline Foam::autoPtr<Foam::constAnIsoSolidTransport<Thermo> >
Foam::constAnIsoSolidTransport<Thermo>::New
(
const dictionary& dict
)
{
return autoPtr<constAnIsoSolidTransport<Thermo> >
(
new constAnIsoSolidTransport<Thermo>(dict)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class thermo>
inline Foam::scalar Foam::constAnIsoSolidTransport<thermo>::
kappa(const scalar T) const
template<class Thermo>
inline Foam::scalar Foam::constAnIsoSolidTransport<Thermo>::
kappa(const scalar p, const scalar T) const
{
return mag(kappa_);
}
template<class thermo>
inline Foam::vector Foam::constAnIsoSolidTransport<thermo>::
Kappa(const scalar T) const
template<class Thermo>
inline Foam::vector Foam::constAnIsoSolidTransport<Thermo>::
Kappa(const scalar p, const scalar T) const
{
return kappa_;
}
template<class Thermo>
inline Foam::scalar Foam::constAnIsoSolidTransport<Thermo>::
mu(const scalar p, const scalar T) const
{
notImplemented
(
"Foam::scalar Foam::constAnIsoSolidTransport<Thermo>mu::"
"("
" const scalar p, const scalar T"
") const"
);
return scalar(0);
}
template<class Thermo>
inline Foam::vector Foam::constAnIsoSolidTransport<Thermo>::
alphah(const scalar p, const scalar T) const
{
return kappa_/this->Cpv(p, T);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class thermo>
inline Foam::constAnIsoSolidTransport<thermo>&
Foam::constAnIsoSolidTransport<thermo>::operator=
template<class Thermo>
inline Foam::constAnIsoSolidTransport<Thermo>&
Foam::constAnIsoSolidTransport<Thermo>::operator=
(
const constAnIsoSolidTransport<thermo>& ct
const constAnIsoSolidTransport<Thermo>& ct
)
{
kappa_ = ct.kappa_;
@ -80,10 +116,10 @@ Foam::constAnIsoSolidTransport<thermo>::operator=
}
template<class thermo>
inline void Foam::constAnIsoSolidTransport<thermo>::operator+=
template<class Thermo>
inline void Foam::constAnIsoSolidTransport<Thermo>::operator+=
(
const constAnIsoSolidTransport<thermo>& ct
const constAnIsoSolidTransport<Thermo>& ct
)
{
scalar molr1 = this->nMoles();
@ -95,10 +131,10 @@ inline void Foam::constAnIsoSolidTransport<thermo>::operator+=
}
template<class thermo>
inline void Foam::constAnIsoSolidTransport<thermo>::operator-=
template<class Thermo>
inline void Foam::constAnIsoSolidTransport<Thermo>::operator-=
(
const constAnIsoSolidTransport<thermo>& ct
const constAnIsoSolidTransport<Thermo>& ct
)
{
scalar molr1 = this->nMoles();
@ -113,16 +149,16 @@ inline void Foam::constAnIsoSolidTransport<thermo>::operator-=
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class thermo>
inline Foam::constAnIsoSolidTransport<thermo> Foam::operator*
template<class Thermo>
inline Foam::constAnIsoSolidTransport<Thermo> Foam::operator*
(
const scalar s,
const constAnIsoSolidTransport<thermo>& ct
const constAnIsoSolidTransport<Thermo>& ct
)
{
return constAnIsoSolidTransport<thermo>
return constAnIsoSolidTransport<Thermo>
(
s*static_cast<const thermo&>(ct),
s*static_cast<const Thermo&>(ct),
ct.kappa_
);
}

18
src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.C
View File

@ -28,26 +28,26 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class thermo>
Foam::constIsoSolidTransport<thermo>::constIsoSolidTransport
template<class Thermo>
Foam::constIsoSolidTransport<Thermo>::constIsoSolidTransport
(
const dictionary& dict
)
:
thermo(dict),
Thermo(dict),
kappa_(readScalar(dict.subDict("transport").lookup("kappa")))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class thermo>
void Foam::constIsoSolidTransport<thermo>::constIsoSolidTransport::write
template<class Thermo>
void Foam::constIsoSolidTransport<Thermo>::constIsoSolidTransport::write
(
Ostream& os
) const
{
thermo::write(os);
Thermo::write(os);
dictionary dict("transport");
dict.add("kappa", kappa_);
@ -57,14 +57,14 @@ void Foam::constIsoSolidTransport<thermo>::constIsoSolidTransport::write
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
template<class thermo>
template<class Thermo>
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const constIsoSolidTransport<thermo>& ct
const constIsoSolidTransport<Thermo>& ct
)
{
operator<<(os, static_cast<const thermo&>(ct));
operator<<(os, static_cast<const Thermo&>(ct));
os << tab << ct.kappa_;
os.check

17
src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransport.H
View File

@ -97,6 +97,12 @@ public:
//- Construct from Istream
constIsoSolidTransport(const dictionary& dict);
// Selector from dictionary
inline static autoPtr<constIsoSolidTransport> New
(
const dictionary& dict
);
// Member functions
@ -107,10 +113,17 @@ public:
}
//- Isotropic thermal conductivity [W/mK]
inline scalar kappa(const scalar T) const;
inline scalar kappa(const scalar p, const scalar T) const;
//- Un-isotropic thermal conductivity [W/mK]
inline vector Kappa(const scalar T) const;
inline vector Kappa(const scalar p, const scalar T) const;
//- Dynamic viscosity [kg/ms]
inline scalar mu(const scalar p, const scalar T) const;
//- Thermal diffusivity of enthalpy [kg/ms]
inline scalar alphah(const scalar p, const scalar T) const;
//- Write to Ostream
void write(Ostream& os) const;

39
src/thermophysicalModels/solidSpecie/transport/const/constIsoSolidTransportI.H
View File

@ -49,23 +49,58 @@ inline Foam::constIsoSolidTransport<thermo>::constIsoSolidTransport
{}
template<class Thermo>
inline Foam::autoPtr<Foam::constIsoSolidTransport<Thermo> >
Foam::constIsoSolidTransport<Thermo>::New
(
const dictionary& dict
)
{
return autoPtr<constIsoSolidTransport<Thermo> >
(
new constIsoSolidTransport<Thermo>(dict)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class thermo>
inline Foam::scalar Foam::constIsoSolidTransport<thermo>::
kappa(const scalar T) const
kappa(const scalar p, const scalar T) const
{
return kappa_;
}
template<class thermo>
inline Foam::vector Foam::constIsoSolidTransport<thermo>::
Kappa(const scalar T) const
Kappa(const scalar p, const scalar T) const
{
return vector(kappa_, kappa_, kappa_);
}
template<class thermo>
inline Foam::scalar Foam::constIsoSolidTransport<thermo>::
mu(const scalar p, const scalar T) const
{
notImplemented
(
"Foam::scalar Foam::constIsoSolidTransport<thermo>mu::"
"("
" const scalar p, const scalar T"
") const"
);
return scalar(0);
}
template<class thermo>
inline Foam::scalar Foam::constIsoSolidTransport<thermo>::
alphah(const scalar p, const scalar T) const
{
return kappa_/this->Cpv(p, T);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class thermo>

17
src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.H
View File

@ -103,10 +103,15 @@ public:
const exponentialSolidTransport&
);
//- Construct from dictionary
exponentialSolidTransport(const dictionary&);
// Selector from dictionary
inline static autoPtr<exponentialSolidTransport> New
(
const dictionary& dict
);
// Member functions
@ -117,10 +122,16 @@ public:
}
//- Thermal conductivity [W/mK]
inline scalar kappa(const scalar T) const;
inline scalar kappa(const scalar p, const scalar T) const;
//- Thermal conductivity [W/mK]
inline vector Kappa(const scalar T) const;
inline vector Kappa(const scalar p, const scalar T) const;
//- Dynamic viscosity [kg/ms]
inline scalar mu(const scalar p, const scalar T) const;
//- Thermal diffusivity of enthalpy [kg/ms]
inline scalar alphah(const scalar p, const scalar T) const;
//- Write to Ostream
void write(Ostream& os) const;

Some files were not shown because too many files have changed in this diff Show More