From c13421b10a5fcd2808d62e3556866237c6815329 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sun, 20 Mar 2016 19:44:29 +0000
Subject: [PATCH] Matrix: Added (i, j) addressing to allow support for
 addressing blocks of the matrix

This change brings OpenFOAM into line with the standard matrix
addressing in other C++ libraries for better interoperability.
---
 applications/test/Matrix/Test-Matrix.C        | 54 ++++++++++------
 applications/test/ODE/Test-ODE.C              | 34 +++++-----
 .../test/simpleMatrix/Test-simpleMatrix.C     | 20 +++---
 .../surfaceFeatureExtract.C                   | 10 +--
 src/ODE/ODESolvers/EulerSI/EulerSI.C          |  6 +-
 .../ODESolvers/Rosenbrock12/Rosenbrock12.C    |  6 +-
 .../ODESolvers/Rosenbrock23/Rosenbrock23.C    |  6 +-
 .../ODESolvers/Rosenbrock34/Rosenbrock34.C    |  6 +-
 src/ODE/ODESolvers/SIBS/SIMPR.C               |  6 +-
 src/ODE/ODESolvers/SIBS/polyExtrapolate.C     |  4 +-
 src/ODE/ODESolvers/rodas23/rodas23.C          |  6 +-
 src/ODE/ODESolvers/rodas34/rodas34.C          |  6 +-
 src/ODE/ODESolvers/seulex/seulex.C            | 14 ++---
 .../matrices/DiagonalMatrix/DiagonalMatrix.C  |  2 +-
 src/OpenFOAM/matrices/Matrix/Matrix.C         |  8 +--
 src/OpenFOAM/matrices/Matrix/MatrixI.H        |  4 +-
 .../matrices/SquareMatrix/SquareMatrix.C      |  2 +-
 .../SymmetricSquareMatrix.C                   |  8 +--
 .../SymmetricSquareMatrix.H                   |  7 ---
 .../SymmetricSquareMatrixI.H                  | 36 -----------
 .../matrices/scalarMatrices/SVD/SVD.C         | 62 +++++++++----------
 .../matrices/scalarMatrices/scalarMatrices.C  | 26 ++++----
 .../scalarMatrices/scalarMatricesTemplates.C  | 18 +++---
 .../CentredFitSnGrad/CentredFitSnGradData.C   | 26 ++++----
 .../schemes/FitData/FitData.C                 | 14 ++---
 .../chemistryModel/chemistryModel.C           |  2 +-
 .../EulerImplicit/EulerImplicit.C             |  6 +-
 .../radiationModels/viewFactor/viewFactor.C   | 20 +++---
 .../pyrolysisChemistryModel.C                 |  2 +-
 29 files changed, 197 insertions(+), 224 deletions(-)

diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index 380576d69..3e14e5708 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,15 +36,15 @@ int main(int argc, char *argv[])
 {
     SquareMatrix<scalar> hmm(3);
 
-    hmm[0][0] = -3.0;
-    hmm[0][1] = 10.0;
-    hmm[0][2] = -4.0;
-    hmm[1][0] = 2.0;
-    hmm[1][1] = 3.0;
-    hmm[1][2] = 10.0;
-    hmm[2][0] = 2.0;
-    hmm[2][1] = 6.0;
-    hmm[2][2] = 1.0;
+    hmm(0, 0) = -3.0;
+    hmm(0, 1) = 10.0;
+    hmm(0, 2) = -4.0;
+    hmm(1, 0) = 2.0;
+    hmm(1, 1) = 3.0;
+    hmm(1, 2) = 10.0;
+    hmm(2, 0) = 2.0;
+    hmm(2, 1) = 6.0;
+    hmm(2, 2) = 1.0;
 
     //Info<< hmm << endl << hmm - 2.0*(-hmm) << endl;
     Info<< max(hmm) << endl;
@@ -106,15 +106,15 @@ int main(int argc, char *argv[])
     {
         scalarSquareMatrix squareMatrix(3, 3, 0);
 
-        squareMatrix[0][0] = 4;
-        squareMatrix[0][1] = 12;
-        squareMatrix[0][2] = -16;
-        squareMatrix[1][0] = 12;
-        squareMatrix[1][1] = 37;
-        squareMatrix[1][2] = -43;
-        squareMatrix[2][0] = -16;
-        squareMatrix[2][1] = -43;
-        squareMatrix[2][2] = 98;
+        squareMatrix(0, 0) = 4;
+        squareMatrix(0, 1) = 12;
+        squareMatrix(0, 2) = -16;
+        squareMatrix(1, 0) = 12;
+        squareMatrix(1, 1) = 37;
+        squareMatrix(1, 2) = -43;
+        squareMatrix(2, 0) = -16;
+        squareMatrix(2, 1) = -43;
+        squareMatrix(2, 2) = 98;
 
         const scalarSquareMatrix squareMatrixCopy = squareMatrix;
         Info<< nl << "Square Matrix = " << squareMatrix << endl;
@@ -131,6 +131,22 @@ int main(int argc, char *argv[])
         Info<< "det = " << detDecomposed(squareMatrix, sign) << endl;
     }
 
+    {
+        scalarSquareMatrix squareMatrix(3000, 3000, 1);
+
+        for(label i=0; i<squareMatrix.n(); i++)
+        {
+            squareMatrix(i, i) = 10;
+        }
+
+        scalarField rhs(squareMatrix.n(), 0);
+        rhs[0] = 1;
+        rhs[1] = 2;
+        rhs[2] = 3;
+
+        LUsolve(squareMatrix, rhs);
+    }
+
     Info<< "\nEnd\n" << endl;
 
     return 0;
diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C
index ae2037425..689c634fd 100644
--- a/applications/test/ODE/Test-ODE.C
+++ b/applications/test/ODE/Test-ODE.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,25 +75,25 @@ public:
         dfdx[2] = (2.0/sqr(x))*y[2];
         dfdx[3] = (3.0/sqr(x))*y[3];
 
-        dfdy[0][0] = 0.0;
-        dfdy[0][1] = -1.0;
-        dfdy[0][2] = 0.0;
-        dfdy[0][3] = 0.0;
+        dfdy(0, 0) = 0.0;
+        dfdy(0, 1) = -1.0;
+        dfdy(0, 2) = 0.0;
+        dfdy(0, 3) = 0.0;
 
-        dfdy[1][0] = 1.0;
-        dfdy[1][1] = -1.0/x;
-        dfdy[1][2] = 0.0;
-        dfdy[1][3] = 0.0;
+        dfdy(1, 0) = 1.0;
+        dfdy(1, 1) = -1.0/x;
+        dfdy(1, 2) = 0.0;
+        dfdy(1, 3) = 0.0;
 
-        dfdy[2][0] = 0.0;
-        dfdy[2][1] = 1.0;
-        dfdy[2][2] = -2.0/x;
-        dfdy[2][3] = 0.0;
+        dfdy(2, 0) = 0.0;
+        dfdy(2, 1) = 1.0;
+        dfdy(2, 2) = -2.0/x;
+        dfdy(2, 3) = 0.0;
 
-        dfdy[3][0] = 0.0;
-        dfdy[3][1] = 0.0;
-        dfdy[3][2] = 1.0;
-        dfdy[3][3] = -3.0/x;
+        dfdy(3, 0) = 0.0;
+        dfdy(3, 1) = 0.0;
+        dfdy(3, 2) = 1.0;
+        dfdy(3, 3) = -3.0/x;
     }
 };
 
diff --git a/applications/test/simpleMatrix/Test-simpleMatrix.C b/applications/test/simpleMatrix/Test-simpleMatrix.C
index e5010212e..d67b9d3bd 100644
--- a/applications/test/simpleMatrix/Test-simpleMatrix.C
+++ b/applications/test/simpleMatrix/Test-simpleMatrix.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,15 +35,15 @@ int main(int argc, char *argv[])
 {
     simpleMatrix<vector> hmm(3);
 
-    hmm[0][0] = -3.0;
-    hmm[0][1] = 10.0;
-    hmm[0][2] = -4.0;
-    hmm[1][0] = 2.0;
-    hmm[1][1] = 3.0;
-    hmm[1][2] = 10.0;
-    hmm[2][0] = 2.0;
-    hmm[2][1] = 6.0;
-    hmm[2][2] = 1.0;
+    hmm(0, 0) = -3.0;
+    hmm(0, 1) = 10.0;
+    hmm(0, 2) = -4.0;
+    hmm(1, 0) = 2.0;
+    hmm(1, 1) = 3.0;
+    hmm(1, 2) = 10.0;
+    hmm(2, 0) = 2.0;
+    hmm(2, 1) = 6.0;
+    hmm(2, 2) = 1.0;
 
     hmm.source()[0] = vector(2.0, 1.0, 3.0);
     hmm.source()[1] = vector(1.0, 4.0, 3.0);
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 542c9ba4a..510305f28 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -312,11 +312,11 @@ triSurfacePointScalarField calcCurvature
             scalar x = edgeVectors[i] & faceCoordSys[0];
             scalar y = edgeVectors[i] & faceCoordSys[1];
 
-            T[0][0] += sqr(x);
-            T[1][0] += x*y;
-            T[1][1] += sqr(x) + sqr(y);
-            T[2][1] += x*y;
-            T[2][2] += sqr(y);
+            T(0, 0) += sqr(x);
+            T(1, 0) += x*y;
+            T(1, 1) += sqr(x) + sqr(y);
+            T(2, 1) += x*y;
+            T(2, 2) += sqr(y);
 
             scalar dndx = normalDifferences[i] & faceCoordSys[0];
             scalar dndy = normalDifferences[i] & faceCoordSys[1];
diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.C b/src/ODE/ODESolvers/EulerSI/EulerSI.C
index 91a60461a..0648b2664 100644
--- a/src/ODE/ODESolvers/EulerSI/EulerSI.C
+++ b/src/ODE/ODESolvers/EulerSI/EulerSI.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,10 +67,10 @@ Foam::scalar Foam::EulerSI::solve
     {
         for (label j=0; j<n_; j++)
         {
-            a_[i][j] = -dfdy_[i][j];
+            a_(i, j) = -dfdy_(i, j);
         }
 
-        a_[i][i] += 1.0/dx;
+        a_(i, i) += 1.0/dx;
     }
 
     LUDecompose(a_, pivotIndices_);
diff --git a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C
index f5fd2b94b..fa163aec0 100644
--- a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C
+++ b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,10 +81,10 @@ Foam::scalar Foam::Rosenbrock12::solve
     {
         for (label j=0; j<n_; j++)
         {
-            a_[i][j] = -dfdy_[i][j];
+            a_(i, j) = -dfdy_(i, j);
         }
 
-        a_[i][i] += 1.0/(gamma*dx);
+        a_(i, i) += 1.0/(gamma*dx);
     }
 
     LUDecompose(a_, pivotIndices_);
diff --git a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C
index 495f351e8..3b6c8cf46 100644
--- a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C
+++ b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -94,10 +94,10 @@ Foam::scalar Foam::Rosenbrock23::solve
     {
         for (label j=0; j<n_; j++)
         {
-            a_[i][j] = -dfdy_[i][j];
+            a_(i, j) = -dfdy_(i, j);
         }
 
-        a_[i][i] += 1.0/(gamma*dx);
+        a_(i, i) += 1.0/(gamma*dx);
     }
 
     LUDecompose(a_, pivotIndices_);
diff --git a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C
index db16f6195..5b761b267 100644
--- a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C
+++ b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -139,10 +139,10 @@ Foam::scalar Foam::Rosenbrock34::solve
     {
         for (label j=0; j<n_; j++)
         {
-            a_[i][j] = -dfdy_[i][j];
+            a_(i, j) = -dfdy_(i, j);
         }
 
-        a_[i][i] += 1.0/(gamma*dx);
+        a_(i, i) += 1.0/(gamma*dx);
     }
 
     LUDecompose(a_, pivotIndices_);
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index 0017d4c23..d747b36dc 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/src/ODE/ODESolvers/SIBS/SIMPR.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,9 +46,9 @@ void Foam::SIBS::SIMPR
     {
         for (label j=0; j<n_; j++)
         {
-            a[i][j] = -h*dfdy[i][j];
+            a(i, j) = -h*dfdy(i, j);
         }
-        ++a[i][i];
+        ++a(i, i);
     }
 
     labelList pivotIndices(n_);
diff --git a/src/ODE/ODESolvers/SIBS/polyExtrapolate.C b/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
index 8b1c72d09..3f74d2091 100644
--- a/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
+++ b/src/ODE/ODESolvers/SIBS/polyExtrapolate.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,7 @@ void Foam::SIBS::polyExtrapolate
     {
         for (label j=0; j<n; j++)
         {
-            d[j][0] = yest[j];
+            d(j, 0) = yest[j];
         }
     }
     else
diff --git a/src/ODE/ODESolvers/rodas23/rodas23.C b/src/ODE/ODESolvers/rodas23/rodas23.C
index af4540b96..31d5e8ace 100644
--- a/src/ODE/ODESolvers/rodas23/rodas23.C
+++ b/src/ODE/ODESolvers/rodas23/rodas23.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,10 +85,10 @@ Foam::scalar Foam::rodas23::solve
     {
         for (label j=0; j<n_; j++)
         {
-            a_[i][j] = -dfdy_[i][j];
+            a_(i, j) = -dfdy_(i, j);
         }
 
-        a_[i][i] += 1.0/(gamma*dx);
+        a_(i, i) += 1.0/(gamma*dx);
     }
 
     LUDecompose(a_, pivotIndices_);
diff --git a/src/ODE/ODESolvers/rodas34/rodas34.C b/src/ODE/ODESolvers/rodas34/rodas34.C
index 6a15ad191..5f403238d 100644
--- a/src/ODE/ODESolvers/rodas34/rodas34.C
+++ b/src/ODE/ODESolvers/rodas34/rodas34.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -108,10 +108,10 @@ Foam::scalar Foam::rodas34::solve
     {
         for (label j=0; j<n_; j++)
         {
-            a_[i][j] = -dfdy_[i][j];
+            a_(i, j) = -dfdy_(i, j);
         }
 
-        a_[i][i] += 1.0/(gamma*dx);
+        a_(i, i) += 1.0/(gamma*dx);
     }
 
     LUDecompose(a_, pivotIndices_);
diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C
index 51389d077..f91a55557 100644
--- a/src/ODE/ODESolvers/seulex/seulex.C
+++ b/src/ODE/ODESolvers/seulex/seulex.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,7 +92,7 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
         for (int l=0; l<k; l++)
         {
             scalar ratio = scalar(nSeq_[k])/nSeq_[l];
-            coeff_[k][l] = 1.0/(ratio - 1.0);
+            coeff_(k, l) = 1.0/(ratio - 1.0);
         }
     }
 }
@@ -117,10 +117,10 @@ bool Foam::seulex::seul
     {
         for (label j=0; j<n_; j++)
         {
-            a_[i][j] = -dfdy_[i][j];
+            a_(i, j) = -dfdy_(i, j);
         }
 
-        a_[i][i] += 1.0/dx;
+        a_(i, i) += 1.0/dx;
     }
 
     LUDecompose(a_, pivotIndices_);
@@ -192,13 +192,13 @@ void Foam::seulex::extrapolate
         for (label i=0; i<n_; i++)
         {
             table[j-1][i] =
-                table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
+                table(j, i) + coeff_(k, j)*(table(j, i) - table[j-1][i]);
         }
     }
 
     for (int i=0; i<n_; i++)
     {
-        y[i] = table[0][i] + coeff_[k][0]*(table[0][i] - y[i]);
+        y[i] = table(0, i) + coeff_(k, 0)*(table(0, i) - y[i]);
     }
 }
 
@@ -288,7 +288,7 @@ void Foam::seulex::solve
                 forAll(scale_, i)
                 {
                     scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
-                    err += sqr((y[i] - table_[0][i])/scale_[i]);
+                    err += sqr((y[i] - table_(0, i))/scale_[i]);
                 }
                 err = sqrt(err/n_);
                 if (err > 1.0/SMALL || (k > 1 && err >= errOld))
diff --git a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
index 69844cd40..4345ed240 100644
--- a/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
+++ b/src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
@@ -42,7 +42,7 @@ Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a)
 {
     forAll(*this, i)
     {
-        this->operator[](i) = a[i][i];
+        this->operator[](i) = a(i, i);
     }
 }
 
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 829156693..6ccd1c79e 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -156,7 +156,7 @@ Form Foam::Matrix<Form, Type>::T() const
     {
         for (label j=0; j<n(); j++)
         {
-            At[j][i] = A[i][j];
+            At(j, i) = A(i, j);
         }
     }
 
@@ -238,7 +238,7 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
             << abort(FatalError);
 
         // Return in error to keep compiler happy
-        return a[0][0];
+        return a(0, 0);
     }
 }
 
@@ -270,7 +270,7 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
             << abort(FatalError);
 
         // Return in error to keep compiler happy
-        return a[0][0];
+        return a(0, 0);
     }
 }
 
@@ -411,7 +411,7 @@ Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
         {
             for (label l = 0; l < b.m(); l++)
             {
-                ab[i][j] += a[i][l]*b[l][j];
+                ab(i, j) += a(i, l)*b(l, j);
             }
         }
     }
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index b9ecb4c74..532e9df3f 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -136,7 +136,7 @@ inline const Type& Foam::Matrix<Form, Type>::operator()
 ) const
 {
     checki(i);
-    checki(j);
+    checkj(j);
     return v_[i*nCols_ + j];
 }
 
@@ -149,7 +149,7 @@ inline Type& Foam::Matrix<Form, Type>::operator()
 )
 {
     checki(i);
-    checki(j);
+    checkj(j);
     return v_[i*nCols_ + j];
 }
 
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
index e8e7dd0e8..74439c15d 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
@@ -39,7 +39,7 @@ Foam::scalar Foam::detDecomposed
 
     for (label i = 0; i < matrix.m(); ++i)
     {
-        diagProduct *= matrix[i][i];
+        diagProduct *= matrix(i, i);
     }
 
     return sign*diagProduct;
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
index 1ab99b124..28ad1d274 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
@@ -37,7 +37,7 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
 
     for (label i = 0; i < matrix.m(); ++i)
     {
-        inv[i][i] = 1.0/matrix[i][i];
+        inv(i, i) = 1.0/matrix(i, i);
 
         for (label j = 0; j < i; ++j)
         {
@@ -45,10 +45,10 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
 
             for (label k = j; k < i; k++)
             {
-                sum -= matrix[i][k]*inv[k][j];
+                sum -= matrix(i, k)*inv(k, j);
             }
 
-            inv[i][j] = sum/matrix[i][i];
+            inv(i, j) = sum/matrix(i, i);
         }
     }
 
@@ -77,7 +77,7 @@ Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
 
     for (label i = 0; i < matrix.m(); ++i)
     {
-        diagProduct *= matrix[i][i];
+        diagProduct *= matrix(i, i);
     }
 
     return sqr(diagProduct);
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
index a5050b24c..193ee1602 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H
@@ -76,13 +76,6 @@ public:
 
         //- Clone
         inline autoPtr<SymmetricSquareMatrix<Type>> clone() const;
-
-
-        //- Return subscript-checked row of Matrix.
-        inline Type& operator()(const label r, const label c);
-
-        //- Return subscript-checked row of constant Matrix.
-        inline const Type& operator()(const label r, const label c) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
index 00757eabc..a363026aa 100644
--- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrixI.H
@@ -94,40 +94,4 @@ Foam::SymmetricSquareMatrix<Type>::clone() const
 }
 
 
-template<class Type>
-inline Type& Foam::SymmetricSquareMatrix<Type>::operator()
-(
-    const label r,
-    const label c
-)
-{
-    if (r > c)
-    {
-        return this->operator[](r)[c];
-    }
-    else
-    {
-        return this->operator[](c)[r];
-    }
-}
-
-
-template<class Type>
-inline const Type& Foam::SymmetricSquareMatrix<Type>::operator()
-(
-    const label r,
-    const label c
-) const
-{
-    if (r > c)
-    {
-        return this->operator[](r)[c];
-    }
-    else
-    {
-        return this->operator[](c)[r];
-    }
-}
-
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
index bbe1e2148..1bb8ab1ca 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
@@ -60,40 +60,40 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
         {
             for (label k = i; k < Um; k++)
             {
-                scale += mag(U_[k][i]);
+                scale += mag(U_(k, i));
             }
 
             if (scale != 0)
             {
                 for (label k = i; k < Um; k++)
                 {
-                    U_[k][i] /= scale;
-                    s += U_[k][i]*U_[k][i];
+                    U_(k, i) /= scale;
+                    s += U_(k, i)*U_(k, i);
                 }
 
-                scalar f = U_[i][i];
+                scalar f = U_(i, i);
                 g = -sign(Foam::sqrt(s), f);
                 scalar h = f*g - s;
-                U_[i][i] = f - g;
+                U_(i, i) = f - g;
 
                 for (label j = l-1; j < Un; j++)
                 {
                     s = 0;
                     for (label k = i; k < Um; k++)
                     {
-                        s += U_[k][i]*U_[k][j];
+                        s += U_(k, i)*U_(k, j);
                     }
 
                     f = s/h;
                     for (label k = i; k < A.m(); k++)
                     {
-                        U_[k][j] += f*U_[k][i];
+                        U_(k, j) += f*U_(k, i);
                     }
                 }
 
                 for (label k = i; k < Um; k++)
                 {
-                    U_[k][i] *= scale;
+                    U_(k, i) *= scale;
                 }
             }
         }
@@ -106,15 +106,15 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
         {
             for (label k = l-1; k < Un; k++)
             {
-                scale += mag(U_[i][k]);
+                scale += mag(U_(i, k));
             }
 
             if (scale != 0)
             {
                 for (label k=l-1; k < Un; k++)
                 {
-                    U_[i][k] /= scale;
-                    s += U_[i][k]*U_[i][k];
+                    U_(i, k) /= scale;
+                    s += U_(i, k)*U_(i, k);
                 }
 
                 scalar f = U_[i][l-1];
@@ -124,7 +124,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 
                 for (label k = l-1; k < Un; k++)
                 {
-                    rv1[k] = U_[i][k]/h;
+                    rv1[k] = U_(i, k)/h;
                 }
 
                 for (label j = l-1; j < Um; j++)
@@ -132,17 +132,17 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                     s = 0;
                     for (label k = l-1; k < Un; k++)
                     {
-                        s += U_[j][k]*U_[i][k];
+                        s += U_(j, k)*U_(i, k);
                     }
 
                     for (label k = l-1; k < Un; k++)
                     {
-                        U_[j][k] += s*rv1[k];
+                        U_(j, k) += s*rv1[k];
                     }
                 }
                 for (label k = l-1; k < Un; k++)
                 {
-                    U_[i][k] *= scale;
+                    U_(i, k) *= scale;
                 }
             }
         }
@@ -158,7 +158,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
             {
                 for (label j = l; j < Un; j++)
                 {
-                    V_[j][i] = (U_[i][j]/U_[i][l])/g;
+                    V_(j, i) = (U_(i, j)/U_(i, l))/g;
                 }
 
                 for (label j=l; j < Un; j++)
@@ -166,23 +166,23 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                     s = 0;
                     for (label k = l; k < Un; k++)
                     {
-                        s += U_[i][k]*V_[k][j];
+                        s += U_(i, k)*V_(k, j);
                     }
 
                     for (label k = l; k < Un; k++)
                     {
-                        V_[k][j] += s*V_[k][i];
+                        V_(k, j) += s*V_(k, i);
                     }
                 }
             }
 
             for (label j = l; j < Un;j++)
             {
-                V_[i][j] = V_[j][i] = 0.0;
+                V_(i, j) = V_(j, i) = 0.0;
             }
         }
 
-        V_[i][i] = 1;
+        V_(i, i) = 1;
         g = rv1[i];
         l = i;
     }
@@ -194,7 +194,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 
         for (label j = l; j < Un; j++)
         {
-            U_[i][j] = 0.0;
+            U_(i, j) = 0.0;
         }
 
         if (g != 0)
@@ -206,31 +206,31 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                 s = 0;
                 for (label k = l; k < Um; k++)
                 {
-                    s += U_[k][i]*U_[k][j];
+                    s += U_(k, i)*U_(k, j);
                 }
 
-                scalar f = (s/U_[i][i])*g;
+                scalar f = (s/U_(i, i))*g;
 
                 for (label k = i; k < Um; k++)
                 {
-                    U_[k][j] += f*U_[k][i];
+                    U_(k, j) += f*U_(k, i);
                 }
             }
 
             for (label j = i; j < Um; j++)
             {
-                U_[j][i] *= g;
+                U_(j, i) *= g;
             }
         }
         else
         {
             for (label j = i; j < Um; j++)
             {
-                U_[j][i] = 0.0;
+                U_(j, i) = 0.0;
             }
         }
 
-        ++U_[i][i];
+        ++U_(i, i);
     }
 
     for (label k = Un-1; k >= 0; k--)
@@ -272,9 +272,9 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                     for (label j = 0; j < Um; j++)
                     {
                         scalar y = U_[j][mn];
-                        scalar z = U_[j][i];
+                        scalar z = U_(j, i);
                         U_[j][mn] = y*c + z*s;
-                        U_[j][i] = z*c - y*s;
+                        U_(j, i) = z*c - y*s;
                     }
                 }
             }
@@ -287,7 +287,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                 {
                     S_[k] = -z;
 
-                    for (label j = 0; j < Un; j++) V_[j][k] = -V_[j][k];
+                    for (label j = 0; j < Un; j++) V_(j, k) = -V_(j, k);
                 }
                 break;
             }
@@ -383,7 +383,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
     {
         for (label j = 0; j < A.n(); j++)
         {
-            diff = mag(A[i][j] - SVDA[i][j]);
+            diff = mag(A(i, j) - SVDA(i, j));
             if (diff > maxDiff) maxDiff = diff;
         }
     }
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index 46e972f72..0317fc27a 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
@@ -84,7 +84,7 @@ void Foam::LUDecompose
             scalar sum = matrixi[j];
             for (label k=0; k<i; k++)
             {
-                sum -= matrixi[k]*matrix[k][j];
+                sum -= matrixi[k]*matrix(k, j);
             }
             matrixi[j] = sum;
         }
@@ -99,7 +99,7 @@ void Foam::LUDecompose
 
             for (label k=0; k<j; k++)
             {
-                sum -= matrixi[k]*matrix[k][j];
+                sum -= matrixi[k]*matrix(k, j);
             }
 
             matrixi[j] = sum;
@@ -138,7 +138,7 @@ void Foam::LUDecompose
 
             for (label i=j+1; i<m; i++)
             {
-                matrix[i][j] *= rDiag;
+                matrix(i, j) *= rDiag;
             }
         }
     }
@@ -155,7 +155,7 @@ void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix)
     {
         for (label k = j + 1; k < size; k++)
         {
-            matrix[j][k] = 0.0;
+            matrix(j, k) = 0.0;
         }
     }
 
@@ -169,18 +169,18 @@ void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix)
 
             for (label i = 0; i < k; i++)
             {
-                s += matrix[i][k]*matrix[i][j];
+                s += matrix(i, k)*matrix(i, j);
             }
 
-            s = (matrix[j][k] - s)/matrix[k][k];
+            s = (matrix(j, k) - s)/matrix(k, k);
 
-            matrix[k][j] = s;
-            matrix[j][k] = s;
+            matrix(k, j) = s;
+            matrix(j, k) = s;
 
             d += sqr(s);
         }
 
-        d = matrix[j][j] - d;
+        d = matrix(j, j) - d;
 
         if (d < 0.0)
         {
@@ -190,7 +190,7 @@ void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix)
                 << abort(FatalError);
         }
 
-        matrix[j][j] = sqrt(d);
+        matrix(j, j) = sqrt(d);
     }
 }
 
@@ -232,9 +232,9 @@ void Foam::multiply
                 scalar ab = 0;
                 for (label j = 0; j < A.n(); j++)
                 {
-                    ab += A[i][j]*B[j][l];
+                    ab += A(i, j)*B(j, l);
                 }
-                ans[i][g] += C[l][g] * ab;
+                ans(i, g) += C(l, g) * ab;
             }
         }
     }
@@ -273,7 +273,7 @@ void Foam::multiply
         {
             for (label l = 0; l < C.m(); l++)
             {
-                ans[i][g] += C[l][g] * A[i][l]*B[l];
+                ans(i, g) += C(l, g) * A(i, l)*B[l];
             }
         }
     }
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
index 077fa60ce..c6acc739e 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
@@ -47,7 +47,7 @@ void Foam::solve
         // Swap entries around to find a good pivot
         for (label j=i+1; j<m; j++)
         {
-            if (mag(tmpMatrix[j][i]) > largestCoeff)
+            if (mag(tmpMatrix(j, i)) > largestCoeff)
             {
                 iMax = j;
                 largestCoeff = mag(tmpMatrix[iMax][i]);
@@ -58,13 +58,13 @@ void Foam::solve
         {
             for (label k=i; k<m; k++)
             {
-                Swap(tmpMatrix[i][k], tmpMatrix[iMax][k]);
+                Swap(tmpMatrix(i, k), tmpMatrix[iMax][k]);
             }
             Swap(sourceSol[i], sourceSol[iMax]);
         }
 
         // Check that the system of equations isn't singular
-        if (mag(tmpMatrix[i][i]) < 1e-20)
+        if (mag(tmpMatrix(i, i)) < 1e-20)
         {
             FatalErrorInFunction
                 << "Singular Matrix"
@@ -74,12 +74,12 @@ void Foam::solve
         // Reduce to upper triangular form
         for (label j=i+1; j<m; j++)
         {
-            sourceSol[j] -= sourceSol[i]*(tmpMatrix[j][i]/tmpMatrix[i][i]);
+            sourceSol[j] -= sourceSol[i]*(tmpMatrix(j, i)/tmpMatrix(i, i));
 
             for (label k=m-1; k>=i; k--)
             {
-                tmpMatrix[j][k] -=
-                    tmpMatrix[i][k]*tmpMatrix[j][i]/tmpMatrix[i][i];
+                tmpMatrix(j, k) -=
+                    tmpMatrix(i, k)*tmpMatrix(j, i)/tmpMatrix(i, i);
             }
         }
     }
@@ -91,10 +91,10 @@ void Foam::solve
 
         for (label k=j+1; k<m; k++)
         {
-            ntempvec += tmpMatrix[j][k]*sourceSol[k];
+            ntempvec += tmpMatrix(j, k)*sourceSol[k];
         }
 
-        sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix[j][j];
+        sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix(j, j);
     }
 }
 
@@ -257,7 +257,7 @@ void Foam::multiply
         {
             for (label l = 0; l < B.m(); l++)
             {
-                ans[i][j] += A[i][l]*B[l][j];
+                ans(i, j) += A(i, l)*B(l, j);
             }
         }
     }
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C
index 8820a78e6..b1bf85aa8 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C
@@ -122,8 +122,8 @@ void Foam::CentredFitSnGradData<Polynomial>::calcFit
     // Additional weighting for constant and linear terms
     for (label i = 0; i < B.m(); i++)
     {
-        B[i][0] *= wts[0];
-        B[i][1] *= wts[0];
+        B(i, 0) *= wts[0];
+        B(i, 1) *= wts[0];
     }
 
     // Set the fit
@@ -137,14 +137,14 @@ void Foam::CentredFitSnGradData<Polynomial>::calcFit
 
         for (label i=0; i<stencilSize; i++)
         {
-            coeffsi[i] = wts[1]*wts[i]*svd.VSinvUt()[1][i]/scale;
+            coeffsi[i] = wts[1]*wts[i]*svd.VSinvUt()(1, i)/scale;
         }
 
         goodFit =
         (
-            mag(wts[0]*wts[0]*svd.VSinvUt()[0][0] - wLin)
+            mag(wts[0]*wts[0]*svd.VSinvUt()(0, 0) - wLin)
           < this->linearLimitFactor()*wLin)
-         && (mag(wts[0]*wts[1]*svd.VSinvUt()[0][1] - (1 - wLin)
+         && (mag(wts[0]*wts[1]*svd.VSinvUt()(0, 1) - (1 - wLin)
         ) < this->linearLimitFactor()*(1 - wLin))
          && coeffsi[0] < 0 && coeffsi[1] > 0
          && mag(coeffsi[0] + deltaCoeff) < 0.5*deltaCoeff
@@ -163,10 +163,10 @@ void Foam::CentredFitSnGradData<Polynomial>::calcFit
                 << "    deltaCoeff " << deltaCoeff << nl
                 << "    sing vals " << svd.S() << nl
                 << "Components of goodFit:\n"
-                << "    wts[0]*wts[0]*svd.VSinvUt()[0][0] = "
-                << wts[0]*wts[0]*svd.VSinvUt()[0][0] << nl
-                << "    wts[0]*wts[1]*svd.VSinvUt()[0][1] = "
-                << wts[0]*wts[1]*svd.VSinvUt()[0][1]
+                << "    wts[0]*wts[0]*svd.VSinvUt()(0, 0) = "
+                << wts[0]*wts[0]*svd.VSinvUt()(0, 0) << nl
+                << "    wts[0]*wts[1]*svd.VSinvUt()(0, 1) = "
+                << wts[0]*wts[1]*svd.VSinvUt()(0, 1)
                 << " dim = " << this->dim() << endl;
 
             wts[0] *= 10;
@@ -174,14 +174,14 @@ void Foam::CentredFitSnGradData<Polynomial>::calcFit
 
             for (label j = 0; j < B.n(); j++)
             {
-                B[0][j] *= 10;
-                B[1][j] *= 10;
+                B(0, j) *= 10;
+                B(1, j) *= 10;
             }
 
             for (label i = 0; i < B.m(); i++)
             {
-                B[i][0] *= 10;
-                B[i][1] *= 10;
+                B(i, 0) *= 10;
+                B(i, 1) *= 10;
             }
         }
     }
diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
index bc86903cb..eb6da09f1 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
@@ -196,8 +196,8 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
     // Additional weighting for constant and linear terms
     for (label i = 0; i < B.m(); i++)
     {
-        B[i][0] *= wts[0];
-        B[i][1] *= wts[0];
+        B(i, 0) *= wts[0];
+        B(i, 1) *= wts[0];
     }
 
     // Set the fit
@@ -214,7 +214,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
 
         for (label i=0; i<stencilSize; i++)
         {
-            coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
+            coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()(0, i);
             if (mag(coeffsi[i]) > maxCoeff)
             {
                 maxCoeff = mag(coeffsi[i]);
@@ -269,14 +269,14 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
 
             for (label j = 0; j < B.n(); j++)
             {
-                B[0][j] *= 10;
-                B[1][j] *= 10;
+                B(0, j) *= 10;
+                B(1, j) *= 10;
             }
 
             for (label i = 0; i < B.m(); i++)
             {
-                B[i][0] *= 10;
-                B[i][1] *= 10;
+                B(i, 0) *= 10;
+                B(i, 1) *= 10;
             }
         }
     }
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 7c2d9beac..1b424d1cc 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -341,7 +341,7 @@ void Foam::chemistryModel<CompType, ThermoType>::jacobian
     {
         for (label j=0; j<nEqns(); j++)
         {
-            dfdc[i][j] = 0.0;
+            dfdc(i, j) = 0.0;
         }
     }
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
index 15a256f66..2b52c0a7c 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -151,7 +151,7 @@ void Foam::EulerImplicit<ChemistryModel>::solve
         scalar d = 0;
         for (label j=0; j<nSpecie; j++)
         {
-            d -= RR[i][j]*c[j];
+            d -= RR(i, j)*c[j];
         }
 
         if (d < -SMALL)
@@ -172,7 +172,7 @@ void Foam::EulerImplicit<ChemistryModel>::solve
     // Add the diagonal and source contributions from the time-derivative
     for (label i=0; i<nSpecie; i++)
     {
-        RR[i][i] += 1.0/deltaT;
+        RR(i, i) += 1.0/deltaT;
         RR.source()[i] = c[i]/deltaT;
     }
 
diff --git a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
index 71e34ab46..2c8a7e2ba 100644
--- a/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
+++ b/src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
@@ -209,12 +209,12 @@ void Foam::radiation::viewFactor::initialise()
                 scalar sumF = 0.0;
                 for (label j=0; j<totalNCoarseFaces_; j++)
                 {
-                    sumF += Fmatrix_()[i][j];
+                    sumF += Fmatrix_()(i, j);
                 }
                 scalar delta = sumF - 1.0;
                 for (label j=0; j<totalNCoarseFaces_; j++)
                 {
-                    Fmatrix_()[i][j] *= (1.0 - delta/(sumF + 0.001));
+                    Fmatrix_()(i, j) *= (1.0 - delta/(sumF + 0.001));
                 }
             }
         }
@@ -541,13 +541,13 @@ void Foam::radiation::viewFactor::calculate()
 
                     if (i==j)
                     {
-                        C[i][j] = invEj - (invEj - 1.0)*Fmatrix_()[i][j];
-                        q[i] += (Fmatrix_()[i][j] - 1.0)*sigmaT4 - QrExt[j];
+                        C(i, j) = invEj - (invEj - 1.0)*Fmatrix_()(i, j);
+                        q[i] += (Fmatrix_()(i, j) - 1.0)*sigmaT4 - QrExt[j];
                     }
                     else
                     {
-                        C[i][j] = (1.0 - invEj)*Fmatrix_()[i][j];
-                        q[i] += Fmatrix_()[i][j]*sigmaT4;
+                        C(i, j) = (1.0 - invEj)*Fmatrix_()(i, j);
+                        q[i] += Fmatrix_()(i, j)*sigmaT4;
                     }
 
                 }
@@ -570,11 +570,11 @@ void Foam::radiation::viewFactor::calculate()
                         scalar invEj = 1.0/E[j];
                         if (i==j)
                         {
-                            CLU_()[i][j] = invEj-(invEj-1.0)*Fmatrix_()[i][j];
+                            CLU_()(i, j) = invEj-(invEj-1.0)*Fmatrix_()(i, j);
                         }
                         else
                         {
-                            CLU_()[i][j] = (1.0 - invEj)*Fmatrix_()[i][j];
+                            CLU_()(i, j) = (1.0 - invEj)*Fmatrix_()(i, j);
                         }
                     }
                 }
@@ -598,11 +598,11 @@ void Foam::radiation::viewFactor::calculate()
 
                     if (i==j)
                     {
-                        q[i] += (Fmatrix_()[i][j] - 1.0)*sigmaT4  - QrExt[j];
+                        q[i] += (Fmatrix_()(i, j) - 1.0)*sigmaT4  - QrExt[j];
                     }
                     else
                     {
-                        q[i] += Fmatrix_()[i][j]*sigmaT4;
+                        q[i] += Fmatrix_()(i, j)*sigmaT4;
                     }
                 }
             }
diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
index dee571fe1..54992d7ef 100644
--- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C
@@ -363,7 +363,7 @@ jacobian
     {
         for (label j=0; j<nEqns(); j++)
         {
-            dfdc[i][j] = 0.0;
+            dfdc(i, j) = 0.0;
         }
     }