FEATURE: Block solve reorganisation

applications/solvers/coupled/Allwmake

@@ -8,4 +8,6 @@ wmake blockCoupledScalarTransportFoam
 wmake conjugateHeatFoam
 wmake conjugateHeatSimpleFoam
 wmake pUCoupledFoam
+wmake pUCoupledFullPicardFoam
+wmake pUCoupledSemiPicardFoam
 # ----------------------------------------------------------------- end-of-file

applications/solvers/coupled/blockCoupledScalarTransportFoam/Make/options

@@ -4,5 +4,4 @@ EXE_INC = \
     -I$(LIB_SRC)/VectorN/lnInclude
 
 EXE_LIBS = \
-    -lfiniteVolume \
-    -lVectorN
+    -lfiniteVolume

applications/solvers/coupled/conjugateHeatTransfer/Make/options

@@ -8,6 +8,4 @@ LIB_LIBS = \
     -lfiniteVolume \
     -lspecie \
     -lbasicThermophysicalModels \
-    -lradiation \
-    -lVectorN
-
+    -lradiation

applications/solvers/coupled/pUCoupledFoam/Make/options

@@ -3,13 +3,10 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/VectorN/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -llduSolvers \
-    -lVectorN
-
+    -llduSolvers

applications/solvers/coupled/pUCoupledFoam/UEqn.H

@@ -7,4 +7,4 @@ fvVectorMatrix UEqn
 
 UEqn.relax();
 
-blockMatrixTools::insertEquation(0, UEqn, A, x, b);
+UpEqn.insertEquation(0, UEqn);

applications/solvers/coupled/pUCoupledFoam/calculateCouplingMatrices.H

@@ -1,2 +0,0 @@
-blockVectorMatrix pInU(fvm::grad(p));
-blockVectorMatrix UInp(fvm::div(U));

applications/solvers/coupled/pUCoupledFoam/initializeBlockMatrix.H

@@ -1,9 +0,0 @@
-// Create block matrix
-BlockLduMatrix<vector4> A(mesh);
-
-// Block matrix - create x and b
-vector4Field& x = Up.internalField();
-vector4Field b(x.size(), vector4::zero);
-
-// Set block interfaces properly
-A.interfaces() = Up.boundaryField().blockInterfaces();

applications/solvers/coupled/pUCoupledFoam/pEqn.H

@@ -8,7 +8,7 @@ surfaceScalarField rUAf
 surfaceScalarField presSource
 (
     "presSource",
-    rUAf*(fvc::interpolate(fvc::grad(p)) & mesh.Sf())
+    rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
 );
 
 fvScalarMatrix pEqn
@@ -20,4 +20,4 @@ fvScalarMatrix pEqn
 
 pEqn.setReference(pRefCell, pRefValue);
 
-blockMatrixTools::insertEquation(3, pEqn, A, x, b);
+UpEqn.insertEquation(3, pEqn);

applications/solvers/coupled/pUCoupledFoam/pUCoupledFoam.C

@@ -1,32 +1,33 @@
 /*---------------------------------------------------------------------------*\
   =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | For copyright notice see file Copyright
+    \\  /    A nd           | Copyright held by original author
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
-    This file is part of foam-extend.
+    This file is part of OpenFOAM.
 
-    foam-extend is free software: you can redistribute it and/or modify it
+    OpenFOAM is free software; you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
+    Free Software Foundation; either version 2 of the License, or (at your
     option) any later version.
 
-    foam-extend 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.
+    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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
     pUCoupledFoam
 
 Description
     Steady-state solver for incompressible, turbulent flow, with implicit
-    coupling between pressure and velocity achieved by BlockLduMatrix
+    coupling between pressure and velocity achieved by fvBlockMatrix.
     Turbulence is in this version solved using the existing turbulence
     structure.
 
@@ -39,13 +40,7 @@ Authors
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
-
-#include "VectorNFieldTypes.H"
-#include "volVectorNFields.H"
-#include "blockLduSolvers.H"
-#include "blockVectorNMatrices.H"
-
-#include "blockMatrixTools.H"
+#include "fvBlockMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,12 +54,8 @@ int main(int argc, char *argv[])
 #   include "initContinuityErrs.H"
 #   include "readBlockSolverControls.H"
 
-    // Calculate coupling matrices only once since the mesh doesn't change and
-    // implicit div and grad operators are only dependant on Sf. Actually
-    // coupling terms (div(U) and grad(p)) in blockMatrix do not change,
-    // so they could be inserted only once, resetting other parts of
-    // blockMatrix to zero at the end of each time step. VV, 30/April/2014
-#   include "calculateCouplingMatrices.H"
+    // Calculate div U coupling only once since it is only geometry dependant.
+    BlockLduSystem<vector, scalar> UInp(fvm::div(U));
 
     Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
@@ -73,8 +64,8 @@ int main(int argc, char *argv[])
 
         p.storePrevIter();
 
-        // Initialize block matrix
-#       include "initializeBlockMatrix.H"
+        // Initialize the Up block system (matrix, source and reference to Up)
+        fvBlockMatrix<vector4> UpEqn(Up);
 
         // Assemble and insert momentum equation
 #       include "UEqn.H"
@@ -82,27 +73,22 @@ int main(int argc, char *argv[])
         // Assemble and insert pressure equation
 #       include "pEqn.H"
 
-        // Insert coupling, updating the boundary contributions
-        // Last argument in insertBlockCoupling says if the first location
+        // Calculate grad p coupling matrix. Needs to be here if one uses
+        // gradient schemes with limiters. VV, 9/June/2014
+        BlockLduSystem<vector, vector> pInU(fvm::grad(p));
+
+        // Last argument in insertBlockCoupling says if the column direction
         // should be incremented. This is needed for arbitrary positioning
         // of U and p in the system. This could be better. VV, 30/April/2014
-        blockMatrixTools::insertBlockCoupling(3, 0, UInp, U, A, b, false);
-        blockMatrixTools::insertBlockCoupling(0, 3, pInU, p, A, b, true);
+        UpEqn.insertBlockCoupling(0, 3, pInU, true);
+        UpEqn.insertBlockCoupling(3, 0, UInp, false);
 
         // Solve the block matrix
-        BlockSolverPerformance<vector4> solverPerf =
-            BlockLduSolver<vector4>::New
-            (
-                word("Up"),
-                A,
-                mesh.solutionDict().solver("Up")
-            )->solve(Up, b);
-
-        solverPerf.print();
+        UpEqn.solve();
 
         // Retrieve solution
-        blockMatrixTools::retrieveSolution(0, U.internalField(), Up);
-        blockMatrixTools::retrieveSolution(3, p.internalField(), Up);
+        UpEqn.retrieveSolution(0, U.internalField());
+        UpEqn.retrieveSolution(3, p.internalField());
 
         U.correctBoundaryConditions();
         p.correctBoundaryConditions();

applications/solvers/coupled/pUCoupledFullPicardFoam/Make/files

@@ -0,0 +1,3 @@
+pUCoupledFullPicardFoam.C
+
+EXE = $(FOAM_APPBIN)/pUCoupledFullPicardFoam

applications/solvers/coupled/pUCoupledFullPicardFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -llduSolvers

applications/solvers/coupled/pUCoupledFullPicardFoam/UEqn.H

@@ -0,0 +1,12 @@
+// Momentum equation
+fvVectorMatrix UEqn
+(
+    fvm::div(phi, U)
+  + fvc::div(-phi, U, "div(phi,U)")
+  + turbulence->divDevReff(U)
+);
+
+UEqn.relax();
+
+UpEqn.insertEquation(0, UEqn);
+UpEqn.insertPicardTensor(0, U, phi);

applications/solvers/coupled/pUCoupledFullPicardFoam/createFields.H

@@ -0,0 +1,52 @@
+Info << "Reading field p\n" << endl;
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info << "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "createPhi.H"
+
+singlePhaseTransportModel laminarTransport(U, phi);
+autoPtr<incompressible::RASModel> turbulence
+(
+    incompressible::RASModel::New(U, phi, laminarTransport)
+);
+
+// Block vector field for velocity (first entry) and pressure (second
+// entry).
+Info << "Creating field Up\n" << endl;
+volVector4Field Up
+(
+    IOobject
+    (
+       "Up",
+       runTime.timeName(),
+       mesh,
+       IOobject::NO_READ,
+       IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedVector4("zero", dimless, vector4::zero)
+);

applications/solvers/coupled/pUCoupledFullPicardFoam/pEqn.H

@@ -0,0 +1,23 @@
+// Pressure parts of the continuity equation
+surfaceScalarField rUAf
+(
+    "rUAf",
+    fvc::interpolate(1.0/UEqn.A())
+);
+
+surfaceScalarField presSource
+(
+    "presSource",
+    rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
+);
+
+fvScalarMatrix pEqn
+(
+  - fvm::laplacian(rUAf, p)
+ ==
+  - fvc::div(presSource)
+);
+
+pEqn.setReference(pRefCell, pRefValue);
+
+UpEqn.insertEquation(3, pEqn);

applications/solvers/coupled/pUCoupledFullPicardFoam/pUCoupledFullPicardFoam.C

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Application
+    pUCoupledFullPicardFoam
+
+Description
+    Steady-state solver for incompressible, turbulent flow, with implicit
+    coupling between pressure and velocity achieved by fvBlockMatrix.
+    Turbulence is in this version solved using the existing turbulence
+    structure.
+    Convective term is properly linearised using Picard linearisation and
+    additional block coefficients. This is full, proper variant!
+    VV, 23/July/2014.
+
+Authors
+    Klas Jareteg, Chalmers University of Technology,
+    Vuko Vukcevic, FMENA Zagreb.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+#include "fvBlockMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+#   include "readBlockSolverControls.H"
+
+    // Calculate div U coupling only once since it is only geometry dependant.
+    BlockLduSystem<vector, scalar> UInp(fvm::div(U));
+
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        p.storePrevIter();
+
+        // Initialize the Up block system (matrix, source and reference to Up)
+        fvBlockMatrix<vector4> UpEqn(Up);
+
+        // Assemble and insert momentum equation
+#       include "UEqn.H"
+
+        // Assemble and insert pressure equation
+#       include "pEqn.H"
+
+        // Calculate grad p coupling matrix. Needs to be here if one uses
+        // gradient schemes with limiters. VV, 9/June/2014
+        BlockLduSystem<vector, vector> pInU(fvm::grad(p));
+
+        // Last argument in insertBlockCoupling says if the column direction
+        // should be incremented. This is needed for arbitrary positioning
+        // of U and p in the system. This could be better. VV, 30/April/2014
+        UpEqn.insertBlockCoupling(0, 3, pInU, true);
+        UpEqn.insertBlockCoupling(3, 0, UInp, false);
+
+        // Solve the block matrix
+        UpEqn.solve();
+
+        // Retrieve solution
+        UpEqn.retrieveSolution(0, U.internalField());
+        UpEqn.retrieveSolution(3, p.internalField());
+
+        U.correctBoundaryConditions();
+        p.correctBoundaryConditions();
+
+        phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;
+
+#       include "continuityErrs.H"
+
+        p.relax();
+
+        turbulence->correct();
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}

applications/solvers/coupled/pUCoupledFullPicardFoam/readBlockSolverControls.H

@@ -0,0 +1,3 @@
+label pRefCell = 0;
+scalar pRefValue = 0;
+setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue);

applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/files

@@ -0,0 +1,3 @@
+pUCoupledSemiPicardFoam.C
+
+EXE = $(FOAM_APPBIN)/pUCoupledSemiPicardFoam

applications/solvers/coupled/pUCoupledSemiPicardFoam/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -llduSolvers

applications/solvers/coupled/pUCoupledSemiPicardFoam/UEqn.H

@@ -0,0 +1,11 @@
+// Momentum equation
+fvVectorMatrix UEqn
+(
+    fvm::div(2*phi, U,"div(phi,U)")
+  + fvc::div(-phi, U,"div(phi,U)")
+  + turbulence->divDevReff(U)
+);
+
+UEqn.relax();
+
+UpEqn.insertEquation(0, UEqn);

applications/solvers/coupled/pUCoupledSemiPicardFoam/createFields.H

@@ -0,0 +1,52 @@
+Info << "Reading field p\n" << endl;
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info << "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "createPhi.H"
+
+singlePhaseTransportModel laminarTransport(U, phi);
+autoPtr<incompressible::RASModel> turbulence
+(
+    incompressible::RASModel::New(U, phi, laminarTransport)
+);
+
+// Block vector field for velocity (first entry) and pressure (second
+// entry).
+Info << "Creating field Up\n" << endl;
+volVector4Field Up
+(
+    IOobject
+    (
+       "Up",
+       runTime.timeName(),
+       mesh,
+       IOobject::NO_READ,
+       IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedVector4("zero", dimless, vector4::zero)
+);

applications/solvers/coupled/pUCoupledSemiPicardFoam/pEqn.H

@@ -0,0 +1,23 @@
+// Pressure parts of the continuity equation
+surfaceScalarField rUAf
+(
+    "rUAf",
+    fvc::interpolate(1.0/UEqn.A())
+);
+
+surfaceScalarField presSource
+(
+    "presSource",
+    rUAf*(fvc::interpolate(fvc::grad(p, "grad(pSource)")) & mesh.Sf())
+);
+
+fvScalarMatrix pEqn
+(
+  - fvm::laplacian(rUAf, p)
+ ==
+  - fvc::div(presSource)
+);
+
+pEqn.setReference(pRefCell, pRefValue);
+
+UpEqn.insertEquation(3, pEqn);

applications/solvers/coupled/pUCoupledSemiPicardFoam/pUCoupledSemiPicardFoam.C

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Application
+    pUCoupledSemiPicardFoam
+
+Description
+    Steady-state solver for incompressible, turbulent flow, with implicit
+    coupling between pressure and velocity achieved by fvBlockMatrix.
+    Turbulence is in this version solved using the existing turbulence
+    structure.
+    Convective term is properly linearised using Picard linearisation and
+    additional block coefficients. This is semi variant where the difference
+    between succesive velocity fields is neglected! VV, 23/July/2014.
+
+Authors
+    Klas Jareteg, Chalmers University of Technology,
+    Vuko Vukcevic, FMENA Zagreb.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+#include "fvBlockMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+#   include "readBlockSolverControls.H"
+
+    // Calculate div U coupling only once since it is only geometry dependant.
+    BlockLduSystem<vector, scalar> UInp(fvm::div(U));
+
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        p.storePrevIter();
+
+        // Initialize the Up block system (matrix, source and reference to Up)
+        fvBlockMatrix<vector4> UpEqn(Up);
+
+        // Assemble and insert momentum equation
+#       include "UEqn.H"
+
+        // Assemble and insert pressure equation
+#       include "pEqn.H"
+
+        // Calculate grad p coupling matrix. Needs to be here if one uses
+        // gradient schemes with limiters. VV, 9/June/2014
+        BlockLduSystem<vector, vector> pInU(fvm::grad(p));
+
+        // Last argument in insertBlockCoupling says if the column direction
+        // should be incremented. This is needed for arbitrary positioning
+        // of U and p in the system. This could be better. VV, 30/April/2014
+        UpEqn.insertBlockCoupling(0, 3, pInU, true);
+        UpEqn.insertBlockCoupling(3, 0, UInp, false);
+
+        // Solve the block matrix
+        UpEqn.solve();
+
+        // Retrieve solution
+        UpEqn.retrieveSolution(0, U.internalField());
+        UpEqn.retrieveSolution(3, p.internalField());
+
+        U.correctBoundaryConditions();
+        p.correctBoundaryConditions();
+
+        phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;
+
+#       include "continuityErrs.H"
+
+        p.relax();
+
+        turbulence->correct();
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}

applications/solvers/coupled/pUCoupledSemiPicardFoam/readBlockSolverControls.H

@@ -0,0 +1,3 @@
+label pRefCell = 0;
+scalar pRefValue = 0;
+setRefCell(p, mesh.solutionDict().subDict("blockSolver"), pRefCell, pRefValue);

src/Allwmake

@@ -37,7 +37,6 @@ wmake libso finiteVolume
 wmake libso finiteArea
 
 wmake libso lduSolvers
-wmake libso VectorN
 
 wmake libso tetFiniteElement
 

src/VectorN/Make/files

@@ -1,28 +0,0 @@
-foam/DimensionedTypes/dimensionedVectorTensorN/dimensionedVectorTensorN.C
-
-finiteVolume/fields/fvPatchFields/fvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/genericFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/calculatedFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/emptyFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/transformFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/wedgeFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/coupledFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/processorFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/cyclicFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/fixedValueFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/zeroGradientFvPatchVectorNFields.C
-finiteVolume/fields/fvPatchFields/fixedGradientFvPatchVectorNFields.C
-
-finiteVolume/fields/fvsPatchFields/fvsPatchVectorNFields.C
-finiteVolume/fields/fvsPatchFields/calculatedFvsPatchVectorNFields.C
-finiteVolume/fields/fvsPatchFields/emptyFvsPatchVectorNFields.C
-finiteVolume/fields/fvsPatchFields/wedgeFvsPatchVectorNFields.C
-finiteVolume/fields/fvsPatchFields/coupledFvsPatchVectorNFields.C
-finiteVolume/fields/fvsPatchFields/processorFvsPatchVectorNFields.C
-
-finiteVolume/fields/volFields/volVectorNFields.C
-finiteVolume/fields/surfaceFields/surfaceVectorNFields.C
-
-finiteVolume/interpolation/VectorNSurfaceInterpolationSchemes.C
-
-LIB = $(FOAM_LIBBIN)/libVectorN

src/VectorN/Make/options

@@ -1,4 +0,0 @@
-EXE_INC = \
-    -I$(FOAM_SRC)/finiteVolume/lnInclude
-
-LIB_LIBS =

src/finiteVolume/Make/files

@@ -172,7 +172,19 @@ $(derivedFvPatchFields)/uniformDensityHydrostaticPressure/uniformDensityHydrosta
 $(derivedFvPatchFields)/pulseFixedValue/pulseFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/waveTransmissiveInlet/waveTransmissiveInletFvPatchFields.C
 
-
+fvPatchVectorNFields = $(fvPatchFields)/fvPatchVectorNFields
+$(fvPatchVectorNFields)/fvPatchVectorNFields.C
+$(fvPatchVectorNFields)/genericFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/calculatedFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/emptyFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/transformFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/wedgeFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/coupledFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/processorFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/cyclicFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/fixedValueFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/zeroGradientFvPatchVectorNFields.C
+$(fvPatchVectorNFields)/fixedGradientFvPatchVectorNFields.C
 
 fvsPatchFields = fields/fvsPatchFields
 $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C
@@ -195,8 +207,18 @@ $(constraintFvsPatchFields)/overlapGgi/overlapGgiFvsPatchFields.C
 $(constraintFvsPatchFields)/mixingPlane/mixingPlaneFvsPatchFields.C
 $(constraintFvsPatchFields)/regionCoupling/regionCouplingFvsPatchFields.C
 
+fvsPatchVectorNFields = $(fvsPatchFields)/fvsPatchVectorNFields
+$(fvsPatchVectorNFields)/fvsPatchVectorNFields.C
+$(fvsPatchVectorNFields)/calculatedFvsPatchVectorNFields.C
+$(fvsPatchVectorNFields)/emptyFvsPatchVectorNFields.C
+$(fvsPatchVectorNFields)/wedgeFvsPatchVectorNFields.C
+$(fvsPatchVectorNFields)/coupledFvsPatchVectorNFields.C
+$(fvsPatchVectorNFields)/processorFvsPatchVectorNFields.C
+
 fields/volFields/volFields.C
+fields/volFields/volVectorNFields.C
 fields/surfaceFields/surfaceFields.C
+fields/surfaceFields/surfaceVectorNFields.C
 
 fvMatrices/fvMatrices.C
 fvMatrices/fvScalarMatrix/fvScalarMatrix.C
@@ -224,6 +246,7 @@ $(pointVolInterpolation)/pointVolInterpolation.C
 surfaceInterpolation = interpolation/surfaceInterpolation
 $(surfaceInterpolation)/surfaceInterpolation/surfaceInterpolation.C
 $(surfaceInterpolation)/surfaceInterpolationScheme/surfaceInterpolationSchemes.C
+$(surfaceInterpolation)/VectorNSurfaceInterpolationSchemes/VectorNSurfaceInterpolationSchemes.C
 
 schemes = $(surfaceInterpolation)/schemes
 $(schemes)/linear/linear.C
@@ -315,14 +338,17 @@ $(d2dt2Schemes)/EulerD2dt2Scheme/EulerD2dt2Schemes.C
 $(d2dt2Schemes)/backwardD2dt2Scheme/backwardD2dt2Schemes.C
 
 divSchemes = finiteVolume/divSchemes
-$(divSchemes)/divScheme/divSchemes.C
+$(divSchemes)/divScheme/divSchemes.Ci
+$(divSchemes)/vectorGaussDivScheme/vectorGaussDivScheme.C
 $(divSchemes)/gaussDivScheme/gaussDivSchemes.C
 
 gradSchemes = finiteVolume/gradSchemes
 $(gradSchemes)/gradScheme/gradSchemes.C
+$(gradSchemes)/scalarGaussGrad/scalarGaussGrad.C
 $(gradSchemes)/gaussGrad/gaussGrads.C
 $(gradSchemes)/beGaussGrad/beGaussGrads.C
 $(gradSchemes)/leastSquaresGrad/leastSquaresVectors.C
+$(gradSchemes)/scalarLeastSquaresGrad/scalarLeastSquaresGrad.C
 $(gradSchemes)/leastSquaresGrad/leastSquaresGrads.C
 $(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresVectors.C
 $(gradSchemes)/extendedLeastSquaresGrad/extendedLeastSquaresGrads.C

src/finiteVolume/blockLduSystem/BlockLduSystem.C

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlockLduSystem.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class blockType, class sourceType>
+Foam::BlockLduSystem<blockType, sourceType>::BlockLduSystem
+(
+    const lduMesh& ldu
+)
+:
+    BlockLduMatrix<blockType>(ldu),
+    source_(ldu.lduAddr().size(), pTraits<sourceType>::zero)
+{}
+
+
+template<class blockType, class sourceType>
+Foam::BlockLduSystem<blockType, sourceType>::BlockLduSystem
+(
+    const lduMesh& ldu,
+    const Field<sourceType>& s
+)
+:
+    BlockLduMatrix<blockType>(ldu),
+    source_(s)
+{
+    if (ldu.lduAddr().size() != s.size())
+    {
+        FatalErrorIn
+        (
+            "BlockLduSystem::BlockLduSystem\n"
+            "(\n"
+            "    const lduMesh& ldu,"
+            "    const Field<sourceType>& s,"
+            ")\n"
+        )   << "Sizes of ldu addressing and source field are not the same."
+            << abort(FatalError);
+    }
+}
+
+
+template<class blockType, class sourceType>
+Foam::BlockLduSystem<blockType, sourceType>::BlockLduSystem
+(
+    const BlockLduMatrix<blockType>& bm,
+    const Field<sourceType>& s
+)
+:
+    BlockLduMatrix<blockType>(bm),
+    source_(s)
+{
+    if (this->lduAddr().size() != s.size())
+    {
+        FatalErrorIn
+        (
+            "BlockLduSystem::BlockLduSystem\n"
+            "(\n"
+            "    const BlockLduMatrix<blockType>& bm,"
+            "    const Field<sourceType>& s,"
+            ")\n"
+        )   << "Sizes of block matrix and source field are not the same."
+            << abort(FatalError);
+    }
+}
+
+
+template<class blockType, class sourceType>
+Foam::BlockLduSystem<blockType, sourceType>::BlockLduSystem
+(
+    const BlockLduSystem<blockType, sourceType>& bs
+)
+:
+    BlockLduMatrix<blockType>(bs),
+    source_(bs.source())
+{}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class blockType, class sourceType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const BlockLduSystem<blockType, sourceType>& bs
+)
+{
+    os  << static_cast<const BlockLduMatrix<blockType>&>(bs) << nl
+        << bs.source() << endl;
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/blockLduSystem/BlockLduSystem.H

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::BlockLduSystem
+
+Description
+    BlockLduSystem is a wrapper for BlockLduMatrix with source field. This is
+    the return type of implicit div and grad operators needed for implicitly
+    coupled systems (namely p-U coupling).
+
+Author
+    Vuko Vukcevic, FMENA Zagreb.
+
+SourceFiles
+    BlockLduSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlockLduSystem_H
+#define BlockLduSystem_H
+
+#include "blockVectorNMatrices.H"
+#include "VectorNFieldTypes.H"
+#include "volVectorNFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class blockType, class sourceType>
+class BlockLduSystem;
+
+template<class blockType, class sourceType>
+Ostream& operator<<(Ostream&, const BlockLduSystem<blockType, sourceType>&);
+
+
+/*---------------------------------------------------------------------------*\
+                        Class BlockLduSystem Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class blockType, class sourceType>
+class BlockLduSystem
+:
+    public BlockLduMatrix<blockType>
+{
+    // Private data
+
+        //- Source term
+        Field<sourceType> source_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct given addressing
+        explicit BlockLduSystem(const lduMesh&);
+
+        //- Construct given addressing and source field
+        BlockLduSystem(const lduMesh&, const Field<sourceType>&);
+
+        //- Construct from components
+        BlockLduSystem
+        (
+            const BlockLduMatrix<blockType>&,
+            const Field<sourceType>&
+        );
+
+        //- Construct as copy
+        BlockLduSystem(const BlockLduSystem<blockType, sourceType>&);
+
+
+    // Destructor
+
+       virtual ~BlockLduSystem()
+       {}
+
+
+    // Member functions
+
+        //- Access
+
+            Field<sourceType>& source()
+            {
+                return source_;
+            }
+
+            const Field<sourceType>& source() const
+            {
+                return source_;
+            }
+            
+
+    // Member operators
+
+        void operator=(const BlockLduSystem<blockType, sourceType>&);
+
+        void negate();
+
+        void operator+=(const BlockLduSystem<blockType, sourceType>&);
+        void operator-=(const BlockLduSystem<blockType, sourceType>&);
+
+        void operator*=(const scalarField&);
+        void operator*=(const scalar);
+
+
+    // Ostream operator
+
+        friend Ostream& operator<< <blockType, sourceType>
+        (
+            Ostream&,
+            const BlockLduSystem<blockType, sourceType>&
+        );
+};
+
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "BlockLduSystem.C"
+#   include "BlockLduSystemOperations.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/blockLduSystem/BlockLduSystemOperations.C

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlockLduSystem.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class blockType, class sourceType>
+void Foam::BlockLduSystem<blockType, sourceType>::negate()
+{
+    BlockLduMatrix<blockType>::negate();
+    source_.negate();
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class blockType, class sourceType>
+void Foam::BlockLduSystem<blockType, sourceType>::operator=
+(
+    const BlockLduSystem<blockType, sourceType>& bs
+)
+{
+    if (this == &bs)
+    {
+        FatalErrorIn
+        (
+            "void BlockLduSystem<blockType, sourceType>::operator="
+            "(const BlockLduSystem<blockType, sourceType>& bs)"
+        )   << "attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    BlockLduMatrix<blockType>::operator=(bs);
+    source_ = bs.source();
+}
+
+
+template<class blockType, class sourceType>
+void Foam::BlockLduSystem<blockType, sourceType>::operator+=
+(
+    const BlockLduSystem<blockType, sourceType>& bs
+)
+{
+    BlockLduMatrix<blockType>::operator+=(bs);
+    source_ += bs.source();
+}
+
+template<class blockType, class sourceType>
+void Foam::BlockLduSystem<blockType, sourceType>::operator-=
+(
+    const BlockLduSystem<blockType, sourceType>& bs
+)
+{
+    BlockLduMatrix<blockType>::operator-=(bs);
+    source_ -= bs.source();
+}
+
+template<class blockType, class sourceType>
+void Foam::BlockLduSystem<blockType, sourceType>::operator*=
+(
+    const scalarField& sf
+)
+{
+    BlockLduMatrix<blockType>::operator*=(sf);
+    source_ *= sf;
+}
+
+template<class blockType, class sourceType>
+void Foam::BlockLduSystem<blockType, sourceType>::operator*=
+(
+    const scalar s
+)
+{
+    BlockLduMatrix<blockType>::operator*=(s);
+    source_ *= s;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/blockMatrixTools/blockMatrixTools.C

@@ -862,12 +862,12 @@ void insertEquation
 }
 
 
-template<class blockType1, class blockType2>
+template<class blockType1, class fieldType, class blockType2>
 void insertBlock
 (
     const label loc1,
     const label loc2,
-    const BlockLduMatrix<blockType1>& blockMatrix,
+    const BlockLduSystem<blockType1, fieldType>& blockSystem,
     BlockLduMatrix<blockType2>& A,
     const bool incFirst
 )
@@ -875,9 +875,9 @@ void insertBlock
     // Sanity checks
     {
         const label blockMatrixSize = pTraits<blockType1>::nComponents;
-        const label blockSystemSize = pTraits<blockType2>::nComponents;
+        const label blockMatrixASize = pTraits<blockType2>::nComponents;
 
-        if (blockMatrixSize > blockSystemSize)
+        if (blockMatrixSize > blockMatrixASize)
         {
             FatalErrorIn
             (
@@ -908,7 +908,6 @@ void insertBlock
                 << "member function."
                 << abort(FatalError);
         }
-
     }
 
     const label nCmpts = pTraits<blockType1>::nComponents;
@@ -917,11 +916,11 @@ void insertBlock
 
     // Get references to ldu fields of blockMatrix always as linear
     const typename CoeffField<blockType1>::linearTypeField& bmd =
-        blockMatrix.diag().asLinear();
+        blockSystem.diag().asLinear();
     const typename CoeffField<blockType1>::linearTypeField& bmu =
-        blockMatrix.upper().asLinear();
+        blockSystem.upper().asLinear();
     const typename CoeffField<blockType1>::linearTypeField& bml =
-        blockMatrix.lower().asLinear();
+        blockSystem.lower().asLinear();
 
     // Get references to ldu fields of A matrix always as square
     typename CoeffField<blockType2>::squareTypeField& blockDiag =
@@ -936,13 +935,16 @@ void insertBlock
     {
         forAll(bmd, cellI)
         {
-            blockDiag[cellI](localLoc1, localLoc2) += bmd[cellI](cmptI);
+            blockDiag[cellI](localLoc1, localLoc2) +=
+                bmd[cellI].component(cmptI);
         }
 
         forAll(bmu, faceI)
         {
-            blockUpper[faceI](localLoc1, localLoc2) += bmu[faceI](cmptI);
-            blockLower[faceI](localLoc1, localLoc2) += bml[faceI](cmptI);
+            blockUpper[faceI](localLoc1, localLoc2) +=
+                bmu[faceI].component(cmptI);
+            blockLower[faceI](localLoc1, localLoc2) +=
+                bml[faceI].component(cmptI);
         }
 
         if (incFirst)
@@ -957,256 +959,86 @@ void insertBlock
 }
 
 
-template<class blockType1, class blockType2>
+template<class blockType1, class blockType2, class fieldType>
 void insertBoundaryContributions
 (
     const label loc1,
     const label loc2,
-    const BlockLduMatrix<blockType1>& blockMatrix,
-    GeometricField<scalar, fvPatchField, volMesh>& psi,
+    const BlockLduSystem<blockType1, fieldType>& blockSystem,
     BlockLduMatrix<blockType2>& A,
     Field<blockType2>& b,
     const bool incFirst
 )
 {
-    // Interpolation scheme for pressure gradient
-    tmp<surfaceInterpolationScheme<scalar> >
-    tinterpScheme
-    (
-        surfaceInterpolationScheme<scalar>::New
-        (
-            psi.mesh(),
-            psi.mesh().schemesDict().interpolationScheme
-            (
-                "grad(" + psi.name() + ')'
-            )
-        )
-    );
-
-    surfaceScalarField weights(tinterpScheme().weights(psi));
+    // Need to get reference to fvMesh instead of lduMesh
+    const fvBoundaryMesh& bmesh = refCast<const fvMesh>(A.mesh()).boundary();
 
     const label nCmpts = pTraits<blockType1>::nComponents;
+    const label nSrcCmpts = pTraits<fieldType>::nComponents;
     label localLoc1 = loc1;
     label localLoc2 = loc2;
 
-    // Get references to ldu fields of A matrix always as square
-    typename CoeffField<blockType2>::squareTypeField& blockDiag =
-         A.diag().asSquare();
+    const Field<fieldType>& source = blockSystem.source();
 
-    // Insert boundary contributions directly into the matrix, this assumes
-    // that blockMatrix came from implicit grad or div operators!
-    psi.boundaryField().updateCoeffs();
-    forAll(psi.boundaryField(), patchI)
+    // Insert source from block system to rhs
+    for (label cmptI = 0; cmptI < nSrcCmpts; cmptI++)
     {
-        const fvPatchScalarField& pf = psi.boundaryField()[patchI];
-        const fvPatch& patch = pf.patch();
-        const vectorField& Sf = patch.Sf();
-        const fvsPatchScalarField& pw = weights.boundaryField()[patchI];
+        scalarField sourceCmpt(source.component(cmptI));
 
-        scalarField internalCoeffs(pf.valueInternalCoeffs(pw));
-
-        if (patch.coupled())
+        forAll(b, cellI)
         {
-            typename CoeffField<blockType2>::squareTypeField& pcoupleUpper =
-                A.coupleUpper()[patchI].asSquare();
-            typename CoeffField<blockType2>::squareTypeField& pcoupleLower =
-                A.coupleLower()[patchI].asSquare();
+            b[cellI](localLoc1) += sourceCmpt[cellI];
+        }
 
-            vectorField pcl = -pw*Sf;
-            vectorField pcu = pcl + Sf;
-
-            for (label cmptI = 0; cmptI < nCmpts; cmptI++)
-            {
-                forAll(pf, faceI)
-                {
-                    label cellI = patch.faceCells()[faceI];
-
-                    // Diag contribution
-                    blockDiag[cellI](localLoc1, localLoc2) +=
-                        internalCoeffs[faceI]*Sf[faceI].component(cmptI);
-
-                    // Coupling contributions
-                    pcoupleUpper[faceI](localLoc1, localLoc2) -=
-                          pcu[faceI].component(cmptI);
-                    pcoupleLower[faceI](localLoc1, localLoc2) -=
-                          pcl[faceI].component(cmptI);
-                }
-
-                if (incFirst)
-                {
-                    localLoc1++;
-                }
-                else
-                {
-                    localLoc2++;
-                }
-            }
-
-            // Reset local locations for other patches
-            localLoc1 = loc1;
-            localLoc2 = loc2;
+        if (incFirst)
+        {
+            localLoc1++;
         }
         else
         {
-            scalarField boundaryCoeffs(pf.valueBoundaryCoeffs(pw));
-
-            for (label cmptI = 0; cmptI < nCmpts; cmptI++)
-            {
-                forAll(pf, faceI)
-                {
-                    label cellI = patch.faceCells()[faceI];
-
-                    // Diag contribution
-                    blockDiag[cellI](localLoc1, localLoc2) +=
-                        internalCoeffs[faceI]*Sf[faceI].component(cmptI);
-
-                    // Source contribution
-                    b[cellI](localLoc1) -=
-                        boundaryCoeffs[faceI]*Sf[faceI].component(cmptI);
-                }
-
-                if (incFirst)
-                {
-                    localLoc1++;
-                }
-                else
-                {
-                    localLoc2++;
-                }
-            }
-
-            // Reset local locations for other patches
-            localLoc1 = loc1;
-            localLoc2 = loc2;
+            localLoc2++;
         }
     }
-}
 
+    // Reset local locations for coupling contributions
+    localLoc1 = loc1;
+    localLoc2 = loc2;
 
-template<class blockType1, class blockType2>
-void insertBoundaryContributions
-(
-    const label loc1,
-    const label loc2,
-    const BlockLduMatrix<blockType1>& blockMatrix,
-    GeometricField<vector, fvPatchField, volMesh>& psi,
-    BlockLduMatrix<blockType2>& A,
-    Field<blockType2>& b,
-    const bool incFirst
-)
-{
-    // Interpolation scheme for velocity divergence
-    tmp<surfaceInterpolationScheme<scalar> >
-    tinterpScheme
-    (
-        surfaceInterpolationScheme<scalar>::New
-        (
-            psi.mesh(),
-            psi.mesh().schemesDict().interpolationScheme
-            (
-                "div(" + psi.name() + ')'
-            )
-        )
-    );
-
-    surfaceScalarField weights(tinterpScheme().weights(mag(psi)));
-
-    const label nCmpts = pTraits<blockType1>::nComponents;
-    label localLoc1 = loc1;
-    label localLoc2 = loc2;
-
-    // Get references to ldu fields of A matrix always as square
-    typename CoeffField<blockType2>::squareTypeField& blockDiag =
-         A.diag().asSquare();
-
-    // Insert boundary contributions directly into the matrix, this assumes
-    // that blockMatrix came from implicit grad or div operators!
-    psi.boundaryField().updateCoeffs();
-    forAll(psi.boundaryField(), patchI)
+    // Insert coupling contributions into block matrix
+    forAll(bmesh, patchI)
     {
-        const fvPatchVectorField& pf = psi.boundaryField()[patchI];
-        const fvPatch& patch = pf.patch();
-        const vectorField& Sf = patch.Sf();
-        const fvsPatchScalarField& pw = weights.boundaryField()[patchI];
-
-        vectorField internalCoeffs(pf.valueInternalCoeffs(pw));
-
-        if (patch.coupled())
+        if (bmesh[patchI].coupled())
         {
             typename CoeffField<blockType2>::squareTypeField& pcoupleUpper =
                 A.coupleUpper()[patchI].asSquare();
             typename CoeffField<blockType2>::squareTypeField& pcoupleLower =
                 A.coupleLower()[patchI].asSquare();
 
-            vectorField pcl = -pw*Sf;
-            vectorField pcu = pcl + Sf;
+            const typename CoeffField<blockType1>::linearTypeField& bmcu =
+                blockSystem.coupleUpper()[patchI].asLinear();
+            const typename CoeffField<blockType1>::linearTypeField& bmcl =
+                blockSystem.coupleLower()[patchI].asLinear();
 
-            for (label cmptI = 0; cmptI < nCmpts; cmptI++)
-            {
-                forAll(pf, faceI)
+                for (label cmptI = 0; cmptI < nCmpts; cmptI++)
                 {
-                    label cellI = patch.faceCells()[faceI];
+                    forAll(bmcu, faceI)
+                    {
+                        pcoupleUpper[faceI](localLoc1, localLoc2) +=
+                              bmcu[faceI].component(cmptI);
+                        pcoupleLower[faceI](localLoc1, localLoc2) +=
+                              bmcl[faceI].component(cmptI);
+                    }
 
-                    // Diag contribution
-                    blockDiag[cellI](localLoc1, localLoc2) += cmptMultiply
-                    (
-                        internalCoeffs[faceI], Sf[faceI]
-                    ).component(cmptI);
-
-                    // Coupling contributions
-                    pcoupleLower[faceI](localLoc1, localLoc2) -=
-                        pcl[faceI].component(cmptI);
-                    pcoupleUpper[faceI](localLoc1, localLoc2) -=
-                        pcu[faceI].component(cmptI);
+                    if (incFirst)
+                    {
+                        localLoc1++;
+                    }
+                    else
+                    {
+                        localLoc2++;
+                    }
                 }
 
-                if (incFirst)
-                {
-                    localLoc1++;
-                }
-                else
-                {
-                    localLoc2++;
-                }
-            }
-
-            // Reset local locations for other patches
-            localLoc1 = loc1;
-            localLoc2 = loc2;
-        }
-        else
-        {
-            vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw));
-
-            for (label cmptI = 0; cmptI < nCmpts; cmptI++)
-            {
-                forAll(pf, faceI)
-                {
-                    label cellI = patch.faceCells()[faceI];
-
-                    // Diag contribution
-                    blockDiag[cellI](localLoc1, localLoc2) += cmptMultiply
-                    (
-                        internalCoeffs[faceI], Sf[faceI]
-                    ).component(cmptI);
-
-                    // Source contribution
-                    b[cellI](localLoc1) -= cmptMultiply
-                    (
-                        boundaryCoeffs[faceI], Sf[faceI]
-                    ).component(cmptI);
-                }
-
-                if (incFirst)
-                {
-                    localLoc1++;
-                }
-                else
-                {
-                    localLoc2++;
-                }
-            }
-
             // Reset local locations for other patches
             localLoc1 = loc1;
             localLoc2 = loc2;
@@ -1220,15 +1052,14 @@ void insertBlockCoupling
 (
     const label loc1,
     const label loc2,
-    const BlockLduMatrix<blockType1>& blockMatrix,
-    GeometricField<fieldType, fvPatchField, volMesh>& psi,
+    const BlockLduSystem<blockType1, fieldType>& blockSystem,
     BlockLduMatrix<blockType2>& A,
     Field<blockType2>& b,
     const bool incFirst
 )
 {
-    insertBlock(loc1, loc2, blockMatrix, A, incFirst);
-    insertBoundaryContributions(loc1, loc2, blockMatrix, psi, A, b, incFirst);
+    insertBlock(loc1, loc2, blockSystem, A, incFirst);
+    insertBoundaryContributions(loc1, loc2, blockSystem, A, b, incFirst);
 }
 
 
@@ -1262,7 +1093,7 @@ void insertEquationCoupling
     }
 
     // Get references to fvScalarMatrix fields, updating boundary contributions
-    const scalarField& diag = matrix.D()();
+    scalarField& diag = matrix.D();
     scalarField& source = matrix.source();
     matrix.addBoundarySource(source, false);
 

src/finiteVolume/blockMatrixTools/blockMatrixTools.H

@@ -25,7 +25,11 @@ InNamespace
     Foam::
 
 Description
-    Block matrix insertion and retrieval tools
+    Block matrix insertion and retrieval tools.
+    Note: these will be obsolete since all of the functions are now member
+    functions in fvBlockMatrix class. All top level coupled solvers should
+    be rewritten using fvBlockMatrix - then we can get rid of this.
+    VV, 20/July/2014.
 
 SourceFiles
     blockMatrixTools.C
@@ -36,7 +40,7 @@ SourceFiles
 #ifndef blockMatrixTools_H
 #define blockMatrixTools_H
 
-#include "blockLduMatrices.H"
+#include "BlockLduSystem.H"
 #include "fvMatrices.H"
 #include "surfaceFieldsFwd.H"
 
@@ -218,60 +222,37 @@ namespace blockMatrixTools
             Field<blockType>& b
         );
 
-        // Insert block coupling
-        template<class blockType1, class blockType2>
+        // Insert diagonal, lower and upper of blockMatrix into A
+        template<class blockType1, class fieldType, class blockType2>
         void insertBlock
         (
             const label loc1,
             const label loc2,
-            const BlockLduMatrix<blockType1>& blockMatrix,
+            const BlockLduSystem<blockType1, fieldType>& blockMatrix,
             BlockLduMatrix<blockType2>& A,
             const bool incFirst
         );
 
-        // Insert boundary contributions. These functions update the matrix coeffs
-        // according to boundary conditions. Two functions are needed for
-        // specialization if the field is scalar or vector. If the field is scalar,
-        // it is assumed that the grad(p) was used to obtain the matrix, and the
-        // boundary contributions to the source are vectors. If the field is
-        // vector, assumed is the div(U) and the source is scalar.
-        // Note parameter GeometricField<scalar...>, this is grad specialization.
-        template<class blockType1, class blockType2>
+        // Insert source and coupling coeffs of block system into A and b
+        template<class blockType1, class blockType2, class fieldType>
         void insertBoundaryContributions
         (
             const label loc1,
             const label loc2,
-            const BlockLduMatrix<blockType1>& blockMatrix,
-            GeometricField<scalar, fvPatchField, volMesh>& psi,
+            const BlockLduSystem<blockType1, fieldType>& blockSystem,
             BlockLduMatrix<blockType2>& A,
             Field<blockType2>& b,
             const bool incFirst
         );
 
-        // Note parameter GeometricField<vector...>, this is div specialization.
-        template<class blockType1, class blockType2>
-        void insertBoundaryContributions
-        (
-            const label loc1,
-            const label loc2,
-            const BlockLduMatrix<blockType1>& blockMatrix,
-            GeometricField<vector, fvPatchField, volMesh>& psi,
-            BlockLduMatrix<blockType2>& A,
-            Field<blockType2>& b,
-            const bool incFirst
-        );
-
-        // Insert existing block matrix (obtained by implicit grad/div operator)
-        // into block system, updating the matrix coefficients from boundaries.
-        // To be used ONLY for pressure velocity coupling terms, other scenarios
-        // are not considered at the moment.
+        // Insert existing block system (obtained by implicit grad/div operator)
+        // into block system.
         template<class blockType1, class blockType2, class fieldType>
         void insertBlockCoupling
         (
             const label loc1,
             const label loc2,
-            const BlockLduMatrix<blockType1>& blockMatrix,
-            GeometricField<fieldType, fvPatchField, volMesh>& psi,
+            const BlockLduSystem<blockType1, fieldType>& blockSystem,
             BlockLduMatrix<blockType2>& A,
             Field<blockType2>& b,
             const bool incFirst

src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C

@@ -98,6 +98,35 @@ divScheme<Type>::~divScheme()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+>
+divScheme<Type>::fvmDiv
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    FatalErrorIn
+    (
+        "tmp<BlockLduSystem> divScheme<Type>::fvmDiv\n"
+        "(\n"
+        "    GeometricField<Type, fvPatchField, volMesh>&"
+        ")\n"
+    )   << "Implicit div operator currently defined only for Gauss linear. "
+        << abort(FatalError);
+
+    typedef typename innerProduct<vector, Type>::type DivType;
+
+    tmp<BlockLduSystem<vector, DivType> > tbs
+    (
+        new BlockLduSystem<vector, DivType>(vf.mesh())
+    );
+
+    return tbs;
+}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H

@@ -41,7 +41,7 @@ SourceFiles
 #include "linear.H"
 #include "typeInfo.H"
 #include "runTimeSelectionTables.H"
-#include "blockLduMatrices.H"
+#include "BlockLduSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -152,12 +152,15 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) = 0;
 
-        //- Return the blockLduMatrix<vector> corresponding to the implicit
-        // div discretization. For block coupled system.
-        virtual tmp<blockVectorMatrix> fvmDiv
+        //- Return the BlockLduSystem corresponding to the implicit div
+        // discretization. For block coupled system.
+        virtual tmp
+        <
+            BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+        > fvmDiv
         (
             const GeometricField<Type, fvPatchField, volMesh>&
-        ) const = 0;
+        ) const;
 };
 
 

src/finiteVolume/finiteVolume/divSchemes/divScheme/divSchemes.C

@@ -61,13 +61,13 @@ defineTemplateRunTimeSelectionTable
 
 defineTemplateRunTimeSelectionTable
 (
- divScheme<symmTensor4thOrder>,
+    divScheme<symmTensor4thOrder>,
     Istream
  );
 
 defineTemplateRunTimeSelectionTable
 (
- divScheme<diagTensor>,
+    divScheme<diagTensor>,
     Istream
  );
 

src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C

@@ -69,38 +69,31 @@ gaussDivScheme<Type>::fvcDiv
 
 
 template<class Type>
-tmp<blockVectorMatrix> gaussDivScheme<Type>::fvmDiv
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> gaussDivScheme<Type>::fvmDiv
 (
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
-    tmp<surfaceScalarField> tweights = this->tinterpScheme_().weights(vf);
-    const scalarField& wIn = tweights().internalField();
+   FatalErrorIn
+   (
+       "tmp<BlockLduSystem> fvmDiv\n"
+       "(\n"
+       "    GeometricField<Type, fvPatchField, volMesh>&"
+       ")\n"
+   )   << "Implicit div operator defined only for vector."
+       << abort(FatalError);
 
-    const fvMesh& mesh = vf.mesh();
+   typedef typename innerProduct<vector, Type>::type DivType;
 
-    tmp<blockVectorMatrix> tbm
-    (
-        new blockVectorMatrix
-        (
-           mesh
-        )
-    );
-    blockVectorMatrix& bm = tbm();
+   tmp<BlockLduSystem<vector, DivType> > tbs
+   (
+       new BlockLduSystem<vector, DivType>(vf.mesh())
+   );
 
-    // Grab ldu parts of block matrix as linear always
-    typename CoeffField<vector>::linearTypeField& u = bm.upper().asLinear();
-    typename CoeffField<vector>::linearTypeField& l = bm.lower().asLinear();
-
-    const vectorField& SfIn = mesh.Sf().internalField();
-
-    l = -wIn*SfIn;
-    u = l + SfIn;
-    bm.negSumDiag();
-
-    // Interpolation schemes with corrections not accounted for
-
-    return tbm;
+   return tbs;
 }
 
 

src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.H

@@ -97,9 +97,12 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         );
 
-        //- Return the blockLduMatrix<vector> corresponding to the implicit
-        // div discretization. For block coupled system.
-        tmp<blockVectorMatrix> fvmDiv
+        //- Return the BlockLduSystem corresponding to the implicit div
+        // discretization. For block coupled system.
+        tmp
+        <
+            BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+        > fvmDiv
         (
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const;

src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.C

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Specialisation of gaussDivScheme for vectors. Needed for implicit fvmDiv
+    operator for block coupled systems.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vectorGaussDivScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
+(
+    const GeometricField<vector, fvPatchField, volMesh>& vf
+) const
+{
+    tmp<surfaceScalarField> tweights = this->tinterpScheme_().weights(vf);
+    const scalarField& wIn = tweights().internalField();
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<BlockLduSystem<vector, scalar> > tbs
+    (
+        new BlockLduSystem<vector, scalar>(mesh)
+    );
+    BlockLduSystem<vector, scalar>& bs = tbs();
+    scalarField& source = bs.source();
+
+    // Grab ldu parts of block matrix as linear always
+    typename CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
+    typename CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
+    typename CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
+
+    const vectorField& SfIn = mesh.Sf().internalField();
+
+    l = -wIn*SfIn;
+    u = l + SfIn;
+    bs.negSumDiag();
+
+    // Boundary contributions
+    forAll(vf.boundaryField(), patchI)
+    {
+        const fvPatchVectorField& pf = vf.boundaryField()[patchI];
+        const fvPatch& patch = pf.patch();
+        const vectorField& Sf = patch.Sf();
+        const fvsPatchScalarField& pw = tweights().boundaryField()[patchI];
+        const unallocLabelList& fc = patch.faceCells();
+
+        const vectorField internalCoeffs(pf.valueInternalCoeffs(pw));
+
+        // Diag contribution
+        forAll(pf, faceI)
+        {
+            d[fc[faceI]] += cmptMultiply(internalCoeffs[faceI], Sf[faceI]);
+        }
+
+        if (patch.coupled())
+        {
+            typename CoeffField<vector>::linearTypeField& pcoupleUpper =
+                bs.coupleUpper()[patchI].asLinear();
+            typename CoeffField<vector>::linearTypeField& pcoupleLower =
+                bs.coupleLower()[patchI].asLinear();
+
+            const vectorField pcl = -pw*Sf;
+            const vectorField pcu = pcl + Sf;
+
+            // Coupling  contributions
+            pcoupleLower -= pcl;
+            pcoupleUpper -= pcu;
+        }
+        else
+        {
+            const vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw));
+
+            // Boundary contribution
+            forAll(pf, faceI)
+            {
+                source[fc[faceI]] -= boundaryCoeffs[faceI] & Sf[faceI];
+            }
+        }
+    }
+
+    // Interpolation schemes with corrections not accounted for
+
+    return tbs;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/divSchemes/vectorGaussDivScheme/vectorGaussDivScheme.H

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::fv::vectorGaussDivScheme
+
+Description
+    Specialisation of gaussDivScheme for vectors. Needed for implicit fvmDiv
+    operator for block coupled systems.
+
+SourceFiles
+    vectorGaussDivScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef vectorGaussDivScheme_H
+#define vectorGaussDivScheme_H
+
+#include "gaussDivScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmDiv
+(
+    const GeometricField<vector, fvPatchField, volMesh>&
+) const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/fvm/fvmDiv.C

@@ -98,8 +98,12 @@ div
     return Div;
 }
 
+
 template<class Type>
-tmp<BlockLduMatrix<vector> > div
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> div
 (
     GeometricField<Type, fvPatchField, volMesh>& vf,
     const word& name
@@ -114,7 +118,10 @@ tmp<BlockLduMatrix<vector> > div
 
 
 template<class Type>
-tmp<BlockLduMatrix<vector> > div
+tmp
+<
+    BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+> div
 (
     GeometricField<Type, fvPatchField, volMesh>& vf
 )

src/finiteVolume/finiteVolume/fvm/fvmDiv.H

@@ -39,7 +39,7 @@ SourceFiles
 #include "surfaceFieldsFwd.H"
 #include "surfaceInterpolationScheme.H"
 #include "fvMatrices.H"
-#include "blockLduMatrices.H"
+#include "BlockLduSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -85,14 +85,20 @@ namespace fvm
 
     // Implicit div operator for block systems
     template<class Type>
-    tmp<blockVectorMatrix> div
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    > div
     (
         GeometricField<Type, fvPatchField, volMesh>&,
         const word&
     );
 
     template<class Type>
-    tmp<blockVectorMatrix> div
+    tmp
+    <
+        BlockLduSystem<vector, typename innerProduct<vector, Type>::type>
+    > div
     (
         GeometricField<Type, fvPatchField, volMesh>&
     );