Block solver development

applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C

@@ -95,25 +95,28 @@ int main(int argc, char *argv[])
             // Prepare block system
             BlockLduMatrix<vector2> blockM(mesh);
 
+            //- Transfer the coupled interface list for processor/cyclic/etc.
+            // boundaries
+            blockM.interfaces() = blockT.boundaryField().blockInterfaces();
+
             // Grab block diagonal and set it to zero
             Field<tensor2>& d = blockM.diag().asSquare();
             d = tensor2::zero;
 
             // Grab linear off-diagonal and set it to zero
-            Field<vector2>& u = blockM.upper().asLinear();
             Field<vector2>& l = blockM.lower().asLinear();
+            Field<vector2>& u = blockM.upper().asLinear();
             u = vector2::zero;
             l = vector2::zero;
 
             vector2Field& blockX = blockT.internalField();
-            // vector2Field blockX(mesh.nCells(), vector2::zero);
             vector2Field blockB(mesh.nCells(), vector2::zero);
 
             //- Inset equations into block Matrix
             blockMatrixTools::insertEquation(0, TEqn, blockM, blockX, blockB);
             blockMatrixTools::insertEquation(1, TsEqn, blockM, blockX, blockB);
 
-            //- Add off-diagonal terms and remove from Block source
+            //- Add off-diagonal terms and remove from block source
             forAll(d, i)
             {
                 d[i](0,1) = -alpha.value()*mesh.V()[i];
@@ -123,44 +126,13 @@ int main(int argc, char *argv[])
                 blockB[i][1] -= alpha.value()*blockX[i][0]*mesh.V()[i];
             }
 
-            //- Transfer the coupled interface list for processor/cyclic/etc.
-            // boundaries
-            blockM.interfaces() = blockT.boundaryField().blockInterfaces();
-
-            //- Transfer the coupled interface coefficients
-            forAll(mesh.boundaryMesh(), patchI)
-            {
-                if (blockM.interfaces().set(patchI))
-                {
-                    Field<vector2>& coupledLower =
-                        blockM.coupleLower()[patchI].asLinear();
-
-                    Field<vector2>& coupledUpper =
-                        blockM.coupleUpper()[patchI].asLinear();
-
-                    const scalarField& TLower = TEqn.internalCoeffs()[patchI];
-                    const scalarField& TUpper = TEqn.boundaryCoeffs()[patchI];
-
-                    const scalarField& TsLower =
-                        TsEqn.internalCoeffs()[patchI];
-
-                    const scalarField& TsUpper =
-                        TsEqn.boundaryCoeffs()[patchI];
-
-                    blockMatrixTools::blockInsert(0, TLower, coupledLower);
-                    blockMatrixTools::blockInsert(1, TsLower, coupledLower);
-                    blockMatrixTools::blockInsert(0, TUpper, coupledUpper);
-                    blockMatrixTools::blockInsert(1, TsUpper, coupledUpper);
-                }
-            }
-
             //- Block coupled solver call
             BlockSolverPerformance<vector2> solverPerf =
                 BlockLduSolver<vector2>::New
                 (
-                    word("blockVar"),
+                    blockT.name(),
                     blockM,
-                    mesh.solutionDict().solver("blockVar")
+                    mesh.solutionDict().solver(blockT.name())
                 )->solve(blockX, blockB);
 
             solverPerf.print();

applications/solvers/coupled/blockCoupledScalarTransportFoam/blockMatrixTools.H

@@ -1,130 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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
-
-InNamespace
-    Foam::
-
-Description
-    Block matrix insertion and retrieval tools
-
-SourceFiles
-    blockMatrixTools.C
-
-\*---------------------------------------------------------------------------*/
-
-
-#ifndef blockMatrixTools_H
-#define blockMatrixTools_H
-
-#include "blockLduMatrices.H"
-#include "fvMatrices.H"
-#include "surfaceFieldsFwd.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-              Namespace blockMatrixTools functions Declaration
-\*---------------------------------------------------------------------------*/
-
-namespace blockMatrixTools
-{
-    //- Insert field into block field
-    template<class BlockType>
-    void blockInsert
-    (
-        const direction dir,
-        const scalarField& x,
-        Field<BlockType>& blockX
-    );
-
-    //- Retrieve field from block field
-    template<class BlockType>
-    void blockRetrieve
-    (
-        const direction dir,
-        scalarField& x,
-        const Field<BlockType>& blockX
-    );
-
-    //- Insert matrix diagonal and source into the block system
-    template<class BlockType>
-    void insertDiagSource
-    (
-        const direction dir,
-        const fvScalarMatrix& m,
-        BlockLduMatrix<BlockType>& blockM,
-        Field<BlockType>& blockB
-    );
-
-    // Insert matrix into the block system
-    template<class Type, class BlockType>
-    void insertUpperLower
-    (
-        const direction dir,
-        const fvScalarMatrix& m,
-        BlockLduMatrix<BlockType>& blockM
-    );
-
-    // Insert matrix into the block system
-    template<class BlockType>
-    void insertEquation
-    (
-        const direction dir,
-        const fvScalarMatrix& m,
-        BlockLduMatrix<BlockType>& blockM,
-        Field<BlockType>& blockX,
-        Field<BlockType>& blockB
-    );
-
-    // Update coupling of block system
-    // Subtracts the block-coefficient coupling as specified by the user
-    // from the source, leaving the implicit update given by linearisation
-    template<class BlockType>
-    void updateCoupling
-    (
-        BlockLduMatrix<BlockType>& blockM,
-        Field<BlockType>& x,
-        Field<BlockType>& b
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "blockMatrixTools.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/coupled/blockCoupledScalarTransportFoam/createFields.H

@@ -41,7 +41,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedVector2(word(), dimless, vector2::zero)
+        dimensionedVector2("zero", dimless, vector2::zero)
     );
 
     {
@@ -53,38 +53,6 @@
 
     blockT.write();
 
-// Removed comparison with the reference field.  HJ, 17/Jun/2010
-#if (0)
-    Info<< "Reading field Tref\n" << endl;
-    volScalarField Tref
-    (
-        IOobject
-        (
-            "Tref",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Reading field Tsref\n" << endl;
-    volScalarField Tsref
-    (
-        IOobject
-        (
-            "Tsref",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-#endif
-
     Info<< "Reading field U\n" << endl;
 
     volVectorField U

src/OpenFOAM/Make/files

@@ -642,11 +642,11 @@ matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/blockILUSmoothers.C
 
 matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/blockLduSolvers.C
 matrices/blockLduMatrix/BlockLduSolvers/BlockDiagonal/blockDiagonalSolvers.C
-/*HJ matrices/blockLduMatrix/BlockLduSolvers/Segregated/segregatedSolvers.C */
 matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/blockGaussSeidelSolvers.C
 matrices/blockLduMatrix/BlockLduSolvers/BlockCG/blockCGSolvers.C
 matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/blockBiCGStabSolvers.C
 matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/blockGMRESSolvers.C
+matrices/blockLduMatrix/BlockLduSolvers/Segregated/segregatedSolvers.C
 
 matrices/blockLduMatrix/BlockLduMatrix/BlockLduInterface/BlockLduInterfaceFields.C
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C

@@ -343,7 +343,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const BlockLduMatrix<Type>& ldum)
 {
     if (ldum.diagPtr_)
     {
-        os.writeKeyword("diagonal") << *ldum.diagPtr_ <<  token::END_STATEMENT << nl;
+        os.writeKeyword("diagonal")
+            << *ldum.diagPtr_ <<  token::END_STATEMENT << nl;
     }
     else
     {

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C

@@ -93,7 +93,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance()))
+    if (!stop(solverPerf))
     {
         scalar rho = this->great_;
         scalar rhoOld = rho;

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C

@@ -92,7 +92,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance()))
+    if (!stop(solverPerf))
     {
         scalar rho = this->great_;
         scalar rhoOld = rho;

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C

@@ -128,7 +128,7 @@ Foam::BlockGMRESSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance()))
+    if (!stop(solverPerf))
     {
         // Create the Hesenberg matrix
         scalarSquareMatrix H(nDirs_, 0);

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C

@@ -49,7 +49,7 @@ Foam::BlockGaussSeidelSolver<Type>::BlockGaussSeidelSolver
         dict
     ),
     gs_(matrix),
-    nResSweeps_(readLabel(this->dict().lookup("nResSweeps")))
+    nSweeps_(readLabel(this->dict().lookup("nSweeps")))
 {}
 
 
@@ -86,13 +86,13 @@ Foam::BlockGaussSeidelSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!solverPerf.checkConvergence(this->tolerance(), this->relTolerance()))
+    if (!stop(solverPerf))
     {
         // Iteration loop
 
         do
         {
-            for (label i = 0; i < nResSweeps_; i++)
+            for (label i = 0; i < nSweeps_; i++)
             {
                 gs_.precondition(x, b);
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.H

@@ -64,7 +64,7 @@ class BlockGaussSeidelSolver
         BlockGaussSeidelPrecon<Type> gs_;
 
         //- Number of sweeps before evaluating residual
-        label nResSweeps_;
+        label nSweeps_;
 
 
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C

@@ -82,7 +82,8 @@ Foam::scalar Foam::BlockIterativeSolver<Type>::normFactor
 
     if (BlockLduMatrix<Type>::debug >= 2)
     {
-        Info<< "Iterative solver normalisation factor = " << normFactor << endl;
+        Info<< "Iterative solver normalisation factor = "
+            << normFactor << endl;
     }
 
     return normFactor;

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.C

@@ -65,7 +65,7 @@ void Foam::BlockLduSolver<Type>::readControl
             ")",
             dict
         )   << "Cannot read control " << controlName
-            << abort(FatalIOError);
+            << exit(FatalIOError);
     }
 }
 
@@ -112,6 +112,52 @@ Foam::BlockLduSolver<Type>::New
 {
     word solverName(dict.lookup("solver"));
 
+    if (matrix.diagonal())
+    {
+        return autoPtr<BlockLduSolver<Type> >
+        (
+            new BlockDiagonalSolver<Type>(fieldName, matrix, dict)
+        );
+    }
+    else if (matrix.symmetric() || matrix.asymmetric())
+    {
+        return BlockLduSolver<Type>::New
+        (
+            solverName,
+            fieldName,
+            matrix,
+            dict
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "BlockLduSolver<Type>::New\n"
+            "(\n"
+            "    const word& fieldName,\n"
+            "    BlockLduMatrix<Type>& matrix,\n"
+            "    const dictionary& dict\n"
+            ")"
+        )   << "cannot solve incomplete matrix, "
+               "no diagonal or off-diagonal coefficient"
+            << exit(FatalError);
+
+        return autoPtr<BlockLduSolver<Type> >(NULL);
+    }
+}
+
+
+template<class Type>
+Foam::autoPtr<typename Foam::BlockLduSolver<Type> >
+Foam::BlockLduSolver<Type>::New
+(
+    const word& solverName,
+    const word& fieldName,
+    BlockLduMatrix<Type>& matrix,
+    const dictionary& dict
+)
+{
     if (matrix.diagonal())
     {
         return autoPtr<BlockLduSolver<Type> >
@@ -143,6 +189,7 @@ Foam::BlockLduSolver<Type>::New
             (
                 "BlockLduSolver<Type>::New\n"
                 "(\n"
+                "    const word& solverName,\n"
                 "    const word& fieldName,\n"
                 "    BlockLduMatrix<Type>& matrix,\n"
                 "    const dictionary& dict\n"
@@ -173,6 +220,7 @@ Foam::BlockLduSolver<Type>::New
             (
                 "BlockLduSolver<Type>::New\n"
                 "(\n"
+                "    const word& solverName,\n"
                 "    const word& fieldName,\n"
                 "    BlockLduMatrix<Type>& matrix,\n"
                 "    const dictionary& dict\n"
@@ -189,6 +237,7 @@ Foam::BlockLduSolver<Type>::New
             (
                 "BlockLduSolver<Type>::New\n"
                 "(\n"
+                "    const word& solverName,\n"
                 "    const word& fieldName,\n"
                 "    BlockLduMatrix<Type>& matrix,\n"
                 "    const dictionary& dict\n"
@@ -217,6 +266,7 @@ Foam::BlockLduSolver<Type>::New
         (
             "BlockLduSolver<Type>::New\n"
             "(\n"
+            "    const word& solverName,\n"
             "    const word& fieldName,\n"
             "    BlockLduMatrix<Type>& matrix,\n"
             "    const dictionary& dict\n"

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockLduSolver.H

@@ -172,7 +172,7 @@ public:
 
     // Selectors
 
-        //- Return a new solver
+        //- Return a new solver from a dictionary
         static autoPtr<BlockLduSolver<Type> > New
         (
             const word& fieldName,
@@ -180,6 +180,15 @@ public:
             const dictionary& dict
         );
 
+        //- Return a new solver given type
+        static autoPtr<BlockLduSolver<Type> > New
+        (
+            const word& solverName,
+            const word& fieldName,
+            BlockLduMatrix<Type>& matrix,
+            const dictionary& dict
+        );
+
 
     // Destructor
 

src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.C

@@ -74,21 +74,6 @@ Foam::BlockSolverPerformance<Type> Foam::SegregatedSolver<Type>::solve
         switchingDiag = true;
     }
 
-    // Check switching upper
-    bool switchingUpper = false;
-
-    if (blockMatrix.thereIsUpper())
-    {
-        if (blockMatrix.upper().activeType() == blockCoeffBase::SCALAR)
-        {
-            scalarMatrix_.upper() = blockMatrix.upper().asScalar();
-        }
-        else
-        {
-            switchingUpper = true;
-        }
-    }
-
     // Check switching lower
     bool switchingLower = false;
 
@@ -104,6 +89,21 @@ Foam::BlockSolverPerformance<Type> Foam::SegregatedSolver<Type>::solve
         }
     }
 
+    // Check switching upper
+    bool switchingUpper = false;
+
+    if (blockMatrix.thereIsUpper())
+    {
+        if (blockMatrix.upper().activeType() == blockCoeffBase::SCALAR)
+        {
+            scalarMatrix_.upper() = blockMatrix.upper().asScalar();
+        }
+        else
+        {
+            switchingUpper = true;
+        }
+    }
+
     // Decouple quadratic coupling by multiplying out the square coefficient
     // coupling
     Field<Type>* dBPtr = NULL;
@@ -120,7 +120,7 @@ Foam::BlockSolverPerformance<Type> Foam::SegregatedSolver<Type>::solve
     }
 
     // Prepare solver performance
-    word segSolverName(this->dict().lookup("solver"));
+    word segSolverName(this->dict().lookup("segSolver"));
 
     BlockSolverPerformance<Type> solverPerf
     (
@@ -156,20 +156,21 @@ Foam::BlockSolverPerformance<Type> Foam::SegregatedSolver<Type>::solve
             scalarMatrix_.diag() = blockMatrix.diag().component(cmpt);
         }
 
-        if (switchingUpper)
-        {
-            scalarMatrix_.upper() = blockMatrix.upper().component(cmpt);
-        }
-
         if (switchingLower)
         {
             scalarMatrix_.lower() = blockMatrix.lower().component(cmpt);
         }
 
+        if (switchingUpper)
+        {
+            scalarMatrix_.upper() = blockMatrix.upper().component(cmpt);
+        }
+
         // Call the scalar solver
         BlockSolverPerformance<scalar> scalarPerf =
             blockScalarSolver::New
             (
+                segSolverName,
                 this->fieldName(),
                 scalarMatrix_,
                 this->dict()

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/processorLduInterface.H

@@ -139,7 +139,8 @@ public:
                 UList<Type>&
             ) const;
 
-            //- Raw field receive function with data compression returning field
+            //- Raw field receive function with data compression,
+            //  returning a field
             template<class Type>
             tmp<Field<Type> > compressedReceive
             (

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C

@@ -46,9 +46,14 @@ void Foam::GAMGSolver::agglomerateMatrix(const label fineLevelIndex)
     const labelList& faceRestrictAddr =
         agglomeration_.faceRestrictAddressing(fineLevelIndex);
 
-    // Coarse matrix diagonal initialised by restricting the finer mesh diagonal
+    // Coarse matrix diagonal initialised by restricting the fine mesh diagonal
     scalarField& coarseDiag = coarseMatrix.diag();
-    agglomeration_.restrictField(coarseDiag, fineMatrix.diag(), fineLevelIndex);
+    agglomeration_.restrictField
+    (
+        coarseDiag,
+        fineMatrix.diag(),
+        fineLevelIndex
+    );
 
     // Get reference to fine-level interfaces
     const lduInterfaceFieldPtrsList& fineInterfaces =

src/VectorN/OpenFOAM/BlockLduMatrices/blockVectorNSolvers.C

@@ -43,6 +43,7 @@ License
 #include "BlockCGSolver.H"
 #include "BlockGaussSeidelSolver.H"
 #include "BlockGMRESSolver.H"
+#include "SegregatedSolver.H"
 
 #include "VectorTensorNFields.H"
 #include "ExpandTensorN.H"
@@ -105,23 +106,29 @@ defineTemplateRunTimeSelectionTable                                           \
 typedef BlockDiagonalSolver<type > block##Type##DiagonalSolver;               \
 defineNamedTemplateTypeNameAndDebug(block##Type##DiagonalSolver, 0);          \
                                                                               \
-typedef BlockBiCGStabSolver<type > block##Type##BiCGStabSolver;               \
+typedef BlockGaussSeidelSolver<type > block##Type##GaussSeidelSolver;         \
-makeBlockSolverTypeName(block##Type##BiCGStabSolver);                         \
+makeBlockSolverTypeName(block##Type##GaussSeidelSolver);                      \
-addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, symMatrix);         \
+addSolverToBlockMatrix(Type, block##Type##GaussSeidelSolver, symMatrix);      \
-addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, asymMatrix);        \
+addSolverToBlockMatrix(Type, block##Type##GaussSeidelSolver, asymMatrix);     \
                                                                               \
 typedef BlockCGSolver<type > block##Type##CGSolver;                           \
 makeBlockSolverTypeName(block##Type##CGSolver);                               \
 addSolverToBlockMatrix(Type, block##Type##CGSolver, symMatrix);               \
                                                                               \
-typedef BlockGaussSeidelSolver<type > block##Type##GaussSeidelSolver;         \
+typedef BlockBiCGStabSolver<type > block##Type##BiCGStabSolver;               \
-makeBlockSolverTypeName(block##Type##GaussSeidelSolver);                      \
+makeBlockSolverTypeName(block##Type##BiCGStabSolver);                         \
-addSolverToBlockMatrix(Type, block##Type##GaussSeidelSolver, symMatrix);      \
+addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, symMatrix);         \
+addSolverToBlockMatrix(Type, block##Type##BiCGStabSolver, asymMatrix);        \
                                                                               \
 typedef BlockGMRESSolver<type > block##Type##GMRESSolver;                     \
 makeBlockSolverTypeName(block##Type##GMRESSolver);                            \
 addSolverToBlockMatrix(Type, block##Type##GMRESSolver, symMatrix);            \
-addSolverToBlockMatrix(Type, block##Type##GMRESSolver, asymMatrix);
+addSolverToBlockMatrix(Type, block##Type##GMRESSolver, asymMatrix);           \
+                                                                              \
+typedef SegregatedSolver<type > block##Type##SegregatedSolver;                \
+makeBlockSolverTypeName(block##Type##SegregatedSolver);                       \
+addSolverToBlockMatrix(Type, block##Type##SegregatedSolver, symMatrix);       \
+addSolverToBlockMatrix(Type, block##Type##SegregatedSolver, asymMatrix);
 
 
 forAllVectorNTypes(makeSolver)

src/VectorN/OpenFOAM/Fields/VectorNFieldTypes.H

@@ -32,46 +32,54 @@ Description
 #define VectorNFieldTypes_H
 
 #include "vector2.H"
+#include "vector3.H"
 #include "vector4.H"
 #include "vector6.H"
 #include "vector8.H"
 
 #include "tensor2.H"
+#include "tensor3.H"
 #include "tensor4.H"
 #include "tensor6.H"
 #include "tensor8.H"
 
 #include "diagTensor2.H"
+#include "diagTensor3.H"
 #include "diagTensor4.H"
 #include "diagTensor6.H"
 #include "diagTensor8.H"
 
 #define forAllVectorNTypes(m, args...)              \
     m(vector2, Vector2, args)                       \
+    m(vector3, Vector3, args)                       \
     m(vector4, Vector4, args)                       \
     m(vector6, Vector6, args)                       \
     m(vector8, Vector8, args)
 
 #define forAllTensorNTypes(m, args...)              \
     m(tensor2, Tensor2, args)                       \
+    m(tensor3, Tensor3, args)                       \
     m(tensor4, Tensor4, args)                       \
     m(tensor6, Tensor6, args)                       \
     m(tensor8, Tensor8, args)
 
 #define forAllDiagTensorNTypes(m, args...)          \
     m(diagTensor2, DiagTensor2, args)               \
+    m(diagTensor3, DiagTensor3, args)               \
     m(diagTensor4, DiagTensor4, args)               \
     m(diagTensor6, DiagTensor6, args)               \
     m(diagTensor8, DiagTensor8, args)
 
 #define forAllSphericalTensorNTypes(m, args...)     \
     m(sphericalTensor2, SphericalTensor2, args)     \
+    m(sphericalTensor3, SphericalTensor3, args)     \
     m(sphericalTensor4, SphericalTensor4, args)     \
     m(sphericalTensor6, SphericalTensor6, args)     \
     m(sphericalTensor8, SphericalTensor8, args)
 
 #define forAllVectorTensorNTypes(m, args...)                            \
     m(tensor2, diagTensor2, sphericalTensor2, vector2, scalar, args)    \
+    m(tensor3, diagTensor3, sphericalTensor3, vector3, scalar, args)    \
     m(tensor4, diagTensor4, sphericalTensor4, vector4, scalar, args)    \
     m(tensor6, diagTensor6, sphericalTensor6, vector6, scalar, args)    \
     m(tensor8, diagTensor8, sphericalTensor8, vector8, scalar, args)

src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.C

@@ -51,6 +51,21 @@ void blockInsert
 }
 
 
+template<class BlockType>
+void blockAdd
+(
+    const direction dir,
+    const scalarField& x,
+    Field<BlockType>& blockX
+)
+{
+    forAll (x, i)
+    {
+        blockX[i](dir) += x[i];
+    }
+}
+
+
 template<class BlockType>
 void blockRetrieve
 (
@@ -131,6 +146,46 @@ void insertUpperLower
         return;
     }
 
+    if (m.symmetric() && blockM.symmetric())
+    {
+        Info<< "Both m and blockM are symmetric: inserting only upper triangle"
+            << endl;
+    }
+    else
+    {
+        // Either scalar or block matrix is asymmetric: insert lower triangle
+        const scalarField& lower = m.lower();
+
+        if (blockM.lower().activeType() == blockCoeffBase::UNALLOCATED)
+        {
+            blockM.lower().asScalar() = lower;
+        }
+        else if
+        (
+            blockM.lower().activeType() == blockCoeffBase::SCALAR
+         || blockM.lower().activeType() == blockCoeffBase::LINEAR
+        )
+        {
+            typename CoeffField<BlockType>::linearTypeField& blockLower =
+                blockM.lower().asLinear();
+
+            forAll (lower, i)
+            {
+                blockLower[i](dir) = lower[i];
+            }
+        }
+        else if (blockM.lower().activeType() == blockCoeffBase::SQUARE)
+        {
+            typename CoeffField<BlockType>::squareTypeField& blockLower =
+                blockM.lower().asSquare();
+
+            forAll (lower, i)
+            {
+                blockLower[i](dir, dir) = lower[i];
+            }
+        }
+    }
+
     if (m.hasUpper())
     {
         const scalarField& upper = m.upper();
@@ -178,42 +233,61 @@ void insertUpperLower
             << abort(FatalError);
     }
 
-    if (m.symmetric() && blockM.symmetric())
+    // Insert block interface fields
+    forAll(blockM.interfaces(), patchI)
     {
-        Info<< "Both m and blockM are symmetric: inserting only upper triangle"
+        if (blockM.interfaces().set(patchI))
-            << endl;
-    }
-    else
         {
-        // Either scalar or block matrix is asymmetric: insert lower triangle
+            // Couple upper and lower
-        const scalarField& lower = m.lower();
+            const scalarField& cUpper = m.boundaryCoeffs()[patchI];
+            const scalarField& cLower = m.internalCoeffs()[patchI];
 
-        if (blockM.lower().activeType() == blockCoeffBase::UNALLOCATED)
+            if
+            (
+                blockM.coupleUpper()[patchI].activeType()
+             == blockCoeffBase::UNALLOCATED
+            )
             {
-            blockM.lower().asScalar() = lower;
+                blockM.coupleUpper()[patchI].asScalar() = cUpper;
+                blockM.coupleLower()[patchI].asScalar() = cLower;
             }
             else if
             (
-            blockM.lower().activeType() == blockCoeffBase::SCALAR
+                blockM.coupleUpper()[patchI].activeType()
-         || blockM.lower().activeType() == blockCoeffBase::LINEAR
+             == blockCoeffBase::SCALAR
+             || blockM.coupleUpper()[patchI].activeType()
+             == blockCoeffBase::LINEAR
             )
             {
+                typename CoeffField<BlockType>::linearTypeField& blockUpper =
+                    blockM.coupleUpper()[patchI].asLinear();
+
                 typename CoeffField<BlockType>::linearTypeField& blockLower =
-                blockM.lower().asLinear();
+                    blockM.coupleLower()[patchI].asLinear();
 
-            forAll (lower, i)
+                forAll (cUpper, i)
                 {
-                blockLower[i](dir) = lower[i];
+                    blockUpper[i](dir) = cUpper[i];
+                    blockLower[i](dir) = cLower[i];
                 }
             }
-        else if (blockM.lower().activeType() == blockCoeffBase::SQUARE)
+            else if
+            (
+                blockM.coupleUpper()[patchI].activeType()
+             == blockCoeffBase::SQUARE
+            )
             {
+                typename CoeffField<BlockType>::squareTypeField& blockUpper =
+                    blockM.coupleUpper()[patchI].asSquare();
+
                 typename CoeffField<BlockType>::squareTypeField& blockLower =
-                blockM.lower().asSquare();
+                    blockM.coupleLower()[patchI].asSquare();
 
-            forAll (lower, i)
+                forAll (cUpper, i)
                 {
-                blockLower[i](dir, dir) = lower[i];
+                    blockUpper[i](dir, dir) = cUpper[i];
+                    blockLower[i](dir, dir) = cLower[i];
+                }
             }
         }
     }
@@ -236,6 +310,142 @@ void insertEquation
 }
 
 
+template<class BlockType>
+void insertCouplingDiagSource
+(
+    const direction dirI,
+    const direction dirJ,
+    const fvScalarMatrix& m,
+    BlockLduMatrix<BlockType>& blockM,
+    Field<BlockType>& blockB
+)
+{
+    // Prepare the diagonal and source
+
+    scalarField diag = m.diag();
+    scalarField source = m.source();
+
+    // Add boundary source contribution
+    m.addBoundaryDiag(diag, 0);
+    m.addBoundarySource(source, false);
+
+
+    // Add off-diagonal block coefficients
+    typename CoeffField<BlockType>::squareTypeField& blockDiag =
+        blockM.diag().asSquare();
+
+    // Set off-diagonal coefficient
+    forAll (diag, i)
+    {
+        blockDiag[i](dirI, dirJ) = diag[i];
+    }
+
+    blockAdd(dirI, source, blockB);
+}
+
+
+template<class BlockType>
+void insertCouplingUpperLower
+(
+    const direction dirI,
+    const direction dirJ,
+    const fvScalarMatrix& m,
+    BlockLduMatrix<BlockType>& blockM
+)
+{
+    if (m.diagonal())
+    {
+        // Matrix for insertion is diagonal-only: nothing to do
+        return;
+    }
+
+    if (m.symmetric() && blockM.symmetric())
+    {
+        Info<< "Both m and blockM are symmetric: inserting only upper triangle"
+            << endl;
+    }
+    else
+    {
+        // Either scalar or block matrix is asymmetric: insert lower triangle
+        const scalarField& lower = m.lower();
+
+        typename CoeffField<BlockType>::squareTypeField& blockLower =
+            blockM.lower().asSquare();
+
+        forAll (lower, i)
+        {
+            blockLower[i](dirI, dirJ) = lower[i];
+        }
+    }
+
+    if (m.hasUpper())
+    {
+        const scalarField& upper = m.upper();
+
+        typename CoeffField<BlockType>::squareTypeField& blockUpper =
+            blockM.upper().asSquare();
+
+        forAll (upper, i)
+        {
+            blockUpper[i](dirI, dirJ) = upper[i];
+        }
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void insertCouplingUpperLower\n"
+            "(\n"
+            "    const direction dirI,\n"
+            "    const direction dirJ,\n"
+            "    const fvScalarMatrix& m,\n"
+            "    BlockLduMatrix<BlockType>& blockM\n"
+            ")"
+        )   << "Error in matrix insertion: problem with block structure"
+            << abort(FatalError);
+    }
+
+    // Insert block interface fields
+    forAll(blockM.interfaces(), patchI)
+    {
+        if (blockM.interfaces().set(patchI))
+        {
+            // Couple upper and lower
+            const scalarField& cUpper = m.boundaryCoeffs()[patchI];
+            const scalarField& cLower = m.internalCoeffs()[patchI];
+
+            typename CoeffField<BlockType>::squareTypeField& blockUpper =
+                blockM.coupleUpper()[patchI].asSquare();
+
+            typename CoeffField<BlockType>::squareTypeField& blockLower =
+                blockM.coupleLower()[patchI].asSquare();
+
+            forAll (cUpper, i)
+            {
+                blockUpper[i](dirI, dirJ) = cUpper[i];
+                blockLower[i](dirI, dirJ) = cLower[i];
+            }
+        }
+    }
+}
+
+
+template<class BlockType>
+void insertCoupling
+(
+    const direction dirI,
+    const direction dirJ,
+    const fvScalarMatrix& m,
+    BlockLduMatrix<BlockType>& blockM,
+    Field<BlockType>& blockX,
+    Field<BlockType>& blockB
+)
+{
+    insertCouplingDiagSource(dirI, dirJ, m, blockM, blockB);
+    insertCouplingUpperLower(dirI, dirJ, m, blockM);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace blockMatrixTools

src/VectorN/finiteVolume/blockMatrixTools/blockMatrixTools.H

@@ -0,0 +1,179 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+
+InNamespace
+    Foam::
+
+Description
+    Block matrix insertion and retrieval tools
+
+SourceFiles
+    blockMatrixTools.C
+
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef blockMatrixTools_H
+#define blockMatrixTools_H
+
+#include "blockLduMatrices.H"
+#include "fvMatrices.H"
+#include "surfaceFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Namespace blockMatrixTools functions Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace blockMatrixTools
+{
+    // Field operations
+
+        //- Insert field into block field
+        template<class BlockType>
+        void blockInsert
+        (
+            const direction dir,
+            const scalarField& x,
+            Field<BlockType>& blockX
+        );
+
+        //- Add field into block field
+        template<class BlockType>
+        void blockAdd
+        (
+            const direction dir,
+            const scalarField& x,
+            Field<BlockType>& blockX
+        );
+
+        //- Retrieve field from block field
+        template<class BlockType>
+        void blockRetrieve
+        (
+            const direction dir,
+            scalarField& x,
+            const Field<BlockType>& blockX
+        );
+
+
+    // Diagonal block operations
+
+        //- Insert matrix diagonal and source into the block system
+        template<class BlockType>
+        void insertDiagSource
+        (
+            const direction dir,
+            const fvScalarMatrix& m,
+            BlockLduMatrix<BlockType>& blockM,
+            Field<BlockType>& blockB
+        );
+
+        // Insert upper and lower coefficients matrix into the block system
+        template<class Type, class BlockType>
+        void insertUpperLower
+        (
+            const direction dir,
+            const fvScalarMatrix& m,
+            BlockLduMatrix<BlockType>& blockM
+        );
+
+        // Insert matrix into the block system
+        template<class BlockType>
+        void insertEquation
+        (
+            const direction dir,
+            const fvScalarMatrix& m,
+            BlockLduMatrix<BlockType>& blockM,
+            Field<BlockType>& blockX,
+            Field<BlockType>& blockB
+        );
+
+
+    // Coupling block operations
+
+        //- Insert coupling matrix diagonal and source into the block system
+        template<class BlockType>
+        void insertCouplingDiagSource
+        (
+            const direction dirI,
+            const direction dirJ,
+            const fvScalarMatrix& m,
+            BlockLduMatrix<BlockType>& blockM,
+            Field<BlockType>& blockB
+        );
+
+        // Insert coupling matrix into the block system
+        template<class Type, class BlockType>
+        void insertCouplingUpperLower
+        (
+            const direction dirI,
+            const direction dirJ,
+            const fvScalarMatrix& m,
+            BlockLduMatrix<BlockType>& blockM
+        );
+
+        // Insert coupling matrix into the block system
+        template<class BlockType>
+        void insertCoupling
+        (
+            const direction dir,
+            const fvScalarMatrix& m,
+            BlockLduMatrix<BlockType>& blockM,
+            Field<BlockType>& blockX,
+            Field<BlockType>& blockB
+        );
+
+        // Update coupling of block system
+        // Subtracts the block-coefficient coupling as specified by the user
+        // from the source, leaving the implicit update given by linearisation
+        template<class BlockType>
+        void updateCoupling
+        (
+            BlockLduMatrix<BlockType>& blockM,
+            Field<BlockType>& x,
+            Field<BlockType>& b
+        );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "blockMatrixTools.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/VectorN/finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFieldsFwd.H

@@ -39,21 +39,25 @@ namespace Foam
 template<class Type> class fvsPatchField;
 
 typedef fvsPatchField<vector2> fvsPatchVector2Field;
+typedef fvsPatchField<vector3> fvsPatchVector3Field;
 typedef fvsPatchField<vector4> fvsPatchVector4Field;
 typedef fvsPatchField<vector6> fvsPatchVector6Field;
 typedef fvsPatchField<vector8> fvsPatchVector8Field;
 
 typedef fvsPatchField<tensor2> fvsPatchTensor2Field;
+typedef fvsPatchField<tensor3> fvsPatchTensor3Field;
 typedef fvsPatchField<tensor4> fvsPatchTensor4Field;
 typedef fvsPatchField<tensor6> fvsPatchTensor6Field;
 typedef fvsPatchField<tensor8> fvsPatchTensor8Field;
 
 typedef fvsPatchField<diagTensor2> fvsPatchDiagTensor2Field;
+typedef fvsPatchField<diagTensor3> fvsPatchDiagTensor3Field;
 typedef fvsPatchField<diagTensor4> fvsPatchDiagTensor4Field;
 typedef fvsPatchField<diagTensor6> fvsPatchDiagTensor6Field;
 typedef fvsPatchField<diagTensor8> fvsPatchDiagTensor8Field;
 
 typedef fvsPatchField<sphericalTensor2> fvsPatchSphericalTensor2Field;
+typedef fvsPatchField<sphericalTensor3> fvsPatchSphericalTensor3Field;
 typedef fvsPatchField<sphericalTensor4> fvsPatchSphericalTensor4Field;
 typedef fvsPatchField<sphericalTensor6> fvsPatchSphericalTensor6Field;
 typedef fvsPatchField<sphericalTensor8> fvsPatchSphericalTensor8Field;

src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C

@@ -90,8 +90,9 @@ void Foam::setRefCell
                      "    bool\n"
                      ")",
                     dict
-                )   << "Unable to set reference cell for field " << field.name()
+                )   << "Unable to set reference cell for field "
-                    << nl << "    Reference point " << refPointName
+                    << field.name() << nl
+                    << "    Reference point " << refPointName
                     << " " << refPointi
                     << " found on " << sumHasRef << " domains (should be one)"
                     << nl << exit(FatalIOError);

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -276,6 +276,13 @@ public:
                 return source_;
             }
 
+            //- Access to fvBoundary scalar field containing
+            //  pseudo-matrix coeffs for internal cells
+            const FieldField<Field, Type>& internalCoeffs() const
+            {
+                return internalCoeffs_;
+            }
+
             //- fvBoundary scalar field containing pseudo-matrix coeffs
             //  for internal cells
             FieldField<Field, Type>& internalCoeffs()
@@ -283,6 +290,13 @@ public:
                 return internalCoeffs_;
             }
 
+            //- Access to fvBoundary scalar field containing
+            //  pseudo-matrix coeffs for boundary cells
+            const FieldField<Field, Type>& boundaryCoeffs() const
+            {
+                return boundaryCoeffs_;
+            }
+
             //- fvBoundary scalar field containing pseudo-matrix coeffs
             //  for boundary cells
             FieldField<Field, Type>& boundaryCoeffs()