--- /dev/null
+Index: LESmodels/compressible/SpalartAllmaras/SpalartAllmaras.C
+===================================================================
+--- src/LESmodels/compressible/SpalartAllmaras/SpalartAllmaras.C (Revision 30)
++++ src/LESmodels/compressible/SpalartAllmaras/SpalartAllmaras.C (Revision 784)
+@@ -44,29 +44,38 @@
+
+
+ // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+ tmp<volScalarField> SpalartAllmaras::fv1() const
+ {
+ volScalarField chi3 = pow(nuTilda_/(mu()/rho()), 3);
+-
+ return chi3/(chi3 + pow(Cv1_, 3));
+ }
+
+
+ tmp<volScalarField> SpalartAllmaras::fv2() const
+ {
+ volScalarField chi = nuTilda_/(mu()/rho());
+ return scalar(1) - chi/(scalar(1) + chi*fv1());
+ }
+
+
+ tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& Stilda) const
+ {
+- volScalarField r = nuTilda_/(Stilda*sqr(kappa_*dTilda_));
++ volScalarField r = min
++ (
++ nuTilda_
++ /(
++ max(Stilda, dimensionedScalar("SMALL", Stilda.dimensions(), SMALL))
++ *sqr(kappa_*dTilda_)
++ ),
++ scalar(10.0)
++ );
++ r.boundaryField() == 0.0;
++
+ volScalarField g = r + Cw2_*(pow(r, 6) - r);
+
+ return g*pow((1.0 + pow(Cw3_, 6))/(pow(g, 6) + pow(Cw3_, 6)), 1.0/6.0);
+ }
+
+
+ // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Index: transportModels/interfaceProperties/gammaContactAngle/dynamicGammaContactAngle/dynamicGammaContactAngleFvPatchScalarField.C
+===================================================================
+--- src/transportModels/interfaceProperties/gammaContactAngle/dynamicGammaContactAngle/dynamicGammaContactAngleFvPatchScalarField.C (Revision 30)
++++ src/transportModels/interfaceProperties/gammaContactAngle/dynamicGammaContactAngle/dynamicGammaContactAngleFvPatchScalarField.C (Revision 784)
+@@ -60,15 +60,15 @@
+ const fvPatchFieldMapper& mapper
+ )
+ :
+ gammaContactAngleFvPatchScalarField(gcpsf, p, iF, mapper),
+ theta0_(gcpsf.theta0_),
+ uTheta_(gcpsf.uTheta_),
+ thetaA_(gcpsf.thetaA_),
+- thetaR_(gcpsf.thetaA_)
++ thetaR_(gcpsf.thetaR_)
+ {}
+
+
+ dynamicGammaContactAngleFvPatchScalarField::
+ dynamicGammaContactAngleFvPatchScalarField
+ (
+ const fvPatch& p,
+
+Index: dynamicMesh/meshCut/meshModifiers/multiDirRefinement/multiDirRefinement.C
+===================================================================
+--- src/dynamicMesh/meshCut/meshModifiers/multiDirRefinement/multiDirRefinement.C (Revision 30)
++++ src/dynamicMesh/meshCut/meshModifiers/multiDirRefinement/multiDirRefinement.C (Revision 784)
+@@ -343,16 +343,16 @@
+ }
+ }
+ }
+
+
+ hexRefiner.setRefinement(consistentCells, meshMod);
+
+- // Use inflation
+- autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, true, true);
++ // Change mesh, no inflation
++ autoPtr<mapPolyMesh> morphMapPtr = meshMod.changeMesh(mesh, false, true);
+ const mapPolyMesh& morphMap = morphMapPtr();
+
+ if (morphMap.hasMotionPoints())
+ {
+ mesh.movePoints(morphMap.preMotionPoints());
+ }
+
+@@ -365,47 +365,31 @@
+ {
+ Pout<< "multiDirRefinement : updated mesh at time "
+ << mesh.time().timeName() << endl;
+ }
+
+ hexRefiner.updateMesh(morphMap);
+
+- // Take over split pattern from hex refiner. (should be empty at this
+- // point)
+-
+- // From old cell label to index
+- Map<label> consistentSet(2*consistentCells.size());
+-
+- forAll(consistentCells, i)
+- {
+- consistentSet.insert(consistentCells[i], i);
+- }
+-
+ // Collect all cells originating from same old cell (original + 7 extra)
+
+- addedCells_.setSize(consistentCells.size());
+- forAll(addedCells_, i)
++ forAll(consistentCells, i)
+ {
+- addedCells_[i].setSize(8);
++ addedCells_[consistentCells[i]].setSize(8);
+ }
+- labelList nAddedCells(consistentCells.size(), 0);
++ labelList nAddedCells(addedCells_.size(), 0);
+
+ const labelList& cellMap = morphMap.cellMap();
+
+ forAll(cellMap, cellI)
+ {
+ label oldCellI = cellMap[cellI];
+
+- Map<label>::const_iterator iter = consistentSet.find(oldCellI);
+-
+- if (iter != consistentSet.end())
++ if (addedCells_[oldCellI].size() > 0)
+ {
+- label index = iter();
+-
+- addedCells_[nAddedCells[index]++] = cellI;
++ addedCells_[oldCellI][nAddedCells[oldCellI]++] = cellI;
+ }
+ }
+ }
+
+
+ void Foam::multiDirRefinement::refineAllDirs
+ (
+
+Index: OpenFOAM/db/error/error.H
+===================================================================
+--- src/OpenFOAM/db/error/error.H (Revision 30)
++++ src/OpenFOAM/db/error/error.H (Revision 784)
+@@ -80,14 +80,17 @@
+
+ //- Construct from title string
+ error(const string& title);
+
+ //- Construct from dictionary
+ error(const dictionary& errDict);
+
++ //- Construct as copy
++ error(const error& err);
++
+
+ // Destructor
+
+ ~error();
+
+
+ // Member functions
+Index: OpenFOAM/db/error/error.C
+===================================================================
+--- src/OpenFOAM/db/error/error.C (Revision 30)
++++ src/OpenFOAM/db/error/error.C (Revision 784)
+@@ -76,14 +76,28 @@
+ "cannot open error stream"
+ << endl;
+ ::exit(1);
+ }
+ }
+
+
++error::error(const error& err)
++:
++ messageStream(err),
++ functionName_(err.functionName_),
++ sourceFileName_(err.sourceFileName_),
++ sourceFileLineNumber_(err.sourceFileLineNumber_),
++ abort_(err.abort_),
++ throwExceptions_(err.throwExceptions_),
++ messageStreamPtr_(new OStringStream(*err.messageStreamPtr_))
++{
++ //*messageStreamPtr_ << err.message();
++}
++
++
+ error::~error()
+ {
+ delete messageStreamPtr_;
+ }
+
+
+ OSstream& error::operator()
+@@ -251,14 +265,15 @@
+ return os;
+ }
+
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+ // Global error definitions
+
+-error FatalError ("--> FOAM FATAL ERROR : ");
++error FatalError("--> FOAM FATAL ERROR : ");
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ } // End namespace Foam
+
+ // ************************************************************************* //
++
+Index: OpenFOAM/db/Time/findTimes.C
+===================================================================
+--- src/OpenFOAM/db/Time/findTimes.C (Revision 30)
++++ src/OpenFOAM/db/Time/findTimes.C (Revision 784)
+@@ -52,21 +52,23 @@
+ fileNameList dirEntries(readDir(directory, fileName::DIRECTORY));
+
+ // Initialise instant list
+ instantList Times(dirEntries.size() + 1);
+ label nTimes = 0;
+
+ // Check for "constant"
++ bool haveConstant = false;
+ forAll(dirEntries, i)
+ {
+ if (dirEntries[i] == "constant")
+ {
+ Times[nTimes].value() = 0;
+ Times[nTimes].name() = dirEntries[i];
+ nTimes++;
++ haveConstant = true;
+ break;
+ }
+ }
+
+ // Read and parse all the entries in the directory
+ forAll(dirEntries, i)
+ {
+@@ -80,17 +82,24 @@
+ nTimes++;
+ }
+ }
+
+ // Reset the length of the times list
+ Times.setSize(nTimes);
+
+- if (nTimes > 1)
++ if (haveConstant)
+ {
+- std::sort(&Times[1], Times.end(), instant::less());
++ if (nTimes > 2)
++ {
++ std::sort(&Times[1], Times.end(), instant::less());
++ }
++ }
++ else if (nTimes > 1)
++ {
++ std::sort(&Times[0], Times.end(), instant::less());
+ }
+
+ return Times;
+ }
+
+
+ // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+Index: OpenFOAM/db/dictionary/dictionary.C
+===================================================================
+--- src/OpenFOAM/db/dictionary/dictionary.C (Revision 30)
++++ src/OpenFOAM/db/dictionary/dictionary.C (Revision 784)
+@@ -119,14 +119,20 @@
+ }
+ else
+ {
+ return -1;
+ }
+ }
+
++// Clear the Dictionary
++void Foam::dictionary::clear()
++{
++ IDLList<entry>::clear();
++ hashedEntries_.clear();
++}
+
+ // Find and return entry
+ bool Foam::dictionary::found(const word& keyword) const
+ {
+ return hashedEntries_.found(keyword);
+ }
+
+Index: OpenFOAM/db/dictionary/dictionary.H
+===================================================================
+--- src/OpenFOAM/db/dictionary/dictionary.H (Revision 30)
++++ src/OpenFOAM/db/dictionary/dictionary.H (Revision 784)
+@@ -126,14 +126,16 @@
+
+ //- Return line number of first token in dictionary
+ label startLineNumber() const;
+
+ //- Return line number of last token in dictionary
+ label endLineNumber() const;
+
++ //- Clear the dictionary
++ void clear();
+
+ // Search and lookup
+
+ //- Search dictionary for given keyword
+ bool found(const word& keyword) const;
+
+ //- Find and return an entry data stream
+Index: OpenFOAM/db/IOstreams/StringStreams/OStringStream.H
+===================================================================
+--- src/OpenFOAM/db/IOstreams/StringStreams/OStringStream.H (Revision 30)
++++ src/OpenFOAM/db/IOstreams/StringStreams/OStringStream.H (Revision 784)
+@@ -63,21 +63,41 @@
+ (
+ streamFormat format=ASCII,
+ versionNumber version=currentVersion
+ )
+ :
+ OSstream
+ (
+- *(new std::ostringstream()),
++ *(new std::ostringstream()),
+ "OStringStream.sinkFile",
+ format,
+ version
+ )
+ {}
+
++ //- Construct as copy
++ OStringStream(const OStringStream& oss)
++ :
++ OSstream
++ (
++ *(
++ new std::ostringstream
++ (
++ dynamic_cast<const std::ostringstream&>
++ (
++ oss.stream()
++ ).str()
++ )
++ ),
++ oss.name(),
++ oss.format(),
++ oss.version()
++ )
++ {}
++
+
+ // Destructor
+
+ ~OStringStream()
+ {
+ delete &dynamic_cast<std::ostringstream&>(stream());
+ }
+Index: OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldFunctions.H
+===================================================================
+--- src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldFunctions.H (Revision 30)
++++ src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldFunctions.H (Revision 784)
+@@ -214,32 +214,28 @@
+ #undef BINARY_TYPE_OPERATOR_FT
+
+
+ #define PRODUCT_OPERATOR(product, op, opFunc) \
+ \
+ template \
+ < \
+- template<class> class PatchField, \
+- class pointPatch, \
+ class Type1, \
+ class Type2 \
+ > \
+ inline void opFunc \
+ ( \
+ pointPatchField \
+ <typename product<Type1, Type2>::type>& f, \
+ const pointPatchField<Type1>& f1, \
+ const pointPatchField<Type2>& f2 \
+ ) \
+ {} \
+ \
+ template \
+ < \
+- template<class> class PatchField, \
+- class pointPatch, \
+ class Type, \
+ class Form, \
+ class Cmpt, \
+ int nCmpt \
+ > \
+ inline void opFunc \
+ ( \
+@@ -248,16 +244,14 @@
+ const pointPatchField<Type>& f1, \
+ const VectorSpace<Form,Cmpt,nCmpt>& vs \
+ ) \
+ {} \
+ \
+ template \
+ < \
+- template<class> class PatchField, \
+- class pointPatch, \
+ class Form, \
+ class Cmpt, \
+ int nCmpt, \
+ class Type \
+ > \
+ inline void opFunc \
+ ( \
+@@ -272,112 +266,99 @@
+ PRODUCT_OPERATOR(crossProduct, ^, cross)
+ PRODUCT_OPERATOR(innerProduct, &, dot)
+ PRODUCT_OPERATOR(scalarProduct, &&, dotdot)
+
+ #undef PRODUCT_OPERATOR
+
+
+-template<template<class> class PatchField, class pointPatch>
+-void hdual
++inline void hdual
+ (
+ pointPatchField<vector>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void hdual
++inline void hdual
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<vector>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void diag
++inline void diag
+ (
+ pointPatchField<vector>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void tr
++inline void tr
+ (
+ pointPatchField<scalar>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void dev
++inline void dev
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void dev2
++inline void dev2
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void det
++inline void det
+ (
+ pointPatchField<scalar>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void inv
++inline void inv
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void symm
++inline void symm
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void twoSymm
++inline void twoSymm
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void skew
++inline void skew
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void eigenValues
++inline void eigenValues
+ (
+ pointPatchField<vector>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+-template<template<class> class PatchField, class pointPatch>
+-void eigenVectors
++inline void eigenVectors
+ (
+ pointPatchField<tensor>&,
+ const pointPatchField<tensor>&
+ )
+ {}
+
+
+Index: OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
+===================================================================
+--- src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C (Revision 30)
++++ src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C (Revision 784)
+@@ -131,16 +131,15 @@
+
+ template<>
+ tmp<Field<symmTensor> > transformFieldMask<symmTensor>
+ (
+ const tmp<tensorField>& ttf
+ )
+ {
+- tmp<Field<symmTensor> > ret =
+- transformFieldMask<symmTensor>(ttf());
++ tmp<Field<symmTensor> > ret = transformFieldMask<symmTensor>(ttf());
+ ttf.clear();
+ return ret;
+ }
+
+
+ template<>
+ tmp<Field<symmTensor> > transformFieldMask<symmTensor>
+@@ -153,18 +152,15 @@
+
+ template<>
+ tmp<Field<symmTensor> > transformFieldMask<symmTensor>
+ (
+ const tmp<symmTensorField>& tstf
+ )
+ {
+- tmp<Field<symmTensor> > ret =
+- transformFieldMask<symmTensor>(tstf());
+- tstf.clear();
+- return ret;
++ return tstf;
+ }
+
+
+ // * * * * * * * * * * * * * * * global operators * * * * * * * * * * * * * //
+
+ UNARY_OPERATOR(vector, symmTensor, *, hdual)
+
+Index: OpenFOAM/primitives/pTraits/pTraits.H
+===================================================================
+--- src/OpenFOAM/primitives/pTraits/pTraits.H (Revision 30)
++++ src/OpenFOAM/primitives/pTraits/pTraits.H (Revision 784)
+@@ -19,15 +19,15 @@
+ 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
+- pTraits
++ Foam::pTraits
+
+ Description
+ Traits class for primitives. All primitives need a specialised
+ version of this class
+
+ \*---------------------------------------------------------------------------*/
+
+@@ -49,14 +49,19 @@
+ class pTraits
+ :
+ public primitive
+ {
+
+ public:
+
++ pTraits(const primitive& p)
++ :
++ primitive(p)
++ {}
++
+ pTraits(Istream& is)
+ :
+ primitive(is)
+ {}
+ };
+
+
+Index: OpenFOAM/primitives/Tensor/TensorI.H
+===================================================================
+--- src/OpenFOAM/primitives/Tensor/TensorI.H (Revision 30)
++++ src/OpenFOAM/primitives/Tensor/TensorI.H (Revision 784)
+@@ -358,14 +358,23 @@
+ v1.x()*v2.x(), v1.x()*v2.y(), v1.x()*v2.z(),
+ v1.y()*v2.x(), v1.y()*v2.y(), v1.y()*v2.z(),
+ v1.z()*v2.x(), v1.z()*v2.y(), v1.z()*v2.z()
+ );
+ }
+
+
++//- Division of a vector by a tensor, i.e. dot-product with the tensor inverse
++template <class Cmpt>
++inline typename innerProduct<Vector<Cmpt>, Tensor<Cmpt> >::type
++operator/(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
++{
++ return inv(t) & v;
++}
++
++
+ // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
+
+ //- Return the trace of a tensor
+ template <class Cmpt>
+ inline Cmpt tr(const Tensor<Cmpt>& t)
+ {
+ return t.xx() + t.yy() + t.zz();
+@@ -707,15 +716,15 @@
+ template <class Cmpt>
+ inline Tensor<Cmpt>
+ operator+(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
+ {
+ return Tensor<Cmpt>
+ (
+ t1.xx() + st2.xx(), t1.xy() + st2.xy(), t1.xz() + st2.xz(),
+- t1.yx() + st2.xy(), t1.yy() + st2.yy(), t1.yz() + st2.yy(),
++ t1.yx() + st2.xy(), t1.yy() + st2.yy(), t1.yz() + st2.yz(),
+ t1.zx() + st2.xz(), t1.zy() + st2.yz(), t1.zz() + st2.zz()
+ );
+ }
+
+
+ template <class Cmpt>
+ inline Tensor<Cmpt>
+@@ -733,15 +742,15 @@
+ template <class Cmpt>
+ inline Tensor<Cmpt>
+ operator-(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
+ {
+ return Tensor<Cmpt>
+ (
+ t1.xx() - st2.xx(), t1.xy() - st2.xy(), t1.xz() - st2.xz(),
+- t1.yx() - st2.xy(), t1.yy() - st2.yy(), t1.yz() - st2.yy(),
++ t1.yx() - st2.xy(), t1.yy() - st2.yy(), t1.yz() - st2.yz(),
+ t1.zx() - st2.xz(), t1.zy() - st2.yz(), t1.zz() - st2.zz()
+ );
+ }
+
+
+ //- Inner-product between a spherical tensor and a tensor
+ template <class Cmpt>
+@@ -790,46 +799,30 @@
+ //- Double-dot-product between a spherical tensor and a tensor
+ template <class Cmpt>
+ inline Cmpt
+ operator&&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
+ {
+ return
+ (
+- st1.xx()*t2.xx() + st1.xy()*t2.yx() + st1.xz()*t2.zx() +
+- st1.xx()*t2.xy() + st1.xy()*t2.yy() + st1.xz()*t2.zy() +
+- st1.xx()*t2.xz() + st1.xy()*t2.yz() + st1.xz()*t2.zz() +
+-
+- st1.xy()*t2.xx() + st1.yy()*t2.yx() + st1.yz()*t2.zx() +
+- st1.xy()*t2.xy() + st1.yy()*t2.yy() + st1.yz()*t2.zy() +
+- st1.xy()*t2.xz() + st1.yy()*t2.yz() + st1.yz()*t2.zz() +
+-
+- st1.xz()*t2.xx() + st1.yz()*t2.yx() + st1.zz()*t2.zx() +
+- st1.xz()*t2.xy() + st1.yz()*t2.yy() + st1.zz()*t2.zy() +
+- st1.xz()*t2.xz() + st1.yz()*t2.yz() + st1.zz()*t2.zz()
++ st1.xx()*t2.xx() + st1.xy()*t2.xy() + st1.xz()*t2.xz() +
++ st1.xy()*t2.yx() + st1.yy()*t2.yy() + st1.yz()*t2.yz() +
++ st1.xz()*t2.zx() + st1.yz()*t2.zy() + st1.zz()*t2.zz()
+ );
+ }
+
+
+ //- Double-dot-product between a tensor and a spherical tensor
+ template <class Cmpt>
+ inline Cmpt
+ operator&&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
+ {
+ return
+ (
+ t1.xx()*st2.xx() + t1.xy()*st2.xy() + t1.xz()*st2.xz() +
+- t1.xx()*st2.xy() + t1.xy()*st2.yy() + t1.xz()*st2.yz() +
+- t1.xx()*st2.xz() + t1.xy()*st2.yz() + t1.xz()*st2.zz() +
+-
+- t1.yx()*st2.xx() + t1.yy()*st2.xy() + t1.yz()*st2.xz() +
+ t1.yx()*st2.xy() + t1.yy()*st2.yy() + t1.yz()*st2.yz() +
+- t1.yx()*st2.xz() + t1.yy()*st2.yz() + t1.yz()*st2.zz() +
+-
+- t1.zx()*st2.xx() + t1.zy()*st2.xy() + t1.zz()*st2.xz() +
+- t1.zx()*st2.xy() + t1.zy()*st2.yy() + t1.zz()*st2.yz() +
+ t1.zx()*st2.xz() + t1.zy()*st2.yz() + t1.zz()*st2.zz()
+ );
+ }
+
+
+ template<class Cmpt>
+ class typeOfSum<SymmTensor<Cmpt>, Tensor<Cmpt> >
+
+Index: finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C
+===================================================================
+--- src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C (Revision 30)
++++ src/finiteVolume/fvMatrices/fvMatrix/fvMatrixSolve.C (Revision 784)
+@@ -20,23 +20,18 @@
+
+ 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
+
+ \*---------------------------------------------------------------------------*/
+
+-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+-
+-namespace Foam
+-{
+-
+ // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+ template<class Type>
+-void fvMatrix<Type>::setComponentReference
++void Foam::fvMatrix<Type>::setComponentReference
+ (
+ const label patchi,
+ const label facei,
+ const direction cmpt,
+ const scalar value
+ )
+ {
+@@ -52,15 +47,18 @@
+ *value;
+ }
+ }
+ }
+
+
+ template<class Type>
+-lduMatrix::solverPerformance fvMatrix<Type>::solve(Istream& solverControls)
++Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve
++(
++ Istream& solverControls
++)
+ {
+ if (debug)
+ {
+ Info<< "fvMatrix<Type>::solve(Istream& solverControls) : "
+ "solving fvMatrix<Type>"
+ << endl;
+ }
+@@ -70,15 +68,19 @@
+ "fvMatrix<Type>::solve",
+ psi_.name()
+ );
+
+ scalarField saveDiag = diag();
+
+ Field<Type> source = source_;
+- addBoundarySource(source, false);
++
++ // At this point include the boundary source from the coupled boundaries.
++ // This is corrected for the implict part by updateMatrixInterfaces within
++ // the component loop.
++ addBoundarySource(source);
+
+ typename Type::labelType validComponents
+ (
+ pow
+ (
+ psi_.mesh().directions(),
+ pTraits<typename powProduct<Vector<label>, Type::rank>::type>::zero
+@@ -105,14 +107,35 @@
+ (
+ internalCoeffs_.component(cmpt)
+ );
+
+ lduInterfaceFieldPtrsList interfaces =
+ psi_.boundaryField().interfaces();
+
++ // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
++ // bouCoeffsCmpt for the explicit part of the coupled boundary
++ // conditions
++ initMatrixInterfaces
++ (
++ bouCoeffsCmpt,
++ interfaces,
++ psiCmpt,
++ sourceCmpt,
++ cmpt
++ );
++
++ updateMatrixInterfaces
++ (
++ bouCoeffsCmpt,
++ interfaces,
++ psiCmpt,
++ sourceCmpt,
++ cmpt
++ );
++
+ lduMatrix::solverPerformance solverPerf;
+
+ // Solver call
+ solverPerf = lduMatrix::solver::New
+ (
+ psi_.name() + pTraits<Type>::componentNames[cmpt],
+ *this,
+@@ -140,35 +163,36 @@
+ psi_.correctBoundaryConditions();
+
+ return solverPerfVec;
+ }
+
+
+ template<class Type>
+-autoPtr<typename fvMatrix<Type>::fvSolver> fvMatrix<Type>::solver()
++Foam::autoPtr<typename Foam::fvMatrix<Type>::fvSolver>
++Foam::fvMatrix<Type>::solver()
+ {
+ return solver(psi_.mesh().solver(psi_.name()));
+ }
+
+ template<class Type>
+-lduMatrix::solverPerformance fvMatrix<Type>::fvSolver::solve()
++Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::fvSolver::solve()
+ {
+ return solve(psi_.mesh().solver(psi_.name()));
+ }
+
+
+ template<class Type>
+-lduMatrix::solverPerformance fvMatrix<Type>::solve()
++Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
+ {
+ return solve(psi_.mesh().solver(psi_.name()));
+ }
+
+
+ template<class Type>
+-tmp<Field<Type> > fvMatrix<Type>::residual() const
++Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::residual() const
+ {
+ tmp<Field<Type> > tres(source_);
+ Field<Type>& res = tres();
+
+ addBoundarySource(res);
+
+ // Loop over field components
+@@ -198,12 +222,8 @@
+ );
+ }
+
+ return tres;
+ }
+
+
+-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+-
+-} // End namespace Foam
+-
+ // ************************************************************************* //
+Index: finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
+===================================================================
+--- src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C (Revision 30)
++++ src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C (Revision 784)
+@@ -181,14 +181,15 @@
+ )
+ :
+ fixedValueFvPatchField<Type>(ptf, iF),
+ setAverage_(ptf.setAverage_),
+ referenceCS_(ptf.referenceCS_),
+ nearestVertex_(ptf.nearestVertex_),
+ nearestVertexWeight_(ptf.nearestVertexWeight_),
++ fieldName_(ptf.fieldName_),
+ sampleTimes_(ptf.sampleTimes_),
+ startSampleTime_(ptf.startSampleTime_),
+ startSampledValues_(ptf.startSampledValues_),
+ startAverage_(ptf.startAverage_),
+ endSampleTime_(ptf.endSampleTime_),
+ endSampledValues_(ptf.endSampledValues_),
+ endAverage_(ptf.endAverage_)
+Index: finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C
+===================================================================
+--- src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C (Revision 30)
++++ src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C (Revision 784)
+@@ -323,15 +323,15 @@
+ *vf.oldTime().internalField()*mesh().V0()
+ - coefft00*rho.oldTime().oldTime().internalField()
+ *vf.oldTime().oldTime().internalField()*mesh().V00()
+ )/mesh().V()
+ ),
+ rDeltaT.value()*
+ (
+- coefft*vf.boundaryField() -
++ coefft*rho.boundaryField()*vf.boundaryField() -
+ (
+ coefft0*rho.oldTime().boundaryField()
+ *vf.oldTime().boundaryField()
+ - coefft00*rho.oldTime().oldTime().boundaryField()
+ *vf.oldTime().oldTime().boundaryField()
+ )
+ )
+Index: finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
+===================================================================
+--- src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C (Revision 30)
++++ src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C (Revision 784)
+@@ -150,15 +150,15 @@
+ if (p.coupled())
+ {
+ forAll(pd, patchFacei)
+ {
+ const vector& d = pd[patchFacei];
+
+ dd[faceCells[patchFacei]] +=
+- (pw[patchFacei]*pMagSf[patchFacei]/magSqr(d))*sqr(d);
++ ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))*sqr(d);
+ }
+ }
+ else
+ {
+ forAll(pd, patchFacei)
+ {
+ const vector& d = pd[patchFacei];
+
+Index: turbulenceModels/compressible/kOmegaSST/kOmegaSST.C
+===================================================================
+--- src/turbulenceModels/compressible/kOmegaSST/kOmegaSST.C (Revision 30)
++++ src/turbulenceModels/compressible/kOmegaSST/kOmegaSST.C (Revision 784)
+@@ -237,15 +237,15 @@
+ {
+ y_.correct();
+ divU += fvc::div(mesh_.phi());
+ }
+
+ tmp<volTensorField> tgradU = fvc::grad(U_);
+ volScalarField S2 = magSqr(symm(tgradU()));
+- volScalarField GbyMu = 2*mut_*(tgradU() && dev(symm(tgradU())));
++ volScalarField GbyMu = 2*(tgradU() && dev(symm(tgradU())));
+ volScalarField G = mut_*GbyMu;
+ tgradU.clear();
+
+ # include "kOmegaWallFunctionsI.H"
+
+ volScalarField CDkOmega =
+ (2*alphaOmega2)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
+@@ -280,15 +280,15 @@
+ // Turbulent kinetic energy equation
+ tmp<fvScalarMatrix> kEqn
+ (
+ fvm::ddt(rho_, k_)
+ + fvm::div(phi_, k_)
+ - fvm::laplacian(DkEff(F1), k_)
+ ==
+- min(G, c1*betaStar*k_*omega_)
++ min(G, (c1*betaStar)*rho_*k_*omega_)
+ - fvm::SuSp(2.0/3.0*rho_*divU, k_)
+ - fvm::Sp(rho_*betaStar*omega_, k_)
+ );
+
+ kEqn().relax();
+ solve(kEqn);
+ bound(k_, k0_);
+
projects / ~andy / sunrise / commitdiff