From 98c0f0e1a77226db6c79ccce64a0dea20fe00094 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Wed, 22 Feb 2017 17:02:09 +0100
Subject: [PATCH] [staggeredGrid][localJacobian] Fix erroneous index

* faceIndices start at numCellCenterEquations ( > 0)
* Accessing a Dune::FieldVector<int,1> with an index > 0 does not yield any
  error, warning, segfault etc.
---
 dumux/implicit/staggered/localjacobian.hh | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index 027747031f..28e12e50a9 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -330,7 +330,7 @@ private:
                 partialDeriv /= eps;
 
                 // update the global jacobian matrix with the current partial derivatives
-                this->updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], ccGlobalI_, globalJ, pvIdx, partialDeriv);
+                this->updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], ccGlobalI_, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
 
                 // restore the original faceVars
                 curFaceVars = origFaceVars;
@@ -436,7 +436,7 @@ private:
                     partialDeriv /= eps;
 
                     // update the global jacobian matrix with the current partial derivatives
-                    this->updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx, partialDeriv);
+                    this->updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
 
                     // restore the original faceVars
                     curFaceVars = origFaceVars;
@@ -521,6 +521,10 @@ private:
             // the residual of equation 'eqIdx' at dof 'i'
             // depending on the primary variable 'pvIdx' at dof
             // 'col'.
+
+            assert(pvIdx >= 0);
+            assert(eqIdx < matrix[globalI][globalJ].size());
+            assert(pvIdx < matrix[globalI][globalJ][eqIdx].size());
             matrix[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
             Valgrind::CheckDefined(matrix[globalI][globalJ][eqIdx][pvIdx]);
         }
-- 
GitLab