diff --git a/SimPEG/Examples/FWI/seismic.py b/SimPEG/Examples/FWI/seismic.py
index e76bd42d..e96ba43d 100644
--- a/SimPEG/Examples/FWI/seismic.py
+++ b/SimPEG/Examples/FWI/seismic.py
@@ -1,5 +1,8 @@
-import sys
-sys.path.append('/Users/haber/dropbox/simpegMaster/')
+try:
+    import sys
+    sys.path.append('/Users/haber/dropbox/simpegMaster/')
+except Exception, e:
+    pass
 
 import SimPEG
 import scipy.sparse as sp
@@ -16,12 +19,12 @@ def HelmholtzSol(model,mesh,w,mbc,q,P):
 
     k = np.sqrt(w**2 * mbc)
 
-    D1 = utils.sdiag(1./h1) * utils.ddx(mesh.nCx)
-    D2 = utils.sdiag(1./h2) * utils.ddx(mesh.nCy)
-    
-    L1 = D1.T*D1; L2 = D2.T*D2;
-    print L1.todense()
-    
+    # D1 = utils.sdiag(1./h1) * utils.ddx(mesh.nCx)
+    # D2 = utils.sdiag(1./h2) * utils.ddx(mesh.nCy)
+
+    # L1 = - D1.T*D1
+    # L2 = - D2.T*D2
+
     Av = mesh.aveN2CC
 
     B1 = utils.spzeros(n1+1,n1+1);
@@ -33,21 +36,23 @@ def HelmholtzSol(model,mesh,w,mbc,q,P):
     B2[0,0] = 2*1j*k*h2[0]; B2[-1,-1] = 2*1j*k*h2[-1]
 
     # generate the 2D Laplacian
-    L = sp.kron(sp.identity(n2+1),L1+B1) + sp.kron(L2+B2,sp.identity(n1+1)) 
-
+    # L = sp.kron(sp.identity(n2+1),L1) + sp.kron(L2,sp.identity(n1+1))
+    L = mesh.nodalLaplacian
+    B = sp.kron(sp.identity(n2+1),B1) + sp.kron(B2,sp.identity(n1+1))
+    L = L+B
     #plt.spy(L)
     #plt.show()
     # Generate the Mass matrix
     M = utils.sdiag(Av.T*utils.mkvc(model))
-    A = L - w**2 * M
+    A = - L - w**2 * M
 
     mesh.ForModMat = A
     u = sp.linalg.spsolve(A,q);
     d = P*u
     return u, d
-    
+
 def HelmholtzJmatVec(v,model,u,mesh,w,mbc,P):
- 
+
     Cm = -w**2 * utils.sdiag(u)*mesh.aveN2CC.T
     Cu = mesh.ForModMat
     Cmv = Cm*v;
@@ -55,20 +60,20 @@ def HelmholtzJmatVec(v,model,u,mesh,w,mbc,P):
     return -P*lam
 
 def HelmholtzJTmatVec(v,model,u,mesh,w,mbc,P):
- 
+
     Cm = -w**2 * utils.sdiag(u)*mesh.aveN2CC.T
     Cu = mesh.ForModMat
     Pv = -P.T*v
     z = sp.linalg.spsolve(Cu.T,Pv);
     return Cm.T*z
-     
-    
+
+
 if __name__ == '__main__':
     # odel,mesh,w,mbc,q
     n1 = 128; n2 = 128
     h1 = np.ones(n1); h2 = np.ones(n2);
     P  = sp.identity((n1+1)*(n2+1))
-    
+
     mesh  = TensorMesh([h1,h2])
     model = np.ones(mesh.nC)
     w     = 1
@@ -76,21 +81,21 @@ if __name__ == '__main__':
     q     = np.zeros((mesh.nNx,mesh.nNy))
     q[n1/2,n2/2] = 1.0
     q     = q.reshape(mesh.nN,order = 'F')
-    
+
     u, d = HelmholtzSol(model,mesh,w,mbc,q,P)
     u = u.reshape((mesh.nCx+1,mesh.nCy+1),order = 'F')
-    
+
     plt.imshow(u.real)
     plt.show()
-    
+
     dm = np.random.rand(mesh.nC)*1e-1+2
     u1, d1 = HelmholtzSol(model+dm,mesh,w,mbc,q,P)
 
     dd = HelmholtzJmatVec(dm,model,u,mesh,w,mbc,P)
-    
+
     print np.linalg.norm(d1-d)
     print np.linalg.norm(d1-d-dd)
-    
-    
-    
-    
\ No newline at end of file
+
+
+
+
diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py
index b767ef31..fa3f42e0 100644
--- a/SimPEG/mesh/DiffOperators.py
+++ b/SimPEG/mesh/DiffOperators.py
@@ -97,6 +97,38 @@ class DiffOperators(object):
     _nodalGrad = None
     nodalGrad = property(**nodalGrad())
 
+    def nodalLaplacian():
+        doc = "Construct laplacian operator (nodes to edges)."
+
+        def fget(self):
+            if(self._nodalLaplacian is None):
+                print 'Warning: Laplacian has not been tested rigorously.'
+                # The number of cell centers in each direction
+                n = self.n
+                # Compute divergence operator on faces
+                if(self.dim == 1):
+                    D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
+                    L  = - D1.T*D1
+                elif(self.dim == 2):
+                    D1 = sdiag(1./self.hx) * ddx(n[0])
+                    D2 = sdiag(1./self.hy) * ddx(n[1])
+                    L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
+                    L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
+                    L  = L1 + L2
+                elif(self.dim == 3):
+                    D1 = sdiag(1./self.hx) * ddx(n[0])
+                    D2 = sdiag(1./self.hy) * ddx(n[1])
+                    D3 = sdiag(1./self.hz) * ddx(n[2])
+                    L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
+                    L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
+                    L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
+                    L  = L1 + L2 + L3
+                self._nodalLaplacian = L
+            return self._nodalLaplacian
+        return locals()
+    _nodalLaplacian = None
+    nodalLaplacian = property(**nodalLaplacian())
+
     def setCellGradBC(self, BC):
         """
         Function that sets the boundary conditions for cell-centred derivative operators.
@@ -172,7 +204,6 @@ class DiffOperators(object):
         return locals()
     cellGradx = property(**cellGradx())
 
-
     def cellGrady():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
@@ -193,8 +224,6 @@ class DiffOperators(object):
         return locals()
     cellGrady = property(**cellGrady())
 
-
-
     def cellGradz():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
@@ -212,7 +241,6 @@ class DiffOperators(object):
         return locals()
     cellGradz = property(**cellGradz())
 
-
     def edgeCurl():
         doc = "Construct the 3D curl operator."
 
@@ -255,6 +283,8 @@ class DiffOperators(object):
     _edgeCurl = None
     edgeCurl = property(**edgeCurl())
 
+    # --------------- Averaging ---------------------
+
     def aveF2CC():
         doc = "Construct the averaging operator on cell faces to cell centers."
 
@@ -356,6 +386,8 @@ class DiffOperators(object):
     _aveN2F = None
     aveN2F = property(**aveN2F())
 
+    # --------------- Methods ---------------------
+
     def getMass(self, materialProp=None, loc='e'):
         """ Produces mass matricies.