diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py
index 15dcb09c..3ee037df 100644
--- a/SimPEG/LogicallyOrthogonalMesh.py
+++ b/SimPEG/LogicallyOrthogonalMesh.py
@@ -87,7 +87,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
         def fget(self):
             if(self._vol is None):
                 if self.dim == 2:
-                    A, B, C, D = indexCube('ABCD', np.array([self.nNx, self.nNy]))
+                    A, B, C, D = indexCube('ABCD', self.n+1)
                     normal, area, length = faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
                     self._vol = area
                 elif self.dim == 3:
@@ -97,7 +97,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
                     # T3 = [B D E G]; % mid
                     # T4 = [B C D G]; % cutted edge
                     # T5 = [D E G H]; % cutted edge
-                    A, B, C, D, E, F, G, H = indexCube('ABCDEFGH', np.array([self.nNx, self.nNy, self.nNz]))
+                    A, B, C, D, E, F, G, H = indexCube('ABCDEFGH', self.n+1)
 
                     v1 = volTetra(self.gridN, A, B, D, E)  # cutted edge
                     v2 = volTetra(self.gridN, B, E, F, G)  # cutted edge
@@ -111,13 +111,46 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
     _vol = None
     vol = property(**vol())
 
+    def area():
+        doc = "Face areas."
+
+        def fget(self):
+            if(self._area is None):
+                # Compute areas of cell faces
+                if(self.dim == 2):
+                    xy = self.gridN
+                    length = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
+
+                    A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
+                    area1 = length(xy[B, :] - xy[A, :])
+                    A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
+                    area2 = length(xy[D, :] - xy[A, :])
+                    self._area = np.r_[mkvc(area1), mkvc(area2)]
+                elif(self.dim == 3):
+
+                    A, E, F, B = indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz]))
+                    normal, area1, length = faceInfo(self.gridN, A, E, F, B)
+
+                    A, D, H, E = indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz]))
+                    normal, area2, length = faceInfo(self.gridN, A, D, H, E)
+
+                    A, B, C, D = indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz]))
+                    normal, area3, length = faceInfo(self.gridN, A, B, C, D)
+
+                    self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)]
+            return self._area
+        return locals()
+    _area = None
+    area = property(**area())
+
 
 if __name__ == '__main__':
     nc = 5
     h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
+    nc = 7
     h2 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
     h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
-    dee3 = False
+    dee3 = True
     if dee3:
         X, Y, Z = ndgrid(h1, h2, h3, vector=False)
         M = LogicallyOrthogonalMesh([X, Y, Z])
@@ -127,4 +160,5 @@ if __name__ == '__main__':
 
     # print M.r(M.gridCC, format='M')
     # print M.gridN[:, 0]
-    print np.sum(M.vol)
+    print M.nE
+    print M.area
diff --git a/SimPEG/utils.py b/SimPEG/utils.py
index eff11748..9847751f 100644
--- a/SimPEG/utils.py
+++ b/SimPEG/utils.py
@@ -124,14 +124,16 @@ def volTetra(xyz, A, B, C, D):
     return V/6
 
 
-def indexCube(nodes, nN):
+def indexCube(nodes, gridSize, n=None):
     """
     Returns the index of nodes on the mesh.
 
 
     Input:
-       nodes  - string of which nodes to return. e.g. 'ABCD'
-       nN     - size of the nodal grid
+       nodes     - string of which nodes to return. e.g. 'ABCD'
+       gridSize  - size of the nodal grid
+       n         - number of nodes each i,j,k direction: [ni,nj,nk]
+
 
     Output:
        index  - index in the order asked e.g. 'ABCD' --> (A,B,C,D)
@@ -171,20 +173,21 @@ def indexCube(nodes, nN):
     """
 
     assert type(nodes) == str, "Nodes must be a str variable: e.g. 'ABCD'"
-    assert type(nN) == np.ndarray, "Number of nodes must be an ndarray"
+    assert type(gridSize) == np.ndarray, "Number of nodes must be an ndarray"
     nodes = nodes.upper()
     # Make sure that we choose from the possible nodes.
-    possibleNodes = 'ABCD' if nN.size == 2 else 'ABCDEFGH'
+    possibleNodes = 'ABCD' if gridSize.size == 2 else 'ABCDEFGH'
     for node in nodes:
         assert node in possibleNodes, "Nodes must be chosen from: '%s'" % possibleNodes
-    dim = nN.size
-    nC = nN - 1
+    dim = gridSize.size
+    if n is None:
+        n = gridSize - 1
 
     if dim == 2:
-        ij = ndgrid(np.arange(nC[0]), np.arange(nC[1]))
+        ij = ndgrid(np.arange(n[0]), np.arange(n[1]))
         i, j = ij[:, 0], ij[:, 1]
     elif dim == 3:
-        ijk = ndgrid(np.arange(nC[0]), np.arange(nC[1]), np.arange(nC[2]))
+        ijk = ndgrid(np.arange(n[0]), np.arange(n[1]), np.arange(n[2]))
         i, j, k = ijk[:, 0], ijk[:, 1], ijk[:, 2]
     else:
         raise Exception('Only 2 and 3 dimensions supported.')
@@ -195,9 +198,9 @@ def indexCube(nodes, nN):
     for node in nodes:
         shift = nodeMap[node]
         if dim == 2:
-            out += (sub2ind(nN, np.c_[i+shift[0], j+shift[1]]).flatten(), )
+            out += (sub2ind(gridSize, np.c_[i+shift[0], j+shift[1]]).flatten(), )
         elif dim == 3:
-            out += (sub2ind(nN, np.c_[i+shift[0], j+shift[1], k+shift[2]]).flatten(), )
+            out += (sub2ind(gridSize, np.c_[i+shift[0], j+shift[1], k+shift[2]]).flatten(), )
 
     return out