ENH: porosity model updates for moving meshes

This commit is contained in:
andy 2014-04-11 14:32:31 +01:00
parent 6459a1d179
commit 36fdb47877
12 changed files with 219 additions and 108 deletions

View file

@ -81,9 +81,10 @@ void Foam::porosityModels::DarcyForchheimer::calcTranformModelData()
{
forAll (cellZoneIDs_, zoneI)
{
D_[zoneI].setSize(1, tensor::zero);
F_[zoneI].setSize(1, tensor::zero);
D_[zoneI].setSize(1);
F_[zoneI].setSize(1);
D_[zoneI][0] = tensor::zero;
D_[zoneI][0].xx() = dXYZ_.value().x();
D_[zoneI][0].yy() = dXYZ_.value().y();
D_[zoneI][0].zz() = dXYZ_.value().z();
@ -91,6 +92,7 @@ void Foam::porosityModels::DarcyForchheimer::calcTranformModelData()
D_[zoneI][0] = coordSys_.R().transformTensor(D_[zoneI][0]);
// leading 0.5 is from 1/2*rho
F_[zoneI][0] = tensor::zero;
F_[zoneI][0].xx() = 0.5*fXYZ_.value().x();
F_[zoneI][0].yy() = 0.5*fXYZ_.value().y();
F_[zoneI][0].zz() = 0.5*fXYZ_.value().z();
@ -104,25 +106,65 @@ void Foam::porosityModels::DarcyForchheimer::calcTranformModelData()
{
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
D_[zoneI].setSize(cells.size(), tensor::zero);
F_[zoneI].setSize(cells.size(), tensor::zero);
D_[zoneI].setSize(cells.size());
F_[zoneI].setSize(cells.size());
forAll(cells, i)
{
D_[zoneI][i] = tensor::zero;
D_[zoneI][i].xx() = dXYZ_.value().x();
D_[zoneI][i].yy() = dXYZ_.value().y();
D_[zoneI][i].zz() = dXYZ_.value().z();
// leading 0.5 is from 1/2*rho
F_[zoneI][i] = tensor::zero;
F_[zoneI][i].xx() = 0.5*fXYZ_.value().x();
F_[zoneI][i].yy() = 0.5*fXYZ_.value().y();
F_[zoneI][i].zz() = 0.5*fXYZ_.value().z();
}
D_[zoneI] = coordSys_.R().transformTensor(D_[zoneI], cells);
F_[zoneI] = coordSys_.R().transformTensor(F_[zoneI], cells);
const coordinateRotation& R = coordSys_.R(mesh_, cells);
D_[zoneI] = R.transformTensor(D_[zoneI], cells);
F_[zoneI] = R.transformTensor(F_[zoneI], cells);
}
}
if (debug && mesh_.time().outputTime())
{
volTensorField Dout
(
IOobject
(
typeName + ":D",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedTensor("0", dXYZ_.dimensions(), tensor::zero)
);
volTensorField Fout
(
IOobject
(
typeName + ":F",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedTensor("0", fXYZ_.dimensions(), tensor::zero)
);
UIndirectList<tensor>(Dout, mesh_.cellZones()[cellZoneIDs_[0]]) = D_[0];
UIndirectList<tensor>(Fout, mesh_.cellZones()[cellZoneIDs_[0]]) = F_[0];
Dout.write();
Fout.write();
}
}

View file

@ -137,14 +137,16 @@ void Foam::porosityModels::fixedCoeff::calcTranformModelData()
{
forAll (cellZoneIDs_, zoneI)
{
alpha_[zoneI].setSize(1, tensor::zero);
beta_[zoneI].setSize(1, tensor::zero);
alpha_[zoneI].setSize(1);
beta_[zoneI].setSize(1);
alpha_[zoneI][0] = tensor::zero;
alpha_[zoneI][0].xx() = alphaXYZ_.value().x();
alpha_[zoneI][0].yy() = alphaXYZ_.value().y();
alpha_[zoneI][0].zz() = alphaXYZ_.value().z();
alpha_[zoneI][0] = coordSys_.R().transformTensor(alpha_[zoneI][0]);
beta_[zoneI][0] = tensor::zero;
beta_[zoneI][0].xx() = betaXYZ_.value().x();
beta_[zoneI][0].yy() = betaXYZ_.value().y();
beta_[zoneI][0].zz() = betaXYZ_.value().z();
@ -157,24 +159,26 @@ void Foam::porosityModels::fixedCoeff::calcTranformModelData()
{
const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
alpha_[zoneI].setSize(cells.size(), tensor::zero);
beta_[zoneI].setSize(cells.size(), tensor::zero);
alpha_[zoneI].setSize(cells.size());
beta_[zoneI].setSize(cells.size());
forAll(cells, i)
{
alpha_[zoneI][i] = tensor::zero;
alpha_[zoneI][i].xx() = alphaXYZ_.value().x();
alpha_[zoneI][i].yy() = alphaXYZ_.value().y();
alpha_[zoneI][i].zz() = alphaXYZ_.value().z();
beta_[zoneI][i] = tensor::zero;
beta_[zoneI][i].xx() = betaXYZ_.value().x();
beta_[zoneI][i].yy() = betaXYZ_.value().y();
beta_[zoneI][i].zz() = betaXYZ_.value().z();
}
alpha_[zoneI] =
coordSys_.R().transformTensor(alpha_[zoneI], cells);
const coordinateRotation& R = coordSys_.R(mesh_, cells);
beta_[zoneI] = coordSys_.R().transformTensor(beta_[zoneI], cells);
alpha_[zoneI] = R.transformTensor(alpha_[zoneI], cells);
beta_[zoneI] = R.transformTensor(beta_[zoneI], cells);
}
}
}

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -157,6 +157,4 @@ Foam::tmp<Foam::vectorField> Foam::cartesianCS::globalToLocal
}
// ************************************************************************* //

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -130,6 +130,12 @@ public:
Rtr_ = sphericalTensor::I;
}
//- Update the rotation for a list of cells
virtual void updateCells(const polyMesh&, const labelList&)
{
// do nothing
}
//- Return local-to-global transformation tensor
virtual const tensor& R() const
{

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -126,6 +126,12 @@ public:
Rtr_ = sphericalTensor::I;
}
//- Update the rotation for a list of cells
virtual void updateCells(const polyMesh&, const labelList&)
{
// do nothing
}
//- Return local-to-global transformation tensor
virtual const tensor& R() const
{

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -50,10 +50,9 @@ void Foam::axesRotation::calcTransform
const axisOrder& order
)
{
vector a = axis1 / mag(axis1);
vector a = axis1/mag(axis1);
vector b = axis2;
// Absorb minor nonorthogonality into axis2
b = b - (b & a)*a;
if (mag(b) < SMALL)
@ -64,36 +63,49 @@ void Foam::axesRotation::calcTransform
<< abort(FatalError);
}
b = b / mag(b);
vector c = a ^ b;
b = b/mag(b);
vector c = a^b;
tensor Rtr;
switch (order)
{
case e1e2:
{
Rtr = tensor(a, b, c);
break;
}
case e2e3:
{
Rtr = tensor(c, a, b);
break;
}
case e3e1:
{
Rtr = tensor(b, c, a);
break;
}
default:
FatalErrorIn("axesRotation::calcTransform()")
<< "programmer error" << endl
{
FatalErrorIn
(
"axesRotation::calcTransform"
"("
"const vector&,"
"const vector&,"
"const axisOrder&"
")"
)
<< "Unhandled axes specifictation" << endl
<< abort(FatalError);
Rtr = tensor::zero;
break;
}
}
// the global -> local transformation
// the global->local transformation
Rtr_ = Rtr;
// the local -> global transformation
// the local->global transformation
R_ = Rtr.T();
}
@ -154,15 +166,13 @@ Foam::axesRotation::axesRotation(const tensor& R)
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::vector Foam::axesRotation::transform(const vector& st) const
const Foam::tensorField& Foam::axesRotation::Tr() const
{
return (R_ & st);
}
Foam::vector Foam::axesRotation::invTransform(const vector& st) const
{
return (Rtr_ & st);
notImplemented
(
"const Foam::tensorField& axesRotation::Tr() const"
);
return *reinterpret_cast<const tensorField*>(0);
}
@ -175,6 +185,12 @@ Foam::tmp<Foam::vectorField> Foam::axesRotation::transform
}
Foam::vector Foam::axesRotation::transform(const vector& st) const
{
return (R_ & st);
}
Foam::tmp<Foam::vectorField> Foam::axesRotation::invTransform
(
const vectorField& st
@ -184,13 +200,9 @@ Foam::tmp<Foam::vectorField> Foam::axesRotation::invTransform
}
const Foam::tensorField& Foam::axesRotation::Tr() const
Foam::vector Foam::axesRotation::invTransform(const vector& st) const
{
notImplemented
(
"const Foam::tensorField& axesRotation::Tr() const"
);
return *reinterpret_cast<const tensorField*>(0);
return (Rtr_ & st);
}
@ -225,13 +237,15 @@ Foam::tmp<Foam::tensorField> Foam::axesRotation::transformTensor
notImplemented
(
"tmp<Foam::tensorField> axesRotation::transformTensor "
" const tensorField& st,"
" const labelList& cellMap "
"("
"const tensorField&,"
"const labelList&"
") const"
);
return tmp<tensorField>(NULL);
}
Foam::tmp<Foam::symmTensorField> Foam::axesRotation::transformVector
(
const vectorField& st

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -99,11 +99,7 @@ public:
axesRotation();
//- Construct from 2 axes
axesRotation
(
const vector& axis,
const vector& dir
);
axesRotation(const vector& axis, const vector& dir);
//- Construct from dictionary
axesRotation(const dictionary&);
@ -135,6 +131,12 @@ public:
Rtr_ = sphericalTensor::I;
}
//- Update the rotation for a list of cells
virtual void updateCells(const polyMesh&, const labelList&)
{
// do nothing
}
//- Return local-to-global transformation tensor
virtual const tensor& R() const
{

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -55,6 +55,7 @@ Description
#include "dictionary.H"
#include "runTimeSelectionTables.H"
#include "objectRegistry.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -135,6 +136,13 @@ public:
//- Reset rotation to an identity rotation
virtual void clear() = 0;
//- Update the rotation for a list of cells
virtual void updateCells
(
const polyMesh& mesh,
const labelList& cells
) = 0;
//- Return local-to-global transformation tensor
virtual const tensor& R() const = 0;

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -48,12 +48,33 @@ namespace Foam
);
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::localAxesRotation::init(const objectRegistry& obr)
{
const polyMesh& mesh = refCast<const polyMesh>(obr);
const vectorField& cc = mesh.cellCentres();
tensorField& R = Rptr_();
forAll(cc, cellI)
{
vector dir = cc[cellI] - origin_;
dir /= mag(dir) + VSMALL;
const axesRotation ar(e3_, dir);
R[cellI] = ar.R();
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::localAxesRotation::localAxesRotation
(
const dictionary& dict,
const objectRegistry& orb
const objectRegistry& obr
)
:
Rptr_(),
@ -69,27 +90,21 @@ Foam::localAxesRotation::localAxesRotation
// rotation axis
dict.lookup("e3") >> e3_;
const polyMesh& mesh = refCast<const polyMesh>(orb);
const polyMesh& mesh = refCast<const polyMesh>(obr);
Rptr_.reset
(
new tensorField(mesh.nCells())
);
init(dict, orb);
Rptr_.reset(new tensorField(mesh.nCells()));
init(obr);
}
Foam::localAxesRotation::localAxesRotation
(
const dictionary& dict
)
Foam::localAxesRotation::localAxesRotation(const dictionary& dict)
:
Rptr_(),
origin_(),
e3_()
{
FatalErrorIn("localAxesRotation(const dictionary&)")
<< " localAxesRotation can not be constructed from dictionary "
<< " localAxesRotation can not be constructed from dictionary "
<< " use the construtctor : "
"("
" const dictionary&, const objectRegistry&"
@ -109,6 +124,23 @@ void Foam::localAxesRotation::clear()
}
void Foam::localAxesRotation::updateCells
(
const polyMesh& mesh,
const labelList& cells
)
{
forAll(cells, i)
{
label cellI = cells[i];
vector dir = mesh.cellCentres()[cellI] - origin_;
dir /= mag(dir) + VSMALL;
Rptr_()[cellI] = axesRotation(e3_, dir).R();
}
}
Foam::vector Foam::localAxesRotation::transform(const vector& st) const
{
notImplemented
@ -212,14 +244,16 @@ Foam::tmp<Foam::tensorField> Foam::localAxesRotation::transformTensor
<< abort(FatalError);
}
const tensorField Rtr(Rptr_().T());
const tensorField& R = Rptr_();
const tensorField Rtr(R.T());
tmp<tensorField> tt(new tensorField(cellMap.size()));
tensorField& t = tt();
forAll(cellMap, i)
{
const label cellI = cellMap[i];
t[i] = Rptr_()[cellI] & st[i] & Rtr[cellI];
t[i] = R[cellI] & st[i] & Rtr[cellI];
}
return tt;
}
@ -239,9 +273,10 @@ Foam::tmp<Foam::symmTensorField> Foam::localAxesRotation::transformVector
tmp<symmTensorField> tfld(new symmTensorField(Rptr_->size()));
symmTensorField& fld = tfld();
const tensorField& R = Rptr_();
forAll(fld, i)
{
fld[i] = transformPrincipal(Rptr_()[i], st[i]);
fld[i] = transformPrincipal(R[i], st[i]);
}
return tfld;
}
@ -260,25 +295,6 @@ Foam::symmTensor Foam::localAxesRotation::transformVector
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::localAxesRotation::init
(
const dictionary& dict,
const objectRegistry& obr
)
{
const polyMesh& mesh = refCast<const polyMesh>(obr);
forAll(mesh.cellCentres(), cellI)
{
vector dir = mesh.cellCentres()[cellI] - origin_;
dir /= mag(dir) + VSMALL;
Rptr_()[cellI] = axesRotation(e3_, dir).R();
}
}
void Foam::localAxesRotation::write(Ostream& os) const
{
os.writeKeyword("e3") << e3() << token::END_STATEMENT << nl;

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -74,8 +74,8 @@ class localAxesRotation
// Private members
//- Init transformation tensor
void init(const dictionary& dict, const objectRegistry& obr);
//- Init transformation tensor field
void init(const objectRegistry& obr);
public:
@ -108,6 +108,9 @@ public:
//- Reset rotation to an identity rotation
virtual void clear();
//- Update the rotation for a list of cells
virtual void updateCells(const polyMesh& mesh, const labelList& cells);
//- Return local-to-global transformation tensor
virtual const tensor& R() const
{

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -127,7 +127,6 @@ Foam::coordinateSystem::coordinateSystem
origin_(point::zero),
R_()
{
const entry* entryPtr = dict.lookupEntryPtr(typeName_(), false, false);
// non-dictionary entry is a lookup into global coordinateSystems
@ -141,7 +140,7 @@ Foam::coordinateSystem::coordinateSystem
if (debug)
{
Info<< "coordinateSystem::coordinateSystem"
"(const dictionary&, const objectRegistry&):"
"(const objectRegistry&, const dictionary&):"
<< nl << "using global coordinate system: "
<< key << "=" << index << endl;
}
@ -151,7 +150,7 @@ Foam::coordinateSystem::coordinateSystem
FatalErrorIn
(
"coordinateSystem::coordinateSystem"
"(const dictionary&, const objectRegistry&)"
"(const objectRegistry&, const dictionary&):"
) << "could not find coordinate system: " << key << nl
<< "available coordinate systems: " << lst.toc() << nl << nl
<< exit(FatalError);
@ -344,7 +343,11 @@ void Foam::coordinateSystem::init
{
if (debug)
{
Pout<< "coordinateSystem::operator=(const dictionary&) : "
Pout<< "coordinateSystem::operator="
"("
"const dictionary&, "
"const objectRegistry&"
") : "
<< "assign from " << rhs << endl;
}

View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,10 +25,9 @@ Class
Foam::coordinateSystem
Description
Base class for other coordinate
system specifications.
Base class for other coordinate system specifications.
All systems are defined by an origin point and a coordinateRotation.
All systems are defined by an origin point and a co-ordinate rotation.
\verbatim
coordinateSystem
@ -37,8 +36,8 @@ Description
origin (0 0 0);
coordinateRotation
{
type STARCDRotation;
rotation (0 0 90);
type localAxesRotation;
e3 (0 0 1);
}
}
\endverbatim
@ -49,16 +48,17 @@ Description
3) localAxesRotation
4) EulerCoordinateRotation
Type of coordinates:
Type of co-ordinates:
1) cartesian
2) cylindricalCS
2) cylindrical
See Also
coordinateSystem and coordinateSystem::New
SourceFiles
coordinateSystem.C
newCoordinateSystem.C
coordinateSystemNew.C
\*---------------------------------------------------------------------------*/
#ifndef coordinateSystem_H
@ -231,7 +231,6 @@ public:
// Member Functions
// Access
//- Return name
@ -258,18 +257,28 @@ public:
return origin_;
}
//- Return const reference to coordinate rotation
//- Return const reference to co-ordinate rotation
const coordinateRotation& R() const
{
return R_();
}
//- Return non const reference to coordinate rotation
//- Return non const reference to co-ordinate rotation
coordinateRotation& R()
{
return R_();
}
//- Update and return the co-ordinate roation for a list of cells
const coordinateRotation& R
(
const polyMesh& mesh,
const labelList& cells
)
{
R_->updateCells(mesh, cells);
return R_();
}
//- Return as dictionary of entries
// \param [in] ignoreType drop type (cartesian, cylindrical, etc)