diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py
index 5bb5482b..15dcb09c 100644
--- a/SimPEG/LogicallyOrthogonalMesh.py
+++ b/SimPEG/LogicallyOrthogonalMesh.py
@@ -88,7 +88,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators):  # , LOMGrid
             if(self._vol is None):
                 if self.dim == 2:
                     A, B, C, D = indexCube('ABCD', np.array([self.nNx, self.nNy]))
-                    normal, area, length = faceInfo(np.c_[self.gridN, np.zeros((self.nC, 1))], A, B, C, D)
+                    normal, area, length = faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
                     self._vol = area
                 elif self.dim == 3:
                     # Each polyhedron can be decomposed into 5 tetrahedrons
@@ -117,8 +117,14 @@ if __name__ == '__main__':
     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.r(M.gridCC, format='M')
-    print M.gridN[:, 0]
+    dee3 = False
+    if dee3:
+        X, Y, Z = ndgrid(h1, h2, h3, vector=False)
+        M = LogicallyOrthogonalMesh([X, Y, Z])
+    else:
+        X, Y = ndgrid(h1, h2, vector=False)
+        M = LogicallyOrthogonalMesh([X, Y])
+
+    # print M.r(M.gridCC, format='M')
+    # print M.gridN[:, 0]
+    print np.sum(M.vol)
diff --git a/SimPEG/utils.py b/SimPEG/utils.py
index 37441e8e..eff11748 100644
--- a/SimPEG/utils.py
+++ b/SimPEG/utils.py
@@ -247,9 +247,9 @@ def faceInfo(xyz, A, B, C, D, average=True):
     DA = xyz[A, :] - xyz[D, :]
 
     def cross(X, Y):
-        return np.c_[X[1, :]*Y[2, :] - X[2, :]*Y[1, :],
-                     X[2, :]*Y[0, :] - X[0, :]*Y[2, :],
-                     X[0, :]*Y[1, :] - X[1, :]*Y[0, :]]
+        return np.c_[X[:, 1]*Y[:, 2] - X[:, 2]*Y[:, 1],
+                     X[:, 2]*Y[:, 0] - X[:, 0]*Y[:, 2],
+                     X[:, 0]*Y[:, 1] - X[:, 1]*Y[:, 0]]
 
     nA = cross(AB, DA)
     nB = cross(BC, AB)
@@ -257,7 +257,7 @@ def faceInfo(xyz, A, B, C, D, average=True):
     nD = cross(DA, CD)
 
     length = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
-    normalize = lambda x: x/np.kron(np.ones((1, x.shape[1]), length(x), 1))
+    normalize = lambda x: x/np.kron(np.ones((1, x.shape[1])), mkvc(length(N), 2))
     if average:
         # average the normals at each vertex.
         N = (nA + nB + nC + nD)/4  # this is intrinsically weighted by area