From 100380aefe509adbfe0b76f43d09814b07feb1c3 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Mon, 31 Oct 2016 10:27:03 +0100
Subject: [PATCH] j_at components are now correct (still open: conditional for
 gradient=zero)

---
 finitedifferences.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/finitedifferences.py b/finitedifferences.py
index 5e6b0ac56..4038f34b0 100644
--- a/finitedifferences.py
+++ b/finitedifferences.py
@@ -61,7 +61,7 @@ def discretizeStaggered(term, symbolsToFieldDict, coordinate, coordinateOffset,
     """
     Expects term that contains given symbols and gradient components of these symbols and replaces them
     by field accesses. Gradients in coordinate direction  are replaced by staggered version at cell boundary.
-    Symbols themselves are replaced by interpolated version at boundary.
+    Symbols themselves and gradients in other directions are replaced by interpolated version at cell face.
 
     :param term: input term where symbols and gradients are replaced
     :param symbolsToFieldDict: mapping of symbols to Field
@@ -96,16 +96,16 @@ def discretizeStaggered(term, symbolsToFieldDict, coordinate, coordinateOffset,
         gradient = grad(symbols)[coordinate]
         substitutions.update({s: (field[offset](i) + field(i)) / 2 for i, s in enumerate(symbols)})
         substitutions.update({g: (field[offset](i) - field(i)) / dx * coordinateOffset for i, g in enumerate(gradient)})
-        #for d in range(dim):
-        #    if d == coordinate:
-        #        continue
-        #    up, down = __upDownOffsets(d, dim)
-        #    for i, s in enumerate(symbols):
-        #        centerGrad = (field[up](i) - field[down](i)) / (2 * dx)
-        #        neighborGrad = (field[up+offset](i) - field[down+offset](i)) / (2 * dx)
-        #        substitutions[grad(s)[d]] = (centerGrad + neighborGrad) / 2
-
-    return term.subs(substitutions)#, substitutions
+        for d in range(dim):
+            if d == coordinate:
+                continue
+            up, down = __upDownOffsets(d, dim)
+            for i, s in enumerate(symbols):
+                centerGrad = (field[up](i) - field[down](i)) / (2 * dx)
+                neighborGrad = (field[up+offset](i) - field[down+offset](i)) / (2 * dx)
+                substitutions[grad(s)[d]] = (centerGrad + neighborGrad) / 2
+
+    return term.subs(substitutions)
 
 
 def discretizeDivergence(vectorTerm, symbolsToFieldDict, dx):
-- 
GitLab