From dc9a92c83a5a36c44a1be1bb717c66bfded34ea5 Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@visiblegeology.com>
Date: Wed, 24 Jul 2013 12:25:31 -0700
Subject: [PATCH] Cell centered grid working in LOM

---
 SimPEG/DiffOperators.py           | 42 +++++++++++++++++++++++++++++++
 SimPEG/LogicallyOrthogonalMesh.py | 22 +++++++++++++---
 2 files changed, 60 insertions(+), 4 deletions(-)

diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py
index c5daa96d..5d8356bd 100644
--- a/SimPEG/DiffOperators.py
+++ b/SimPEG/DiffOperators.py
@@ -245,6 +245,48 @@ class DiffOperators(object):
     _edgeAve = None
     edgeAve = property(**edgeAve())
 
+    def nodalAve():
+        doc = "Construct the averaging operator on cell nodes to cell centers."
+
+        def fget(self):
+            if(self._nodalAve is None):
+                # The number of cell centers in each direction
+                n = self.n
+                if(self.dim == 1):
+                    self._nodalAve = av(n[0])
+                elif(self.dim == 2):
+                    self._nodalAve = sp.hstack((sp.kron(av(n[1]), av(n[0])),
+                                                sp.kron(av(n[1]), av(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._nodalAve = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])),
+                                                kron3(av(n[2]), av(n[1]), av(n[0])),
+                                                kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
+            return self._nodalAve
+        return locals()
+    _nodalAve = None
+    nodalAve = property(**nodalAve())
+
+    def nodalVectorAve():
+        doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension seperate."
+
+        def fget(self):
+            if(self._nodalVectorAve is None):
+                # The number of cell centers in each direction
+                n = self.n
+                if(self.dim == 1):
+                    self._nodalVectorAve = av(n[0])
+                elif(self.dim == 2):
+                    self._nodalVectorAve = sp.block_diag((sp.kron(av(n[1]), av(n[0])),
+                                                          sp.kron(av(n[1]), av(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._nodalVectorAve = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])),
+                                                          kron3(av(n[2]), av(n[1]), av(n[0])),
+                                                          kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
+            return self._nodalVectorAve
+        return locals()
+    _nodalVectorAve = None
+    nodalVectorAve = property(**nodalVectorAve())
+
     def getEdgeMass(self, materialProp=None):
         """mass matix for products of edge functions w'*M(materialProp)*e"""
         if(materialProp is None):
diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py
index b530fcb1..b49d595c 100644
--- a/SimPEG/LogicallyOrthogonalMesh.py
+++ b/SimPEG/LogicallyOrthogonalMesh.py
@@ -1,7 +1,7 @@
 import numpy as np
 from BaseMesh import BaseMesh
 from DiffOperators import DiffOperators
-from utils import mkvc
+from utils import mkvc, ndgrid
 
 
 class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
@@ -17,7 +17,9 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
             assert type(nodes_i) == np.ndarray, ("nodes[%i] is not a numpy array." % i)
             assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i)
 
-        super(LogicallyOrthogonalMesh, self).__init__(np.array(nodes[0].shape), x0)
+        assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
+
+        super(LogicallyOrthogonalMesh, self).__init__(np.array(nodes[0].shape)-1, x0)
 
         assert len(nodes[0].shape) == len(self.x0), "Dimension mismatch. x0 != len(h)"
 
@@ -31,7 +33,8 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
 
         def fget(self):
             if self._gridCC is None:
-                self._gridCC = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorCCz] if not x is None])
+                ccV = (self.nodalVectorAve*mkvc(self.gridN))
+                self._gridCC = ccV.reshape((-1, self.dim), order='F')
             return self._gridCC
         return locals()
     _gridCC = None  # Store grid by default
@@ -42,7 +45,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
 
         def fget(self):
             if self._gridN is None:
-                self._gridN = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorNz] if not x is None])
+                raise Exception("Someone deleted this. I blame you.")
             return self._gridN
         return locals()
     _gridN = None  # Store grid by default
@@ -113,3 +116,14 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
         return locals()
     _gridEz = None  # Store grid by default
     gridEz = property(**gridEz())
+
+if __name__ == '__main__':
+    nc = 5
+    h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
+    h2 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
+    h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
+    h = [h1, h2, h3]
+    X, Y, Z = ndgrid(h1, h2, h3, vector=False)
+    M = LogicallyOrthogonalMesh([X, Y, Z])
+    print M.gridCC[:,0]
+    print M.gridN[:,0]