diff --git a/SimPEG/getDiffop.py b/SimPEG/getDiffop.py
new file mode 100644
index 00000000..6debc348
--- /dev/null
+++ b/SimPEG/getDiffop.py
@@ -0,0 +1,198 @@
+import numpy as np
+from scipy import sparse
+from utils import mkvc
+from sputils import ddx, sdiag, speye, kron3, spzeros, av
+
+def getvol(h):   
+    """Construct cell volumes of the 3D model as 1d array."""        
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # Compute cell volumes                        
+    v12 = h1.T*h2       
+    V = mkvc(v12.reshape(-1,1)*h3)
+    
+    return V
+        
+def getarea(h):   
+    """Construct face areas of the 3D model as 1d array."""        
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)       
+    # Compute areas of cell faces
+    area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) 
+    area2 = h1.T*mkvc(np.ones((n2+1,1))*h3)
+    area3 = h1.T*mkvc(h2.T*np.ones(n3+1))               
+    area = np.concatenate((mkvc(area1), mkvc(area2), mkvc(area3)), axis=0)
+    
+    return area
+
+def getlength_e(h):   
+    """Construct edge legnths of the 3D model as 1d array."""    
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)       
+    # Compute areas of cell faces
+    l1 = h1.T*mkvc(np.ones((n2+1,1))*np.ones(n3+1))
+    l2 = np.ones((n1+1,1))*mkvc(h2.T*np.ones(n3+1))
+    l3 = np.ones((n1+1,1))*mkvc(np.ones((n2+1,1))*h3)               
+    #l = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3)))
+    l = np.concatenate((mkvc(l1), mkvc(l2), mkvc(l3)), axis=0)
+    
+    return l
+
+def getDivMatrix(h):    
+    """Construct the 3D divergence operator on Faces."""
+           
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)  
+        
+    # Compute areas of cell faces
+    S = getarea(h)            
+        
+    # Compute cell volumes                        
+    V = getvol(h)
+        
+    # Compute divergence operator on faces                       
+    d1 = ddx(n1)
+    d2 = ddx(n2)
+    d3 = ddx(n3)
+    D1 = kron3(speye(n3), speye(n2), d1)
+    D2 = kron3(speye(n3), d2, speye(n1))
+    D3 = kron3(d3, speye(n2), speye(n1))       
+        
+    D = sparse.hstack((D1, D2, D3), format="csr")
+    return sdiag(1/V)*D*sdiag(S)
+
+def getGradMatrix(h):    
+    """Construct the 3D nodal gradient operator."""
+           
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)  
+        
+    # Compute lengths of cell edges
+    L = getlength_e(h)                
+    
+    # Compute divergence operator on faces                       
+    d1 = ddx(n1)
+    d2 = ddx(n2)
+    d3 = ddx(n3)
+    D1 = kron3(speye(n3+1), speye(n2+1), d1)
+    D2 = kron3(speye(n3+1), d2, speye(n1+1))
+    D3 = kron3(d3, speye(n2+1), speye(n1+1))       
+        
+    G = sparse.vstack((D1, D2, D3), format="csr")
+    return sdiag(1/L)*G
+    
+def getCurlMatrix(h):    
+    """Construct the 3D curl operator."""
+           
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)  
+        
+    # Compute lengths of cell edges
+    L = getlength_e(h)                
+
+    # Compute areas of cell faces
+    S = getarea(h)            
+    
+    # Compute divergence operator on faces                       
+    d1 = ddx(n1)
+    d2 = ddx(n2)
+    d3 = ddx(n3)
+    
+    D32 = kron3(d3, speye(n2), speye(n1+1))
+    D23 = kron3(speye(n3), d2, speye(n1+1))
+    D31 = kron3(d3, speye(n2+1), speye(n1))
+    D13 = kron3(speye(n3), speye(n2+1), d1)
+    D21 = kron3(speye(n3+1), d2, speye(n1))            
+    D12 = kron3(speye(n3+1), speye(n2), d1)      
+    
+    O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
+    O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
+    O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
+    
+    C = sparse.vstack((sparse.hstack((O1,-D32, D23)), 
+                       sparse.hstack((D31,O2, -D13)), 
+                       sparse.hstack((-D21,D12, O3))), format="csr") 
+    
+    return sdiag(1/S)*(C*sdiag(L))
+
+def getAverageMatrixF(h):
+    """Construct the 3D averaging operator on cell faces."""    
+
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)  
+    
+    av1 = av(n1)
+    av2 = av(n2)
+    av3 = av(n3)
+    
+    AvF = sparse.hstack(kron3(speye(n3), speye(n2), av1),
+                        kron3(speye(n3), av2, speye(n3)),
+                        kron3(av3, speye(n2), speye(n3)), format="csr")
+    return AvF                      
+        
+def getAverageMatrixE(h):    
+    """Construct the 3D averaging operator on cell edges."""        
+    
+    # Cell sizes in each direction 
+    h1 = h[0]
+    h2 = h[1]
+    h3 = h[2]
+
+    # The number of cell centers in each direction
+    n1 = np.size(h1)
+    n2 = np.size(h2)
+    n3 = np.size(h3)      
+
+    av1 = av(n1)
+    av2 = av(n2)
+    av3 = av(n3)
+    
+    AvE = sparse.hstack(kron3(av3, av2, speye(n1)),
+                        kron3(av3, speye(n2), av1),
+                        kron3(speye(n3), av2, av1), format="csr")
+    return AvE
\ No newline at end of file
diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py
index 204e71cc..241d9067 100644
--- a/SimPEG/sputils.py
+++ b/SimPEG/sputils.py
@@ -3,16 +3,24 @@ from scipy import sparse
 
 def ddx(n):
     """Define 1D derivatives"""
-    return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1)
+    return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1, format="csr")
     
 def sdiag(h):
     """Sparse diagonal matrix"""
-    return sparse.spdiags(h, 0, np.size(h), np.size(h))    
+    return sparse.spdiags(h, 0, np.size(h), np.size(h), format="csr")    
 
 def speye(n):
     """Sparse identity"""
-    return sparse.identity(n)
+    return sparse.identity(n, format="csr")
 
 def kron3(A, B, C):
     """Two kron prods"""
-    return sparse.kron(sparse.kron(A, B), C)  
\ No newline at end of file
+    return sparse.kron(sparse.kron(A, B), C, format="csr")
+
+def spzeros(n1, n2):
+    """spzeros"""
+    return sparse.coo_matrix((n1, n2)).tocsr()
+
+def av(n):
+    """Define 1D averaging operator"""
+    return sparse.spdiags((0.5*np.ones((n+1,1))*[1,1]).T, [0,1], n, n+1, format="csr")
\ No newline at end of file
diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py
new file mode 100644
index 00000000..8175ab2f
--- /dev/null
+++ b/SimPEG/tests/test_curl.py
@@ -0,0 +1,48 @@
+import numpy as np
+
+import sys
+sys.path.append('../')
+from TensorMesh import TensorMesh
+from getDiffop import getCurlMatrix
+
+
+err=0.
+print '>> Test Curl operator'
+for i in range(4):
+    icount=i+1
+    nc = 2**icount
+    # Define the mesh    
+    h1 = np.ones((1,nc))/nc
+    h2 = np.ones((1,nc))/nc
+    h3 = np.ones((1,nc))/nc
+    h = [h1, h2, h3]
+    x0 = np.zeros((3, 1))
+    M = TensorMesh(h, x0)
+    #n = M.plotGrid()
+
+    # Generate DIV matrix
+    CURL = getCurlMatrix(h)
+    #Test function
+    fun = lambda x: np.cos(x)  # i (cos(y)) + j (cos(z)) + k (cos(x))
+    sol = lambda x: np.sin(x)  # i (sin(z)) + j (sin(x)) + k (sin(y))
+    
+    Ex = fun(M.gridEx[:,1])
+    Ey = fun(M.gridEy[:,2])
+    Ez = fun(M.gridEz[:,0])
+    E = np.concatenate((Ex,Ey,Ez))     
+
+    Fx = sol(M.gridFx[:,2])
+    Fy = sol(M.gridFy[:,0])
+    Fz = sol(M.gridFz[:,1])
+    curlE_anal = np.concatenate((Fx,Fy,Fz))     
+
+    curlE = CURL*E    
+    err = np.linalg.norm((curlE-curlE_anal), np.inf)
+
+    if icount == 1:
+        print 'h       |   inf norm   | error ratio'        
+        print '---------------------------------------'                
+        print '%6.4f  |  %8.2e  |'% (h1[0,0], err)
+    else:
+        print '%6.4f  |  %8.2e  |  %6.4f' % (h1[0,0], err, err_old/err)
+    err_old = err            
\ No newline at end of file
diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py
index bc0c1a9d..3ab8f6e2 100644
--- a/SimPEG/tests/test_div.py
+++ b/SimPEG/tests/test_div.py
@@ -5,16 +5,18 @@ sys.path.append('../')
 from TensorMesh import TensorMesh
 from getDIV import getDivMatrix, getarea, getvol
 
-# Define the mesh
+
 err=0.
+print '>> Test face Divergence operator'
 for i in range(4):
-    icount=i+1;
-    nc = 2*icount;
-    h1 = np.pi/nc*np.ones((1,nc))
-    h2 = np.pi/nc*np.ones((1,nc))
-    h3 = np.pi/nc*np.ones((1,nc))
+    icount=i+1
+    nc = 2**icount
+    # Define the mesh        
+    h1 = np.ones((1,nc))/nc
+    h2 = np.ones((1,nc))/nc
+    h3 = np.ones((1,nc))/nc
     h = [h1, h2, h3]
-    x0 = -np.pi/2*np.ones((3, 1))
+    x0 = np.zeros((3, 1))
     M = TensorMesh(h, x0)
     #n = M.plotGrid()
 
@@ -34,12 +36,13 @@ for i in range(4):
      
     area = getarea(h)
     vol = getvol(h)
-    err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
+    #err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2)
+    err = np.linalg.norm((divF-divF_anal), np.inf)
+
     if icount == 1:
-        err1 = err
-        print 'h       |   2 norm   | error ratio'        
+        print 'h       |   inf norm   | error ratio'        
         print '---------------------------------------'                
         print '%6.4f  |  %8.2e  |'% (h1[0,0], err)
     else:
-        print '%6.4f  |  %8.2e  |  %6.4f' % (h1[0,0], err, err1/err)
-        
+        print '%6.4f  |  %8.2e  |  %6.4f' % (h1[0,0], err, err_old/err)
+    err_old = err    
diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py
new file mode 100644
index 00000000..a433c28d
--- /dev/null
+++ b/SimPEG/tests/test_grad.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+import sys
+sys.path.append('../')
+from TensorMesh import TensorMesh
+from getDiffop import getGradMatrix
+
+
+err=0.
+print '>> Test nodal Gradient operator'    
+for i in range(4):
+    icount=i+1
+    nc = 2**icount
+    # Define the mesh    
+    h1 = np.ones((1,nc))/nc
+    h2 = np.ones((1,nc))/nc
+    h3 = np.ones((1,nc))/nc
+    h = [h1, h2, h3]
+    x0 = np.zeros((3, 1))
+    M = TensorMesh(h, x0)
+    #n = M.plotGrid()
+
+    # Generate DIV matrix
+    GRAD = getGradMatrix(h)
+    #Test function
+    fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z))    
+    sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z))
+    
+    phi = fun(M.gridN[:,0], M.gridN[:,1], M.gridN[:,2])
+    gradE = GRAD*phi
+
+    Ex = sol(M.gridEx[:,0])
+    Ey = sol(M.gridEy[:,1])
+    Ez = sol(M.gridEz[:,2])
+
+    gradE_anal = np.concatenate((Ex,Ey,Ez))     
+    err = np.linalg.norm((gradE-gradE_anal), np.inf)
+
+    if icount == 1:
+        print 'h       |   inf norm   | error ratio'        
+        print '---------------------------------------'                
+        print '%6.4f  |  %8.2e  |'% (h1[0,0], err)
+    else:
+        print '%6.4f  |  %8.2e  |  %6.4f' % (h1[0,0], err, err_old/err)
+    err_old = err    
+