Commit 5dd93c11 authored by mattijs's avatar mattijs
Browse files

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

parents 992344ff 0f3e4160
......@@ -296,7 +296,7 @@ int main(int argc, char *argv[])
(
"density",
"scalar",
"Specify density,"
"Specify density, "
"kg/m3 for solid properties, kg/m2 for shell properties"
);
......
......@@ -128,7 +128,20 @@ Foam::plane::plane(const vector& normalVector)
:
unitVector_(normalVector),
basePoint_(vector::zero)
{}
{
scalar magUnitVector(mag(unitVector_));
if (magUnitVector > VSMALL)
{
unitVector_ /= magUnitVector;
}
else
{
FatalErrorIn("plane::plane(const point&, const vector&)")
<< "plane normal has zero length"
<< abort(FatalError);
}
}
// Construct from point and normal vector
......@@ -146,8 +159,8 @@ Foam::plane::plane(const point& basePoint, const vector& normalVector)
else
{
FatalErrorIn("plane::plane(const point&, const vector&)")
<< "plane normal has got zero length"
<< abort(FatalError);
<< "plane normal has zero length"
<< abort(FatalError);
}
}
......@@ -217,8 +230,8 @@ Foam::plane::plane(const dictionary& dict)
"plane::plane(const dictionary&)",
dict
)
<< "Invalid plane type: " << planeType
<< abort(FatalIOError);
<< "Invalid plane type: " << planeType
<< abort(FatalIOError);
}
}
......@@ -238,7 +251,7 @@ Foam::plane::plane(Istream& is)
else
{
FatalErrorIn("plane::plane(Istream& is)")
<< "plane normal has got zero length"
<< "plane normal has zero length"
<< abort(FatalError);
}
}
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "expDirectionMixedFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF
)
:
fvPatchField<Type>(p, iF),
refValue_(p.size()),
refGrad_(p.size()),
valueFraction_(p.size())
{}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>& ptf,
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fvPatchField<Type>(ptf, p, iF, mapper),
refValue_(ptf.refValue_, mapper),
refGrad_(ptf.refGrad_, mapper),
valueFraction_(ptf.valueFraction_, mapper)
{}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const fvPatch& p,
const DimensionedField<Type, volMesh>& iF,
const dictionary& dict
)
:
fvPatchField<Type>(p, iF, dict),
refValue_("refValue", dict, p.size()),
refGrad_("refGradient", dict, p.size()),
valueFraction_("valueFraction", dict, p.size())
{
evaluate();
}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>& ptf
)
:
fvPatchField<Type>(ptf),
refValue_(ptf.refValue_),
refGrad_(ptf.refGrad_),
valueFraction_(ptf.valueFraction_)
{}
template<class Type>
expDirectionMixedFvPatchField<Type>::expDirectionMixedFvPatchField
(
const expDirectionMixedFvPatchField<Type>& ptf,
const DimensionedField<Type, volMesh>& iF
)
:
fvPatchField<Type>(ptf, iF),
refValue_(ptf.refValue_),
refGrad_(ptf.refGrad_),
valueFraction_(ptf.valueFraction_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void expDirectionMixedFvPatchField<Type>::autoMap
(
const fvPatchFieldMapper& m
)
{
fvPatchField<Type>::autoMap(m);
refValue_.autoMap(m);
refGrad_.autoMap(m);
valueFraction_.autoMap(m);
}
// Reverse-map the given fvPatchField onto this fvPatchField
template<class Type>
void expDirectionMixedFvPatchField<Type>::rmap
(
const fvPatchField<Type>& ptf,
const labelList& addr
)
{
fvPatchField<Type>::rmap(ptf, addr);
const expDirectionMixedFvPatchField<Type>& edmptf =
refCast<const expDirectionMixedFvPatchField<Type> >(ptf);
refValue_.rmap(edmptf.refValue_, addr);
refGrad_.rmap(edmptf.refGrad_, addr);
valueFraction_.rmap(edmptf.valueFraction_, addr);
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::snGrad() const
{
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
return
valueFraction_*
(mixedValue - patchInternalField())*patch().deltaCoeffs()
+ (1.0 - valueFraction_)*refGrad_;
}
template<class Type>
void expDirectionMixedFvPatchField<Type>::evaluate(const Pstream::commsTypes)
{
if (!updated())
{
updateCoeffs();
}
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
Field<Type>::operator=
(
valueFraction_*mixedValue + (1.0 - valueFraction_)*gradValue
);
fvPatchField<Type>::evaluate();
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::valueInternalCoeffs
(
const tmp<scalarField>&
) const
{
return Type(pTraits<Type>::one)*(1.0 - valueFraction_);
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::valueBoundaryCoeffs
(
const tmp<scalarField>&
) const
{
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
return
valueFraction_*mixedValue
+ (1.0 - valueFraction_)*refGrad_/patch().deltaCoeffs();
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::
gradientInternalCoeffs() const
{
return -Type(pTraits<Type>::one)*valueFraction_*patch().deltaCoeffs();
}
template<class Type>
tmp<Field<Type> > expDirectionMixedFvPatchField<Type>::
gradientBoundaryCoeffs() const
{
const vectorField& nHat = patch().faceNormals();
Field<Type> gradValue =
patchInternalField() + refGrad_/patch().deltaCoeffs();
Field<Type> mixedValue =
nHat*(nHat & refValue_)
+ gradValue - nHat*(nHat & gradValue);
return
valueFraction_*patch().deltaCoeffs()*mixedValue
+ (1.0 - valueFraction_)*refGrad_;
}
template<class Type>
void expDirectionMixedFvPatchField<Type>::write(Ostream& os) const
{
fvPatchField<Type>::write(os);
refValue_.writeEntry("refValue", os);
refGrad_.writeEntry("refGradient", os);
valueFraction_.writeEntry("valueFraction", os);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2007-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
momentOfInertia
Description
Reimplementation of volInt.c by Brian Mirtich.
* mirtich@cs.berkeley.edu *
* http://www.cs.berkeley.edu/~mirtich *
-------------------------------------------------------------------------------
*/
#include "momentOfInertia.H"
//#include "pyramidPointFaceRef.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//Foam::tensor Foam::momentOfInertia
//(
// const pointField& points,
// const faceList& faces,
// const cell& cFaces,
// const point& cc
//)
//{
// tensor t(tensor::zero);
//
// forAll(cFaces, i)
// {
// const face& f = faces[cFaces[i]];
//
// scalar pyrVol = pyramidPointFaceRef(f, cc).mag(points);
//
// vector pyrCentre = pyramidPointFaceRef(f, cc).centre(points);
//
// vector d = pyrCentre - cc;
//
// t.xx() += pyrVol*(sqr(d.y()) + sqr(d.z()));
// t.yy() += pyrVol*(sqr(d.x()) + sqr(d.z()));
// t.zz() += pyrVol*(sqr(d.x()) + sqr(d.y()));
//
// t.xy() -= pyrVol*d.x()*d.y();
// t.xz() -= pyrVol*d.x()*d.z();
// t.yz() -= pyrVol*d.y()*d.z();
// }
//
// // Symmetric
// t.yx() = t.xy();
// t.zx() = t.xz();
// t.zy() = t.yz();
//
// return t;
//}
#define sqr(x) ((x)*(x))
#define pow3(x) ((x)*(x)*(x))
// compute various integrations over projection of face
void Foam::compProjectionIntegrals
(
const pointField& points,
const face& f,
const direction A,
const direction B,
scalar& P1,
scalar& Pa,
scalar& Pb,
scalar& Paa,
scalar& Pab,
scalar& Pbb,
scalar& Paaa,
scalar& Paab,
scalar& Pabb,
scalar& Pbbb
)
{
P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;
forAll(f, i)
{
scalar a0 = points[f[i]][A];
scalar b0 = points[f[i]][B];
scalar a1 = points[f[(i+1) % f.size()]][A];
scalar b1 = points[f[(i+1) % f.size()]][B];
scalar da = a1 - a0;
scalar db = b1 - b0;
scalar a0_2 = a0 * a0;
scalar a0_3 = a0_2 * a0;
scalar a0_4 = a0_3 * a0;
scalar b0_2 = b0 * b0;
scalar b0_3 = b0_2 * b0;
scalar b0_4 = b0_3 * b0;
scalar a1_2 = a1 * a1;
scalar a1_3 = a1_2 * a1;
scalar b1_2 = b1 * b1;
scalar b1_3 = b1_2 * b1;
scalar C1 = a1 + a0;
scalar Ca = a1*C1 + a0_2;
scalar Caa = a1*Ca + a0_3;
scalar Caaa = a1*Caa + a0_4;
scalar Cb = b1*(b1 + b0) + b0_2;
scalar Cbb = b1*Cb + b0_3;
scalar Cbbb = b1*Cbb + b0_4;
scalar Cab = 3*a1_2 + 2*a1*a0 + a0_2;
scalar Kab = a1_2 + 2*a1*a0 + 3*a0_2;
scalar Caab = a0*Cab + 4*a1_3;
scalar Kaab = a1*Kab + 4*a0_3;
scalar Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3;
scalar Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3;
P1 += db*C1;
Pa += db*Ca;
Paa += db*Caa;
Paaa += db*Caaa;
Pb += da*Cb;
Pbb += da*Cbb;
Pbbb += da*Cbbb;
Pab += db*(b1*Cab + b0*Kab);
Paab += db*(b1*Caab + b0*Kaab);
Pabb += da*(a1*Cabb + a0*Kabb);
}
P1 /= 2.0;
Pa /= 6.0;
Paa /= 12.0;
Paaa /= 20.0;
Pb /= -6.0;
Pbb /= -12.0;
Pbbb /= -20.0;
Pab /= 24.0;
Paab /= 60.0;
Pabb /= -60.0;
}
void Foam::compFaceIntegrals
(
const pointField& points,
const face& f,
const vector& n,
const scalar w,
const direction A,
const direction B,
const direction C,
scalar& Fa,
scalar& Fb,
scalar& Fc,
scalar& Faa,
scalar& Fbb,
scalar& Fcc,
scalar& Faaa,
scalar& Fbbb,
scalar& Fccc,
scalar& Faab,
scalar& Fbbc,
scalar& Fcca
)
{
scalar P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb;
compProjectionIntegrals
(
points,
f,
A,
B,
P1,
Pa,
Pb,
Paa,
Pab,
Pbb,
Paaa,
Paab,
Pabb,
Pbbb
);
scalar k1 = 1 / n[C];
scalar k2 = k1 * k1;
scalar k3 = k2 * k1;
scalar k4 = k3 * k1;
Fa = k1 * Pa;
Fb = k1 * Pb;
Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1);
Faa = k1 * Paa;
Fbb = k1 * Pbb;
Fcc = k3 * (sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb
+ w*(2*(n[A]*Pa + n[B]*Pb) + w*P1));
Faaa = k1 * Paaa;
Fbbb = k1 * Pbbb;
Fccc = -k4 * (pow3(n[A])*Paaa + 3*sqr(n[A])*n[B]*Paab
+ 3*n[A]*sqr(n[B])*Pabb + pow3(n[B])*Pbbb
+ 3*w*(sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb)
+ w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1));
Faab = k1 * Paab;
Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb);
Fcca = k3 * (sqr(n[A])*Paaa + 2*n[A]*n[B]*Paab + sqr(n[B])*Pabb
+ w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa));
}
void Foam::compVolumeIntegrals
(
const pointField& points,
const faceList& faces,
const cell& cFaces,
const vectorField& fNorm,
const scalarField& fW,
scalar& T0,
vector& T1,
vector& T2,
vector& TP
)
{
T0 = 0;
T1 = vector::zero;
T2 = vector::zero;
TP = vector::zero;