From 358411792dac59a44556b240e2daf0942fbdbc68 Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@visiblegeology.com>
Date: Wed, 13 Nov 2013 19:37:14 -0800
Subject: [PATCH] Cleaned up averaging. Put it in an AveExtrap function.

---
 SimPEG/mesh/DiffOperators.py | 25 +++++++------------------
 SimPEG/utils/__init__.py     |  2 +-
 SimPEG/utils/sputils.py      |  7 ++++++-
 3 files changed, 14 insertions(+), 20 deletions(-)

diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py
index 2c94a0cf..7176b8fd 100644
--- a/SimPEG/mesh/DiffOperators.py
+++ b/SimPEG/mesh/DiffOperators.py
@@ -1,6 +1,6 @@
 import numpy as np
 from scipy import sparse as sp
-from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av
+from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap
 
 
 def checkBC(bc):
@@ -423,25 +423,14 @@ class DiffOperators(object):
             if(self._aveCC2F is None):
                 n = self.n
                 if(self.dim == 1):
-                    Av = av(n[0]).T
-                    Av = sdiag(1/Av.sum(axis=1))*Av
-                    self._aveCC2F = Av
+                    self._aveCC2F = avExtrap(n[0])
                 elif(self.dim == 2):
-                    Av1 = av(n[0]).T
-                    Av1 = sdiag(1/Av1.sum(axis=1))*Av1
-                    Av2 = av(n[1]).T
-                    Av2 = sdiag(1/Av2.sum(axis=1))*Av2
-                    Av = sp.vstack((sp.kron(speye(n[1]), Av1), sp.kron(Av2, speye(n[0]))), format="csr")
-                    self._aveCC2F = Av
+                    self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
+                                               sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
                 elif(self.dim == 3):
-                    Av1 = av(n[0]).T
-                    Av1 = sdiag(1/Av1.sum(axis=1))*Av1
-                    Av2 = av(n[1]).T
-                    Av2 = sdiag(1/Av2.sum(axis=1))*Av2
-                    Av3 = av(n[2]).T
-                    Av3 = sdiag(1/Av3.sum(axis=1))*Av3
-                    Av = sp.vstack((kron3(speye(n[2]), speye(n[1]), Av1), kron3(speye(n[2]), Av2, speye(n[0])), kron3(Av3, speye(n[1]), speye(n[0]))), format="csr")
-                    self._aveCC2F = Av
+                    self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
+                                               kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
+                                               kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
             return self._aveCC2F
         return locals()
     _aveCC2F = None
diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py
index d2170721..75e7c360 100644
--- a/SimPEG/utils/__init__.py
+++ b/SimPEG/utils/__init__.py
@@ -4,7 +4,7 @@ import lomutils
 import interputils
 import ModelBuilder
 from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
-from sputils import spzeros, kron3, speye, sdiag, ddx, av
+from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap
 from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
 from interputils import interpmat
 from ipythonUtils import easyAnimate as animate
diff --git a/SimPEG/utils/sputils.py b/SimPEG/utils/sputils.py
index e066322b..06eef80d 100644
--- a/SimPEG/utils/sputils.py
+++ b/SimPEG/utils/sputils.py
@@ -29,5 +29,10 @@ def ddx(n):
 
 
 def av(n):
-    """Define 1D averaging operator from cell-centres to nodes."""
+    """Define 1D averaging operator from nodes to cell-centers."""
     return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
+
+def avExtrap(n):
+    """Define 1D averaging operator from cell-centers to nodes."""
+    Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
+    return Av